Entanglement Dynamics in Dispersive Optomechanics: Non-Classicality and Revival

We study entanglement dynamics in a dispersive optomechanical system consisting of two optical modes and a mechanical oscillator inside an optical cavity. The two optical modes interact with the mechanical oscillator, but not directly with each other. The appearance of optical entanglement witnesses non-classicality of the oscillator. We study the dependence of the entanglement dynamics with the optomechanical coupling, the mean cavity photon numbers and temperature. We propose an experimental realization with ultracold atomic ensembles.

Introduction -Entanglement is one of the most striking phenomena of quantum theory [1]. Generating, manipulating and measuring entanglement in systems with many constituents and with a large number of degrees of freedom is one of the challenges of Quantum Information and Metrology [2], and an interesting frontier in Fundamental Physics [3]. In particular, entangling massive objects could open the way to interesting tests of quantum theory [4,5] and experiments aimed at probing gravitational effects of quantum mechanical matter [6][7][8][9][10][11][12].
It is well known that entanglement of massive objects can be realized in quantum cavity optomechanical experiments [13]. For instance, a cavity with a moving end mirror can be used to generate entangled "cat states" of both light [14,15] and matter [16][17][18][19]. Similar systems have also been proposed as an effective nonlinear medium [20][21][22] and squeezing [23][24][25] as well as optical entanglement [26] have been experimentally demonstrated in a variety of set-ups such as cavity cold atomic ensembles [27], dispersive dielectric membranes [28] and silicon micro-resonators [29].
Certifying quantumness of optomechanical systems, however, is a far-from-obvious task. Relations among entanglement and non-classicality measures of quantum states can be used to probe the quantum nature of innaccessible objects such as a harmonic oscillator in an optical cavity [30]. Recurrence of optical squeezing in a cavity with a moving mirror has also been proposed as a witness of non-classicality [31] and it has been shown that when two subsystems locally interact with a third one, but not directly with each other, the appearance of entanglement among those subsystems is sufficient to prove non-classicality of the third party [32]. Building on some of these ideas the present work studies the intricate entanglement dynamics of a dispersive optomechanical system, how to use that dynamics to probe the quantum nature of the oscillator when only optical degrees of freedom are accessible and how to optimize the generated optical entanglement by careful choice of the optomechanical coupling and the number of photons in an experiment.
Considering as possible implementations levitated op- ing a "mirror-in-the-middle" under a harmonic potential. No photon transfer from one cavity to the other is allowed. (b) Schematics of a particle trapped by an optical tweezer coupled to two modes of a cavity. The particle can be considered as a Silica nano-sphere or a cloud of ultracold atoms. When the levitated object is properly positioned, the Hamiltonian describing both systems acquires the same form.
tomechanical systems such as Silica nano-spheres or cold atomic ensembles and a "two-sided" cavity with a moving mirror in the middle, we map how entanglement and entropy appear and evolve among the various optical and mechanical subsystems for different optomechanical coupling strengths and optical field intensities. The appearance of mechanically induced optical entanglement and its subsequent death and revival are generic in these systems, implying quantization of the center-of-mass motion of the moving object. We also point that under certain circumstances entanglement seems to "flow" through different subsystems, and such dynamics can be used to infer non-classicality and entanglement among different components of the system. We consider examples of both non-Gaussian and Gaussian initial quantum states, for which we study the dynamics of concurrence and logarithmic negativity, respectively. An experimental realization using levitated cold atomic ensembles is proposed.
Hamiltonian description -The system we are primarily interested in is shown in Figure 1(a): two optical cavities, of lengths L a and L b , are populated by modes of frequencies ω a and ω b and share a common perfect movable mirror of mass m subject to a harmonic potential of frequency ω m . We refer to this as the "mirror-in-themiddle" configuration. In this system the optical modes never interact directly, except via the dispersive coupling due to the presence of the mechanical mode. Since we are interested in studying entanglement dynamics in optomechanics, we shall assume the cavities can be initialized in particular states and the laser driving-term can be turned off during the course of the experiment. It is also assumed that optical losses are negligible during the time of the experiment and a discussion of the conditions under which this is true and the experimental feasibility is addressed in the experimental proposal.
The Hamiltonian of the system reads [13] where g 0,i = ω i x zpf /L i are the optomechanical couplings, with x zpf = /2mω m the zero point fluctuation of the mirror andâ,b,ĉ (â † ,b † ,ĉ † ) are the annihilation (creation) operators of each optical and mechanical modes, denoted by A, B and C, respectively. Such Hamiltonian can also be implemented using a cavity with a levitated nano-particle [33][34][35] or an ultracold atom cloud [27,36,37] properly positioned along the cavity axis. This is illustrated in Figure 1(b); see the Supplemental Material for details on implementations.
Assuming equal frequencies for the optical modes ω a = ω b , and approximately equal cavity lengths L a ∼ L b we may simplify the notation and directly write g 0 ≡ g 0,a ∼ g 0,b . The unitary evolution operator resulting from equation (1) becomeŝ where we define the dimensionless optomechanical coupling k = g 0 /ω m , r i = ω i /ω m , the scaled time ω m t → t, and the functions η(t) = 1 − e −it and B(t) = −k 2 (t − sin t). Note the evolution operator is comprised of a Kerrlike term, responsible for an effective optical non-linearity [38], as well as an optically-driven displacement operator acting on the mechanical mode.
We expect that a generic separable state evolves into an entangled state by the unitary evolution above. This is not always the case, as can be seen by considering the energy eigenstates of the system where {|n A , m B , C } denotes the number basis and D C (α) the displacement operator acting on the mechanical oscillator, mode C, by a displacement α ∈ C. See the Supplemental Material for a derivation of equation (2) and of the energy eigenstates above. Non-Gaussian states -Consider that the cavities in Figure 1(a) are initially populated by the non-Gaussian state These states can be prepared in an approximate way using photon pair sources and displacement-based detection [39,40] or non-linear light-matter interactions in cavity quantum electrodynamics [41,42]. The "mirrorin-the-middle" is taken to be in the ground state for simplicity; in the next section we shall consider the effects of a finite temperature oscillator. Notice the initial state is separable and hence the appearance of entanglement between modes A and B would evidence the non-classical nature of mode C [32]. Time evolution of (4) in the interaction picture is explicitly given by where ξ(t) = e it η(t). It is expected that the evolved state (5) exhibits entanglement between modes A and B. This can be promptly seen by noticing that coherent states are non-orthogonal and, in the limit of small coupling k, the state assumes the form |Ψ(t) |ϕ AB ⊗ |ϕ C , where and |ϕ C |0 . For times t such that |B(t)| π/2 + 2πn, n ∈ N, the state |ϕ AB becomes maximally entangled. The origin of this entanglement can be heuristically explained by a simple argument: the ground state of the mechanical oscillator is a Gaussian wave packet in the position basis. Each possible position adds-up coherently introducing correlations in the lengths of the left and right cavities in Figure 1(a). This imprints correlations in the phases of the corresponding electromagnetic fields in modes A and B giving rise to entanglement. As long as the dimensionless optomechanical coupling k due to radiation pressure on the middle mirror is sufficiently small, mode C will be approximately unperturbed and therefore, to a good approximation, disentangled from AB. On the other hand, as the coupling strength increases mode C can become significantly entangled with modes A and B; the entanglement of subsystems as a function of k will be addressed in the next section.
To quantitatively evaluate the entanglement in (5) we calculate the three-partite density matrix ρ ABC = |Ψ(t) Ψ(t)|, from which we obtain the reduced state ρ AB = Tr C (ρ ABC ): with C(t) = iB(t) − k 2 |η(t)| 2 /2. Notice that for small values of k the mirror is "weakly entangled" with modes A and B and some of the off-diagonal terms of the reduced density matrix acquire exponentials that alternate between periods of decay and periods of growth. This can be seen as an example of a weak form of decoherence and "non-Markovian" evolution for the partitions of the whole system, in which information about the optical modes leak into correlations with the mirror and is later retrieved. The mirror introduces a "memory" in the system [43]. Non-Markovianity springs from the fact that the mirror is part of the system under study and hence its degrees of freedom are under control.
Since the state ρ AB (t) is restricted to the subspace spanned by {|0 , |1 }, we can use concurrence as a measure of entanglement. As an example, plots of the concurrence C AB (t) and von Neumann entropy S AB (t) of ρ AB are shown in Figure 2 for coupling value k = 0.5. The optical modes exhibit positive concurrence and hence entanglement as a function of time. Moreover, the system exhibits sudden death and birth of entanglement. It is also possible to see that the entropy, which is initially zero, oscillates as a function of time. This is another indication of the non-Markovian nature of the system. A non-zero entropy of AB signals entanglement among the three-partite system ABC. Moreover, note that the maxima of concurrence (entanglement of AB) coincide with the minima of entropy (entanglement of ABC). This suggests that entanglement "flows" (during the limited period of its lifetime) among different partitions of the system. Gaussian states -We now consider a scenario in which initially the optical modes are populated by monochromatic coherent states and the moving object (sphere, cloud of atoms or mirror) is in a thermal state at temperature T where Z = n e − n ωc k B T is the thermal partition function. Note ρ(0) is a Gaussian separable state. Since the dynamical evolution (2) preserves Gaussianity [44], we can conveniently use the Continuous Variable formalism to study the dynamics of entanglement and entropy of the subsystems through their time-dependent covariance matrices [44,45] from which we can calculate exactly the logarithmic negativity, an entanglement monotone [46,47], and the Von Neumann entropy. See Supplemental Material for an overview of the Continuous Variable formalism and details of the calculations. Figure 3(a) shows an example of the logarithmic negativities and entropies for the different partitions of the system: AB (opto-opto), BC (opto-mechanical) and AC (opto-mechanical). Entanglement and entropy exhibit oscillations (death and revivals) of which one full period is shown. Note the initial decrease and subsequent increase in the entropies AC and BC, a characteristic of the non-Markovian nature of the system. As in the non-Gaussian case, the appearance of entanglement between optical modes, given the initially separable state, evidences the non-classicality of the mechanical mode [32].
We observe that the local maxima of logarithmic negativity for subsystems AC and BC coincide with local minima of the respective von Neuman entropies. At the same time, this reduction in entropy is accompanied by an increase in the entropy of modes AB. This suggests that as the entropy "leaks" from the mechanical oscillator into the optical modes, entanglement between modes AC and BC "builds up" and is later transferred to the optical modes. As in the non-Gaussian example, there seems to be evidence for the idea that entanglement flows among the various interacting subsystems. To firmly establish this idea, further investigation is necessary.
We note that the entanglement and entropy dynamics is strongly affected by the dimensionless optomechanical coupling k in a non-trivial way. This can be seen from the correlators of the system, e.g. equation (9), and from Figure 3(b) showing the maximum of the logarithmic negativity for the optical subsystem within one period τ as a function of the coupling k. See the Supplemental Material for further details.
When k < 1/ √ 2, the "low coupling" regime, the overall periodicity of the logarithmic negativities and von Neumann entropies is τ = π/ω m k 2 and observation of entanglement revivals are only possible when π/ω m k 2 κ −1 , where κ −1 is the inverse cavity linewidth, or the approximate photon lifetime in the cavity. This translates into the so-called photon blockade condition g 2 0 /(ω m κ) 1 [36,48]. Furthermore, from Figure 3(b), we observe that max(E AB ) is mostly constant as a function of k. On the other hand, when k ≥ 1/ √ 2, "high coupling" regime, the periodicity is τ = 2π/ω m and observation of entanglement dynamics is conditioned on satisfying 2π/ω m κ −1 . Moreover, max(E AB ) becomes very sensitive to small changes in k, presenting local minima at 2k 2 = N , where N ∈ N * . Knowledge of this behavior can be used when designing an experiment for tailoring the entanglement properties of the system as desired. For instance, if optical entanglement is to be maximized, increasing coupling may not be the best strategy.
To characterise the effect of the oscillator's thermal fluctuations in the the optical entanglement we show a plot of max(E AB ) as a function of temperature in Figure  3(c). The entanglement persists well above the microkelvin temperatures reported in current optomechanical experiments [27,49].
The mean number of photons in the cavity also plays an important role in optical entanglement generation. Figure 3(d) shows a surface plot of max(E AB ) as a function of the coherent state amplitudes α and β, taken to be real numbers for simplicity. We observe an increase in the maximal optical entanglement with mean photon numbers. There is, however, a trade-off: as the number of photons increase, the width of the peaks in E AB (t)and consequently the time during which entanglement is non-vanishing -decreases significantly, a feature which can also be qualitatively understood from the form of the correlators such as the one in (9). Moreover, it is also possible to see that entanglement is optimized when the energy is evenly distributed in the two optical modes, which happens when α = β. Interestingly, as α and β increase the system becomes less sensitive to discrepancies between the mean photon numbers in each optical mode.
Experimental proposal -With increasing advances in the field of quantum cavity optomechanics [49][50][51] experiments in the high coupling and long coherence time regimes are expected. However, most current optomechanical systems have small coupling values and observing entanglement dynamics as described in the present work remains challenging. One notable exception and a promising candidate set-up is optomechanics with ultracold atomic ensembles where a coherent cloud of atoms is trapped within an optical cavity and the collective center of mass coordinate effectively behaves as a quantum mechanical oscillator. Couplings as high as k ≈ 10 have been reported in such ultracold experiments [27], and the system allows wide tunability of the relevant parameters. To observe entanglement dynamics the lifetime of a photon inside the cavity 1/κ needs to be longer than the 'entangling time', the time at which optical entanglement reaches its first local maximum. Table I presents values for the relevant experimental parameters adapted from [27]. With these numbers the resulting cavity photon lifetime is 15.7 × 10 −6 s while the entangling period is 10.5 × 10 −6 s. The corresponding negativity was found to be max(E AB ) ≈ 0.3. This suggests that in these systems the observation of dispersive mechanically-induced optical entanglement could be within reach.
Conclusion -In this Letter we have studied the entanglement and entropy dynamics of a "mirror-in-themiddle" optomechanical system. Implementations using levitated particles have been briefly discussed. We have seen that an initially separable state can evolve to an entangled wavefunction, exhibiting birth and death of entanglement and entropy both for non-Gaussian and Gaussian states. The appearance of entanglement in this setting evidences the non-classical nature of the mechanical oscillator.
The entanglement dynamics is strongly influenced by the system parameters, notably the optomechanical coupling k = g 0 /ω m and the mean number of photons in the experiment. We have shown the existence of two distinct regimes depending on whether k < 1/ √ 2 or k ≥ 1/ √ 2. In addition, we have observed that optical entanglement is maximized when the energy is evenly distributed in the optical modes and that it persists well above microkelvin temperatures of the mechanical mode. This is valuable information when designing an experiment. Optomechanics with ultracold atomic ensembles presents an interesting candidate for implementing the studied entanglement dynamics.
Although a promising candidate, the dispersive Hamiltonian is not the only available platform to untangle the dynamics of entanglement and information flow in optomechanical systems. Exploring alternatives such as coherent scattering [52][53][54] might prove to be a very fruitful approach to observe entanglement and non-classicality in experimental optomechanical systems. Acknowledgments

