Quantum entanglement maintained by virtual excitations in an ultrastrongly-coupled-oscillator system

We study the effect of quantum entanglement maintained by virtual excitations in an ultrastrongly-coupled harmonic-oscillator system. Here, the quantum entanglement is caused by the counterrotating interaction terms and hence it is maintained by the virtual excitations. We obtain the analytical expression for the ground state of the system and analyze the relationship between the average excitation numbers and the ground-state entanglement. We also study the entanglement dynamics between the two oscillators in both the closed- and open-system cases. In the latter case, the quantum master equation is microscopically derived in the normal-mode representation of the coupled-oscillator system. This work will open a route to the study of quantum information processing and quantum physics based on virtual excitations.

www.nature.com/scientificreports/ considered in the quantum Rabi model 8 . In addition, the role of the CR terms in the creation of entanglement between two atoms has been investigated in Ref. 49 . We consider an ultrastrongly-coupled two-harmonic-oscillator system. We study the ground state entanglement of the two oscillators and analyze the average excitation numbers in the system. We also study the entanglement dynamics of the system when it is initially in the zero-excitation state and hence all the excitations are created by the CR terms. The influence of the environment dissipations on the system is analyzed based on a microscopically derived quantum master equation in the normal-mode representation. The rest of this paper is organized as follows. Firstly, we present the physical model of two coupled harmonic oscillators and the Hamiltonian, we also analyze the property of the parity chain in this system. Secondly, we obtain the exact analytical eigensystem of the coupled two-oscillator system. Thirdly, the average virtual excitation numbers are calculated analytically and the quantum entanglement of the ground state is analyzed by calculating the logarithmic negativity. Fourthly, we study the dynamics of the average virtual excitation numbers and quantum entanglement between the two oscillators in both the closed-and open-system cases. Finally, we present a brief conclusion.