I. IMPLEMENTATION USING LEVITATED NANOPARTICLES
Consider the system shown in Figure S1: a nanoparticle of radius r, mass m and refractive index n p is trapped in an harmonic trap of frequencies ω j=x,y,z created by an optical tweezer inside a cavity populated by two optical modes of frequencies ω a and ω b and annihilation/creation operatorsâ/â † andb/b † , respectively. The presence of the particle causes a position dependent shift on the cavity's resonance frequencies, so that the Hamiltonian of this system becomes [35] where x 0 is the position of the center of the trap, x is the particle's displacement, U 0,i = ωâ i α/2 0 V i is the frequency shift when the particle is at an intensity maximum of the cavity, α = 4π 0 r 3 (n 2 p − 1)/(n 2 p + 2) is the polarizability of the particle, V i and k i are the volume and the wavenumber of mode i, respectively, andĉ j /ĉ † j is the phonon annihilation/creation operators along axis j = x, y, z. We will not consider driving terms since we are interested in the dynamics of an isolated cavity without the influence of external systems such as a driving laser, as explained in the main text. The interaction terms may yield linear couplings between the optical modes and the sphere if x 0 , k a and k b are properly chosen. Expanding sin 2 k i (x 0 + x) around x 0 gives Now, consider the particular case in which the frequencies of the optical modes are two consecutive resonance frequencies, and let L = 2n(λ a /2) = (2n + 1)(λ b /2) without loss of generality. Then, if the sphere is placed near the center of the cavity at x 0 = L/2 + λ a /8 ≈ L/2 + λ b /8, we have and Substituting these approximations in equation (S1) and disregarding the motion along the y and z axes, we get where ω a = ω a − U 0,i /2 and g 0 = U 0,a k a x ZP F ≈ U 0,b k b x ZP F is the intended linear coupling, with x ZP F = /2mω x the zero-point fluctuation. This Hamiltonian is formally equivalent to that of a double-sided cavity with a moving "mirror-in-the-middle" analysed in the main text, with the condition g 0,a ≈ g 0,b being met due to the use of consecutive resonance frequencies. As a final remark, note that the modes A and B are not used for cooling of the particle, which is necessary for observing optical entanglement in the system [32]. Cooling of the center-of-mass motion of the particle can be addressed in a number of ways, such as a feedback acting on the trapping laser [55,56], dispersive coupling with a third optical mode [57] or coherent scattering between the trapping beam and the optical cavity [50,51,58].

II. IMPLEMENTATION USING ULTRACOLD ATOMIC ENSEMBLES
In ultracold atom optomechanical experiments, an atomic ensemble is trapped inside an optical cavity. Collective center-of-mass motion of the atoms alters the cavity resonance frequency. It is similar to the dispersive optomechanical experiments with levitated spheres described in the previous section, with the cloud of atoms playing the role of levitated nanoparticle.
The optomechanical coupling between an ultracold atomic cloud and a cavity optical mode with wavelength λ a is [27] where k a is the wavenumber, N is the number of atoms, ∆ ca is the atom-cavity detuning, m is the mass of a single atom, ω m is the mechanical frequency and α 0 = d 2 ω c /2 0 V c is the atom-single photon coupling rate, with d the dipole moment for the transition between the relevant atomic levels and V c the cavity volume. Reported values from Ref. [27] for these quantities are shown in Table S1. As in the previous case of a levitated nanoparticle, we chose sin (2k a z 0 ) = 1. The wavelength λ b should be chosen so that sin (2k b z 0 ) = −1, as to provide g 0,a = g 0,b . As discussed in the main text, for the entanglement dynamics experiment to be feasible, the photon-lifetime τ p must be greater than the entanglement period τ e . One of the main constraints in fulfilling this condition is set by the Finesse of the cavity. Therefore, we look for the minimum value of τ e /τ p by varying the cavity length L, the number of atoms N and the mechanical frequency ω m around the values in Table S1, and then calculate the minimum Finesse necessary to make τ e /τ p < 1 given the optimal values of L, N and ω m .
In doing these calculations, it is necessary to account for the changes in the mode volume, given by V c = πw 2 a L, where is the mode's waist, and the changes in the cavity linewidth with ν F SR = c/2L the cavity's free spectral range. Finally, it is important to make sure that k = g 0 /ω m = N/2, N ∈ N * , since the maximum of logarithmic negativity is significantly smaller for these coupling values. We find that for L = 783 µm, N = 5.43 × 10 5 and ω m = 2π × 95 kHz the dimensionless optomechanical coupling is k = 0.743 and the entanglement period to photon lifetime ratio is τ e /τ p = 3.46. The Finesse should then be increased to 2.01 × 10 6 , so that τ e /τ p 1 and entanglement becomes measurable. In the main text, the proposed Finesse is about 1.5 times larger, so that τ e /τ p = 0.669. Another parameter that could be varied are the radii of the cavity mirrors. Considering 1 cm, 2.5 cm, 5 cm and 10 cm as possible radii, we find the values presented in Table S2. As we can see, the Finesse constraint can be relaxed provided that a smaller radius is used. Overall, the necessary values for N , ω m and F differ by less than one order of magnitude from reported values in the literature.