Results
Model and hamiltonian. We consider an ultrastrong coupling system, in which two harmonic oscillators are ultrastrongly coupled to each other through the so-called "position-position" type interaction (Fig. 1). This system is described by the Hamiltonian where x 1 ( x 2 ) and p 1 ( p 2 ) are, respectively, the coordinate and momentum operators of the first (second) oscillator with the resonance frequency ω 1 ( ω 2 ) and mass µ , the parameter η is the coupling strength between the two oscillators. By expanding the interaction term, Hamiltonian (1) can be expressed as where we introduce the renormalized frequencies and coupling strength as By introducing the following creation and annihilation operators Hamiltonian (2) becomes with C = ( ω a + ω b )/2 being a constant term. Here a † (a) and b † (b) are, respectively, the creation (annihilation) operators of the two oscillators with the corresponding resonance frequencies ω a and ω b . In Eq. (5), the first two terms and the constant term represent the free Hamiltonian of the two oscillators. The parameter g = −ξ/(2µ √ ω a ω b ) denotes the coupling strength between the two oscillators. We note that this interaction includes both the rotating-wave and CR terms. In general, in the case of weak coupling and near resonance, the rotating-wave approximation can be made by discarding the CR terms. In this paper, we consider the ultrastrongcoupling case in which the CR terms cannot be discarded. In the presence of the CR terms, the ground state of the system will include excitations and hence quantum entanglement will exist in the ground state. Note that an ultrastrongly-coupled two-mode system has recently been realized in superconducting circuits 50 . Figure 1. Schematic diagram of the coupled two-harmonic-oscillator system. Two harmonic oscillators with resonance frequencies ω a and ω b are coupled to each other via a "position-position" type interaction with the coupling strength g. The parameters γ a and γ b are the decay rates associated with the heat baths in contacted with the oscillators a and b, respectively. www.nature.com/scientificreports/ In this two-oscillator system, we introduce the parity operator as P = (−1) a † a+b † b , which has the standard properties of a parity operator, such as P 2 = I , P † P = I , and P † = P 7, 51 . The Hamiltonian H in Eq. (5) remains invariant under the transformation P † HP = H , based on the relations P † aP = −a , P † a † P = −a † , P † bP = −b , and P † b † P = −b † . The Hilbert space of the system can be divided into two subspaces with different parities: odd and even. The basis states of the odd-and even-parity subspaces are, respectively, given by and The eigenvalues of the parity operator P corresponding to the odd and even parity states are −1 and 1, respectively.
Eigensystem of the coupled two-oscillator system. To study the quantum entanglement of the eigenstates, we need to diagonalize the Hamiltonian H in Eq. (2). To this end, we introduce the transformation operator 52 where the mixing angle θ is defined by It is obvious that the ground state of the two-oscillator system is |0� A |0� B . To study the virtual excitations in the system, we need to know the eigenstates which are expressed in the representation associated with a † a and b † b . .
(16) H|m� A |n� B = E m,n |m� A |n� B , m, n = 0, 1, 2, . . . ,  www.nature.com/scientificreports/ It can be seen that the superposition components in the ground state are even parity states. This property can be confirmed because the transform U conserves the excitation number and the squeezing operators change the excitation number two by two, without changing the parity.
Ground-state entanglement and quadrature squeezing. We study the ground-state entanglement in this system by calculating the logarithmic negativity. For the two-oscillator system, if the coupling is sufficiently weak, i.e., g ≪ {ω a , ω b } , the interaction Hamiltonian between the two oscillators can be reduced by the RWA as , which conserves the number of excitations. In this case, the ground state of the system is a trivial direct product of two vacuum states |0� a |0� b , which does not contain excitations. In the presence of the CR terms, the |0� a |0� b is not an eigenstate of the system and the ground state will possess excitations. Below, we use numerical method to obtain the ground state of Hamiltonian (5) and calculate the ground state entanglement between the two oscillators. In the presence of the CR terms, the ground state of the two-oscillator system can be expressed as where these superposition coefficients are given by C m,n = a �m| b �n|G� , which should be solved numerically.
The �G|a † a|G� � = 0 and �G|b † b|G� � = 0 reveal that the ground state of the system contains excitations. These excitations in the ground state are called virtual excitations because these excitations cannot be extracted from the system. The effect of the virtual excitations can be seen from the probability amplitudes in the ground state. The distribution of these probability amplitudes can also exhibit the parity of the ground state. As the ground state is an even parity state, and hence these probability amplitudes associated with the odd parity basis states will disappear. In Fig. 2, we show the absolute values of these probability amplitudes |C m,n | . Here we can see that the values of |C m,n | decrease with the increase of m and n and that there is a symmetric relation |C m,n | = |C n,m | . In addition, the values of these odd-parity probability amplitudes C m,n with m + n being an odd number are zero, which is a consequence of the fact that the ground state is an even-parity state.
We also calculate the average excitation numbers a † a and b † b in the ground state |G� as where we have used the formula, In Fig. 3a, we show the average excitation numbers a † a and b † b in the ground state |G� as functions of the scaled coupling strength g/ω r in the degenerate oscillator case ω a = ω b = ω r . These results show that the average excitation numbers of the two modes are identical (two curves overlap each other). This is because the corotating terms conserve the excitations and the CR terms create simultaneously the excitations in the two modes. The The degree of entanglement of the ground state can be quantized by calculating the logarithmic negativity 54,55 . For a bipartite system described by the density matrix ρ , the logarithmic negativity can be defined by where T b denotes the partial transpose of the density matrix ρ of the system with respect to the oscillator b, and the trace norm ρ T b 1 is defined by Using Eqs. (34), (35), and (36), the logarithmic negativity of ground state of the two coupled oscillators can be obtained. In Fig. 3b, we show the logarithmic negativity N as a function of the coupling parameter g/ω r . The curve shows that the degree of entanglement between the two oscillators in the ground state monotonically increases over the entire range of g. This is because the CR terms in Hamiltonian (5) cause the virtual excitations in the ground state of the system and maintain the quantum entanglement between the two oscillators. If the CR terms are discarded, then the ground state of the system becomes a separate state |0� a |0� b .
We also study the quadrature squeezing in the ground state by calculating the fluctuations of the rotated quadrature operators. We introduce the rotated quadrature operators for the two modes as The commutation relation of the above two rotated quadrature operators is According to the uncertainty relation, we have Then the quadrature squeezing appears along the angle θ o if the variances of the rotated quadrature operators satisfy the relation 56 For the ground state |G� given in Eq. (27), the variances of the rotated quadrature operators can be obtained as When we exchange the subscripts a and b in Eq. (41), the expression does not change for a given rotating angle. This means that, in the resonance case ω a = ω b = ω r , the squeezing is the same for the two bosonic modes in  www.nature.com/scientificreports/ the ground state. This point can also be seen from Hamiltonian (5), which is symmetric under the exchange of the subscripts and operators for the two modes in the resonance case. In Fig. 4a, we show the variance �X 2 a (θ a ) as a function of the rotating angle θ a in the resonance case ω a = ω b = ω r . Here we can see that the variance �X 2 a (θ a ) is periodic function of θ a and that the minimum of �X 2 a (θ a ) is obtained at θ a = π/2 and 3π/2 . Note that in the present case (sinh r a cosh r a + sinh r b cosh r b ) < 0 . We also show the variance �X 2 a (π/2) as a function of the coupling strength g/ω r in the resonance case ω a = ω b = ω r , as shown in Fig. 4b. We observe that the squeezing increases with the scaled coupling strength g/ω r . This is because the quadrature squeezing is caused by the CR interaction terms.
Dynamics of quantum entanglement. The phenomenon of quantum entanglement accompanied with virtual excitations can also be seen by analyzing the entanglement dynamics of the system. We consider the case in which the system is initially in the zero-excitation state |0� a |0� b . In the closed-system case, a general state of the system can be written as By substituting Eqs. (5) and (42) into the Schrödinger equation, the equations of motion for these probability amplitudes A m,n (t) are obtained as For the initial state |0� a |0� b , the initial condition of these probability amplitudes reads A m,n (0) = δ m,0 δ n,0 . By numerically solving Eq. (43) under this initial condition, the evolution of these probability amplitudes can be obtained. Using Eqs. (35), (36), and (42), we can calculate numerically the average excitation numbers a † a and b † b and the logarithmic negativity of the state |ψ(t)�.
In Fig. 5a, we show the time evolution of the average excitation numbers a † a and b † b in modes a and b. Here we can see that, similar to the ground state case, the average excitation numbers in the two modes are identical (the two curves overlap each other). In addition, the average excitation numbers experience a periodic oscillation. In Fig. 5b, we show the time dependence of the logarithmic negativity N(t) of the state |ψ(t)� . The curve shows that logarithmic negativity between the two oscillators also experiences a periodic oscillation. Here we choose the initial state of the system as |0� a |0� b , the existence of the CR terms still causes the appearance of virtual excitations, which leads to entanglement between the two oscillators. This result is different from that in the RWA case in which the CR terms are discarded in the two oscillators under the same initial state. When we discard the CR terms and choose the initial state as |0� a |0� b , which is the eigenstate of the corotating interaction term g(a † b + b † a) , the system will stay this state. Then there are no virtual excitations in the system and no quantum entanglement between the two oscillators.
We also study the influence of the environment dissipations on the dynamics of the system. As we consider the ultrastrong-coupling regime of the coupled system, we derive the quantum master equation in the normal-mode representation of these two coupled oscillators. We employ the standard Born-Markov approximation under the condition of weak system-bath couplings and short bath correlation times to derive the quantum master equation. The secular approximation is made by discarding these high-frequency oscillating terms including exp(±iω A t) , exp(±iω B t) , and exp[±i(ω A ± ω B )t] . The quantum master equation in the normal-mode representation of Hamiltonian (5) can be written as   www.nature.com/scientificreports/ Note that the cross terms between the two modes a and b in Eq. (50) are induced by the interaction between the two oscillators. For below calculations, we express the density matrix of the two-oscillator system in the bare-mode representation as with the density matrix elements ρ m,n,j,k (t) = a �m| b �n|ρ s (t)|j� a |k� b . For an initial state |0� a |0� b , the nonzero density matrix element is ρ 0,0,0,0 (0) = 1 . By numerically solving Eq. (47) under the initial condition, the time evolution of the density matrix ρ s (t) can be obtained.
Below we study the dynamics of the average excitation numbers and quantum entanglement in this system. Based on Eq. (47), the expressions of the average excitation numbers a † a(t) and b † b(t) can be expanded as Therefore, the average excitation numbers a † a(t) and b † b(t) can be obtained by solving the equations of motion for these density matrix elements in the number-state representation.
In Fig. 6a, the dynamics of the average excitation numbers a † a(t) and b † b(t) is shown in the open-system case with different time t. We observe that the two excitation numbers a † a(t) and b † b(t) overlap each other and initially experience a large oscillation. With the increase of time t, the oscillation amplitudes of the average excitation numbers decrease gradually. In the long-time limit t ≫ 1/γ a,b , the average excitation numbers will reach steady values due to the dissipations.  www.nature.com/scientificreports/ The entanglement of the density matrix ρ s (t) can be quantified by calculating the logarithmic negativity. In terms of Eqs. (35), (47), and (53), the logarithmic negativity of the state ρ s (t) can be obtained numerically. In Fig. 6b, we show the logarithmic negativity N(t) of the density matrix ρ s (t) versus the time t. The result shows that the logarithmic negativity oscillates very fast due to the free evolution of the system. We also find that the envelope of the logarithmic negativity converges gradually with the evolution time t and eventually reaches a stable value due to the dissipations. The time scale of the oscillation-pattern decay for the logarithmic negativity is very similar to that of the excitations created by the CR interaction terms. In particular, we find that there exists steady-state entanglement due to the presence of the CR interaction terms in this system.
In this work, we consider the ultrastrong-coupling regime and hence the quantum master equation is derived in the normal mode representation. For comparison, we show in Fig. 6c,d the evolution of the average excitation numbers and the logarithmic negativity calculated by solving the phenomenological quantum master equation, which is obtained by adding the dissipators of two free bosonic modes into the Liouville equation, The initial state of the system is the same as that considered in the microscopic quantum master equation. We see from Fig. 6 that, for the average excitation numbers, though these results can approach steady-state values, the envelop and the oscillation amplitude are different for the results obtained with two different quantum master equations. However, for the logarithmic negativity, we find that the difference between the two results exists but is small when g/ω r = 0.2 . We checked the fact that the difference will increase as the increase of the ratio g/ω r . Therefore, the microscopic quantum master equation should be used in the ultrastrongly-coupledoscillator system.

conclusion
In conclusion, we have studied quantum entanglement in an ultrastrongly-coupled two-harmonic-oscillator system. Concretely, we have studied the ground-state entanglement by calculating the logarithmic negativity of the ground state. Here, the quantum entanglement is maintained by the virtual excitations generated by the CR terms and bounded in the ground state. We have also studied the dynamics of quantum entanglement of the system. By microscopically deriving a quantum master equation in the normal-mode representation of the two oscillators, we analyzed the influence of the dissipations on the entanglement dynamics and found that there exists steady-state entanglement in this system.