III. THE UNITARY EVOLUTION
The "mirror-in-the-middle" Hamiltonian can be written aŝ whereâ,b andĉ denote annihilation operators for the optical mode in the left cavity (A), right cavity (B) and the mechanical oscillator (C), respectively. In what follows we will work with the re-scaled time ω m t, denoted henceforth by t, and introduce the dimensionless variables r a = ω a /ω m , r b = ω b /ω m and k = g 0 /ω m . Define the unitary operator [16]Ê The operatorÊ(k) commutes withâ †â andb †b , but not withĉ, Using equation (S10) and its adjoint, we havê Considering the number basis {|n, m, }, we see by the equation above thatÊ(k) |n, m, are the energy eigenstates of the system. Those are of the form:Ê where we denote the displacement operator of the mechanical oscillator byD C (κ) = exp κĉ † − κ * ĉ , where κ is a complex number. The energies corresponding to the eigenstates above are E n,m, = ω m + ω a n + ω b m − ω m k 2 (n − m) 2 . (S14) By exponentiation of equation (S12), the unitary evolution operator is found to bê Next, we recall that e iĉ †ĉ tĉ † e −iĉ †ĉ t = e itĉ † and e iĉ †ĉ tĉ e −iĉ †ĉ t = e −itĉ , from which we derive the following identity whereÂ is an operator that commutes with bothĉ andĉ † . LettingÂ = k(â †â −b †b ) in the identity above, we conclude that We note the commutator with which it becomes straightforward to compute Since this commutes with both (â †â −b †b )(ĉ † e it −ĉe −it ) and (â †â −b †b )(ĉ † −ĉ), equation (S18) can be simplified using the particular form of the Baker-Campbell-Haursdoff formula valid when operatorsÂ andB commute with [Â,B]. Thus, the expression holds for the unitary evolution operatorÛ (t), where η(t) = 1 − e −it and B(t) = −k 2 (t − sin(t)).
In the interaction picture, we evolve states according tô The Heisenberg evolution of the annihilation operatorsâ,b andĉ can then be derived, where ξ(t) = e it η(t) = e it − 1 = −η(t) * . In deriving these expressions we made use of the relations applicable wheneverÂ commutes with bothâ andâ † . For instance, equation (S26) is derived as which leads to the desired result after noting that −η(t) * e −it = η(t).

IV. QUANTIFYING ENTANGLEMENT AND ENTROPY FOR GAUSSIAN STATES
To quantify entanglement for all bipartite subsystems we employed the logarithmic negativity, as it is an entanglement monotone [44] and efficiently computable for Gaussian states. The initial state of the system is assumed to be: where the cavity modes are in a coherent state and the mirror is in a thermal state ρ th ∝ exp − ωmĉ †ĉ k B T . All modes and bipartite subsystems are initially in Gaussian states; moreover, the unitary evolution operator (S22) preserves Gaussianity, so both quantities above can be computed from the covariance matrix of the subsystem [44,45]. Consideringâ 1 andâ 2 to be the annihilation operators of any of the bipartite subsystems, we define their corresponding quadrature field operators asQ i =â i +â † i ,P i = i â † i −â i and the formal vectorx = (Q 1 ,P 1 ,Q 2 ,P 2 ) such that the covariance matrix for this subsystem is defined as: We note that the time-dependent expectation values can be calculated from the initial state ρ 0 using the Heisenberg picture, see equations (S24) -(S26). It is useful to rewrite the covariance matrix in terms of its 2 × 2 submatrices V 1 , V 2 , V 3 : To compute the von Neumann entropy, we first define ∆ ≡ det(V 1 )+det(V 2 )+2 det(V 3 ) and calculate the symplectic eigenvalues of V , given by [45] , the von Neumann entropy S is then given by [44] The logarithmic negativity can be computed as whereν − is the smallest of the symplectic eigenvalues of the covariance matrix calculated with respect to the partial transpose of the reduced density matrix of the subsystem under discussion [45]. In practice, the effect of the partial transposition mentioned above is that the symplectic eigenvalues ofṼ differ from equations (S33) by a sign change in det(V 3 ),ν where σ = det(V 1 ) + det(V 2 ) − 2 det(V 3 ).
To compute the covariance matrices, we proceeded by calculating the averages of all quadratic expressions in a,b,ĉ and their adjoints. By applications of identities already mentioned in the previous section, together with the displacement relationĉD C (κ)(κ) =D C (κ)(κ)(ĉ + κ), we find the following expressions forâ(t) 2 The operatorsâ †â andb †b commute with H, so are constant. All the remaining quadratic expressions follow from the equations above by commutation or by considering the adjoint equations -we then compute the average values of those expressions with respect to the state ρ 0 of the system, and finally the covariance matrices for the three bipartite subsystems.
By Gaussian integrals The following relations are also useful, as well as the analogous expressions forb, with β replacing α. The first one follows from Taylor expanding the exponential, while the second follows from the first by commutingâ † to the left with the adjoint of equation (S27). The identity bellow is trivial to check: Applying the identities above to the expressions for the quadratic terms in the annihilation and creation operators, we compute the matrix elements of all three covariance matrices and therefore quantify entanglement for all bipartite subsystems by calculating the von Neumann entropies and the logarithmic negativities.
Considered the order of magnitude of the experimental parameters used in the main text, ω a = ω b ∼ 10 15 , ω m ∼ 10 5 Hz, the free evolution terms oscillate very rapidly compared to the others and can be approximated on average as constant. Consequently, the general periodicity of the entanglement and entropy are dictated solely by the periodicity of the thermal and coherent terms, given by τ 1 = 2π/ω m and τ 2 = π/ω m k 2 , respectively.
When τ 1 τ 2 ⇒ k 1/ √ 2, the "low coupling" regime, the coherent terms dominate the overall form of the logarithmic negativity and entropy, while the thermal terms describe fine oscillations. Therefore the periodicity of entanglement and entropy is τ = π/ω m k 2 . Numerically, we have observed that the maximum of the logarithmic negativity within one period τ is mostly constant, but since the period of entanglement grows with 1/k 2 as k goes to 0, it becomes potentially much larger than the lifetime of a photon inside the optical cavities and becomes impossible to observe.
When τ 1 τ 2 ⇒ k 1/ √ 2, the "high coupling" regime, the thermal terms dominate the overall form of the logarithmic negativity and entropy, while the coherent terms describe fine oscillations. Therefore the periodicity of entanglement and entropy is τ = 2π/ω m . Numerically, we have observed that the maximum of the logarithmic negativity within one period τ is very sensitive to small changes of k, presenting local minima many orders of magnitude smaller than at its surrounding when τ 1 = N τ 2 ⇒ 2k 2 = N , N ∈ N * . As k grows much bigger than 1, any small change in k can easily throw the logarithmic negativity into a much smaller value than its surroundings.