Instability of Meissner Differential Equation and Its Relation with Photon Excitations and Entanglement in a System of Coupled Quantum Oscillators

: In this work, we investigate the Schrödinger dynamics of photon excitation numbers and entanglement in a system composed by two non-resonant time-dependent coupled oscillators. By considering π periodically pumped parameters (oscillator frequencies and coupling) and using suitable transformations, we show that the quantum dynamics can be determined by two classical Meissner oscillators. We then study analytically the stability of these differential equations and the dynamics of photon excitations and entanglement in the quantum system numerically. Our analysis shows two interesting results, which can be summarized as follows: (i) Classical instability of classical analog of quantum oscillators and photon excitation numbers (expectations (cid:68) N j (cid:69) ) are strongly corre-lated, and (ii) photon excitations and entanglement are connected to each other. These results can be used to shed light on the link between quantum systems and their classical counterparts and provide a nice complement to the existing works studying the dynamics of coupled quantum oscillators.


Introduction
Quantum entanglement is one of the most intriguing concepts of quantum theory [1,2], although it was initially conceptualized and used to refute the basic quantum principles. Indeed, Einstein, Podolsky and Rosen have argued against some results of quantum theory by remarking that the wave functions can be entangled, which entails in their point of view the existence of hidden variables [1,3]. It was later realized that entanglement is an essential component of quantum theory because, without it, it is impossible to interpret and confirm the provisions of the theory [4]. Nowadays, it is widely understood that entanglement plays an important role in quantum information processing protocols, and it is considered a necessary resource to go beyond classical communications and technologies.
During the last years, controlling entanglement in systems of time-dependent coupled harmonic oscillators has been extensively studied [3,[5][6][7][8][9][10][11]. These investigations were also motivated by experiments providing access to the ultra strong coupling (USC) regime [12,13], where the coupling is comparable to the transition frequencies of the coupled quantum systems. In another work, it was shown that it is possible to generate entanglement by phasing control in two [14] and three [15] isotropic harmonic oscillators sinusoidally coupled to each other by J(t) = J 0 cos(ωt) and weakly coupled to a harmonic bath. It was also found that the survival of entanglement for large simulation times is due to the instability of decoupled from the bath normal oscillator modes [3]. The connection between normal mode instability and entanglement growth was also demonstrated in References [16,17], again by using sinusoidal coupling. More recently, it was shown in Reference [18] that the ground state |G of two time-independent resonant oscillators contains, in the USC regime, virtual (cannot be absorbed) excitations VE, i.e., G|N 1 |G = G|N 2 |G = 0, with N 1 and N 2 being the photon number operators corresponding to the coupled oscillators. Note that the term VE is used here to denote virtual photons that exist in the ground state but cannot be detected directly [18], although there are some proposals in order to probe them [19,20]. This is not to be confused with the usage of the term in several other works [21][22][23][24][25] in which qubits detuned from a cavity where they are embedded and interact through the virtual exchange of photons with the cavity. In the weak coupling regime, the ground state of the two coupled oscillatos is empty and does not contain VE, but as the coupling increases and reaches USC it becomes populated [13]. VEs were considered as a consequence of the counter-rotating terms appearing in the Hamiltonian of coupled oscillators, and the maintenance of entanglement between oscillators was attributed to the presence of these excitations.
Motivated by the above studies, in the present work we elucidate the interplay between the instability of classical oscillator equations, which can be associated with the system of time-dependent coupled quantum oscillators under study, the virtual quantum excitations and quantum entanglement. Specifically, we address the question of how classical instabilities affect photon excitations and consequently result in entanglement in the ground state of two time-dependent non-resonant coupled harmonic oscillators. We make these investigations by using a pair of harmonic oscillators connected by a periodically quenched coupling parameter J(t) = J 0 Θ(t) and with time-dependent frequencies where Θ is a π-periodic square wave with alternating value ±1, and is the quench amplitude. Under this choice of time-dependent functions, the resolution of Schrödinger dynamics allows us to find Ermakov equations [26], and the utilization of suitable transformations results in two Meissner differential equations [27] corresponding to classical counterparts of the decoupled Hamiltonian. Then, we study the instabilities of the derived differential equations and obtain stability-instability diagrams, which enable us to investigate the link between two strongly different features, VEs and classical instabilities. Additionally, by computing logarithmic negativity, we establish the relation between entanglement and VE.
The layout of our paper is as follows. In Section 2, we present our coupled oscillator model and show how to exactly decouple the Hamiltonian system by using suitable transformations. The instabilities of Ermakov and Meissner equations emerging from the analysis are discussed in Section 3. In Section 4, we quantify entanglement by using logarithmic negativity and calculate the VE in both oscillators by averaging the number operators over the ground state. In Section 5, we present numerical results for various values of system parameters and discuss the interplay between these phenomena. Finally, Section 6 concludes this work.

Model and Integrability
The model under study is displayed in Figure 1 and consists of two time-dependent non-resonant coupled harmonic oscillators described by the following Hamiltonian: where positionsx j and momentap k and j, k = 1, 2 satisfy the canonical commutation relations x j ,p k = iδ jk ; ω j (t) is the time-dependent frequencies of the oscillators; and J(t) is the time-dependent coupling parameter between them. For the sake of simplicity, we assume unit masses (a transformation that may be applied to the more general case and make this assumption valid [28,29]) andh = 1. Since the time-dependent Hamiltonian (1) involves an interaction term, a straightforward diagonalization is not an easy task. In order to overcome this, we introduce the following time-dependent rotation by the angle α(t) defined below: whereL z =x 1p2 −x 2p1 is the angular momentum operator. The transformed Hamiltonian is given by the following: and after some algebra we obtain the following: where the frequencies Ω j are defined as follows.
For the boundness of this Hamiltonian, it is required that Ω 2 j > 0. From Equation (4), it becomes obvious that the separation of variables is possible forα(t) = 0, and by using Equation (2), this condition results in the following relation: which has been also used in different occasions, see for instance References [28,30]. The original Hamiltonian is, thus, canonically equivalent to the following decoupled Hamiltonian: which can be easily solved.

Time-Dependent Schrödinger Equation
The commutation relation H 1 ,H 2 = 0 implies that the solutions of the time-dependent Schrödinger equation corresponding to the decoupled Hamiltonian (7) have the following form:Ψ where eachφ(x j ; t) satisfies the corresponding single-particle Schrödinger equation.
This equation was studied earlier in References [31][32][33], where the general solution was derived as a superposition of orthonormal expanding modesψ(x j ; t) = ∑ |p n j (t)| 2 = 1. For a single mode and Ω 2 j (0) > 0, the following is the case: with what follows: where H n j is the Hermite polynomial, and the scaled frequencies are defined as j (t) = The functions h j (t) are solutions of the following Ermakov equations (dots stand for time derivatives hereafter):ḧ with initial conditions h j (0) = 1,ḣ j (0) = 0. It is worth it to note that the energy spectrum is time-independent.
On the other hand, the energy expectation value is time-dependent The eigenfunctions of the decoupled Hamiltonian (7) are listed as follows.
In the forthcoming analysis, we concentrate on the following vacuum solution: where the time-dependent parameters involved are read as follows.
These results will be used to compute the occupation numbers and used to discuss entanglement through logarithmic negativity. Note that Solution (16) can also be obtained by using the method presented in Reference [34].

Classical Stability of the Associated Differential Equations
As we have observed above, the ground state (17) depends strongly on the solutions h 1 and h 2 of the Ermakov equations (12). For this reason, it becomes important to study the classical stability of these differential equations. We begin by performing the following transformation on h 1 and h 2 [35]: in order to obtain the following system: describing two decoupled classical parametric oscillators with time-dependent frequencies Ω 1 (t) and Ω 2 (t). The solutions h j (t) of the Ermakov equations can be expressed as follows (see Reference [26]): where x j and y j are independent solutions of Equations (20a) and (20b), respectively, satisfying the initial conditions x 1 (0) = y 1 (0) = 1,ẋ 1 (0) =ẏ 1 (0) = 0 and x 2 (0) = y 2 (0) = 0, while the Wronskians W[x 1 , x 2 ] = x 1ẋ2 − x 2ẋ1 and W[y 1 , y 2 ] = y 1ẏ2 − y 2ẏ1 are both constant. From this, it becomes obvious that the stability of the system (20a) and (20b) determines that of the Ermakov system.
In order to proceed further, we consider a particular shape for the time-dependent frequencies and coupling parameter [36,37], specifically the following modulation: where the quench function Θ(t) is a π-periodic square wave defined by the folloinwg.
ω 0 , J 0 and are the parametric frequency, coupling amplitude and quench amplitude, respectively. With this choice, condition (6) is satisfied, while the frequencies Ω j (5) become the following.
Note that the use of piecewise constant coupling (22) is physically motivated from pulsed laser cooling [38,39] and entanglement generation [8] in cavity optomechanics, where the coupling between the optical cavity mode and the mechanical resonator is modulated by piecewise constant pulses.
Since the modulation in Equation (24) is periodic, Equations (20a) and (20b) are actually Hill equations. This type of equations plays a vital role in the modeling of a variety of dynamical phenomena, which have been studied for more than a hundred years [40]. Investigating the stability of their solutions has a paramount role in the quantum dynamics of open systems [3,14,15]. Since the periodic modulation Θ(t) is a square wave, Equations (20a) and (20b) actually reduce to Meissner equations [37], which is a type of lossless Hill equation that admits an analytical treatment [27]. A Meissner oscillator can be observed, for instance, as an LC-oscillator with parametric frequency ω 2 = (LC) −1 where the capacitance C is pumped by a voltage V(t) such that C(t) = C 0 + C p Θ(t), C p < C 0 [36] or as a charged pendulum in an alternating piece-wise constant homogeneous electric field [37].
Before discussing the stability of these equations, we point out that the allowed values of parameters ω 0 , J 0 and are constrained by the requirement Ω 2 j > 0 expressing the boundedness of the Hamiltonian. Using Equation (24), this condition becomes ω 4 0 > 2 + J 2 0 , which in (J 0 , ) space corresponds to an open disc centered at the origin and with radius ω 2 0 , as displayed in Figure 2a for ω 2 0 = 1. In Figure 2b, we display the range of parameters result in the bounded Hamiltonian in the space (J 0 , ω 2 0 ) and for fixed = 1. These plots show that the regions corresponding to unbounded Hamiltonians occupy larger portions in parameter space, something which limits our choices. Note that the edges of boundness present great importance, for example, in generating entanglement and for inverse engineering of time-dependent coupled harmonic oscillators [41].
We now turn our attention to the stability of Meissner equations under consideration. The stability of a π-periodic second order differential equation is determined by the state transition matrix corresponding to the first period of evolution [0, π] [40]. The reason is that the solution of such an equation at t = mπ + τ, with m positive integer and τ ∈ [0, π), is . In order to find this transition matrix, we rewrite System (20a) and (20b) by exploiting the piecewise constant nature of function Θ(t). Since Equations (25a) and (25b) are decoupled, we treat them independently. Their Wronskians matrices [36] for t ∈ [0, π] are given by the following: where Ω 0 j := Ω j (0) = ω 2 0 ± 2 + J 2 0 . The state transition matrix in the interval [t 1 , t 2 ] is the following [36]: with W being the adjoint matrix. By making use of Equation (27), one can easily show that the following is the case: with τ ∈ [0, π]. After straightforward algebra and using det W(t) = 1, we obtain the state transition matrix for each of Equations (20a) and (20b).
From Equation (23) we have τ = π/2, in which case T X = T Y := T . Consequently, Equations (29a) and (29b) result in the same stability condition |Tr T | < 2 [36], correspond-ing to |λ i | < 1 for the eigenvalues of the transition matrix, which can be evaluated to obtain the following: where Λ is the dimensionless parameter.
As a matter of convenience, we define the stability function as follows.
Now, it is clear that the classical system is stable if S > 0, otherwise it is unstable.
In Figure 3, we show the numerically obtained stability diagram in the ( , ω 2 0 ) space, for two different values of the coupling, J 0 = 0.2 ( Figure 3a) and J 0 = 0 (Figure 3b). From this, it becomes clear that the stability diagram is very sensitive to the coupling, since a small change in the latter can result in important configuration variation. We notice also that the instability region increases as long as the coupling parameter increases, and this observation motivates us further to investigate the relation between classical instabilities, VE and entanglement.

Entanglement and Dynamics Effects
According to the Peres-Horodecki criterion [42,43], the necessary and sufficient conditions for the separability of two Gaussian mode states include the positivity of the partially transposed state. Since the vacuum state (17) is pure and symmetric, then separability can be realized by setting the coupling term (18c) equal to zero, namely having the following condition: which is necessary and sufficient for separability. Obviously, the above condition can be met in two subordinate cases: (i) for J 0 = 0, where the dynamics cannot generate entanglement and the oscillator states remain separable; and (ii) for the Wronskian W[h 1 (t), h 2 (t)] = 0 and 1 (t) = 2 (t), where the geometrical interpretation of W is the difference of the rectangular phase space areas S 1 = h 1ḣ2 and S 2 = h 2ḣ1 [41]. The last case indicates that dynamics can extinct entanglement and proper engineering of the initial and final normal frequencies Ω j (0), Ω j (t f ) and j = 1, 2 is required in order to avoid this . Regarding our case, the vacuum state is completely described by the marginal purities µ j and j = 1, 2.
Since our state is pure and Gaussian, all quantum correlations can be derived from its second moment, that is, the covariance matrix V (t). The covariance matrix can be transformed via a local symplectic transformation S = S 1 ⊕ S 2 to a particular form called standard form V s f (t) [42].
By performing the partial transposition prescription, det(A) −→ det(A) and det(C) −→ − det(C), the minimal symplectic eigenvalue of the partially transposed covariance matrix V is as follows: After some manipulation, we find the following: where the symplectic invariant is ∆(Ṽ ) = 2(det(A) − det(C)). Entanglement is quantified through logarithmic negativity: which is monotonically increasing with |A 12 | 2 . From the expression (18c) of the latter quantity, it appears that the main contributions of the time-dependence in the Hamiltonian is the emergence of an imaginary part, that is,ḣ 2 /h 2 −ḣ 1 /h 1 , as well as the initial normal mode scaling, i.e., Ω j (0) −→ j (t) = Ω j (0)/h 2 j (t).

Virtual Excitations
By using phase space prescription [44], we will analyze VE by computing the expectation value of photon numbers a + j a j in the vacuum state. The case of two resonant time-independent coupled oscillators was analyzed in [18]. where the creation a + j,0 and annihilation a j,0 operators are simply mapped as follows.
(a + j,0 ) + = a j,0 = However, for a time-dependent Hamiltonian, the situation is not that obvious because the realization can be performed as follows [45,46]: where a j , a + k = δ jk , η j (t) = ω j (0) ν 2 j and the functions ν j satisfy the Ermakov equations: with initial conditions ν j (0) = 1,ν j (0) = 0. After some straightforward algebra, one can compute the expectation of photon number operators N j = a + j a j and find the following: where and k = 1, 2 with k = j. We emphasize that the time-dependent Hamiltonian (1) generates important effects such as the scalings η j (t) = ω j (0) ν 2 j and j (t) = Ω j (0) h 2 j (t) , as well as shiftings in positions (43a) and momenta (43b), which are due to the dilatation functions ν j (t) and h j (t). In addition, the appearance of the term x jpj t is purely a consequence of time dependence: observe the factorν j in Equation (43c), which of course disappears in the time-independent case. In the latter time-independent case, we have ν j = h j = 1 andν j =ḣ j = 0, and the VE becomes the following.
We observe that the study in Reference [18] was confined to the resonant case because it was not easy to disentangle the rotation operator (2) in the frame of creation a + j and annihilation a j representation. Now, it becomes clear from our analysis how the phase space picture can be used to overcome such situations. In addition, it is interesting to note that in the time-dependent case for ω 1 = ω 2 ( = 0) and J 0 = 0, the vacuum state (17) may contain excitations, i.e., N j = 0, even if the oscillators are decoupled, while the solution for the wavefunctions is given by formulas analogous to Equations (10) and (11) since ω j is time dependent (see Equation (22) with = 0). Such phenomena do not exist in the frame of resonant and time-independent oscillators.

Results and Discussions
For a numerical study of classical instability effects on generation of photon excitations and entanglement between oscillators, it is convenient to define the dimensionless parameters of the following: where Ω 0 is an arbitrary frequency. Before presenting and discussing the main numerical results, we also provide analytical expressions for the dilatation functions h j and ν j in the first period [0, π]. By using Equation (21), one can solve Equation (12) to obtain the following. (47b) Note that the solutions for ν j can be immediately obtained from Equations (47a) and (47b) with the substitution 2 + J 2 0 → . We also observe that in odd half periods the dynamics will be freezed, while the solutions will evolve dynamically in even half periods due to the periodic quench and the continuity of solutions h j and ν j .

Instability Versus VE
In this paragraph, we investigate the relation between instability of the oscillator dynamics and the VE generated in the ground state. First, we use Equation (32) to plot in Figure 4 the stability map in the ( , J 0 ) space for fixed ω 2 0 = 1.01. Observe that the system remains stable for a small area centered around the origin = J 0 = 0, where the stability function given in Equation (32) is positive, S > 0. Next, in Figure 5, we plot the dynamics of VE for various values of parameter , J 0 , corresponding to stable and unstable oscillators. We start in Figure 5a by setting the quench parameter and the coupling J 0 both equal to 0, in which case the oscillators are stable, and the VE can not be created. Next, in Figure 5b we increase to 0.01 while maintaining the oscillators as decoupled, J 0 = 0. The dynamics still remains stable while the ground state remains empty of VE. The situation is similar in Figure 5c where we set = 0 while increasing the coupling to the small value J 0 = 0.005 and where both paramaters are set to = J 0 = 0.005 in the similar case in Figure 5d .  Until now, the observed behavior is obvious and normal as long as the oscillators are weakly coupled and the virtual excitations are considered as an intrinsic feature of ultrastrong coupling (0 J 0 < ω 2 0 ). However, what is interesting is that, when we increased to 0.1 and maintained the oscillators weakly coupled with J 0 = 0.005, the dynamics becomes unstable, and the ground state becomes populated, as shown in Figure 5e. This feature does not show up when the oscillators are resonant (ω 1 = ω 2 ) and weakly coupled (0 < J 0 ω 2 0 ). This behavior is due to the time dependence of the Hamiltonian that maintains the out-of-resonance term (ω 2 1 − ω 2 2 ) 2 = 4 2 , which has a contribution to the mode frequencies Ω 2 j similar to that of the coupling J 0 (see Equation (24)). For completeness, in Figure 5f we showed the VE for the resonant case = 0 but with stronger coupling J 0 = 0.1, and we observed that the ground state becomes also populated. We additionally displayed the cases = 0.1, J 0 = 0.1 and = 0, J 0 = 0.2 in Figure 5g,h, respectively, which display VE.
From the plots of Figures 4 and 5, we conclude that, when the oscillators are stable, the ground state is empty of VE. Additionally, as we move away from the stability region, the ground state becomes more populated with a sea of VE. This leads us to consider using the classical instability as a signature of VE in a quantum system. Moreover, it also provides an easier method to detect the generation of excitations by engineering only the classical counterpart of the quantum system. This connection between quantum systems and their classical counterparts may also shed light on the nature of the ground state VE [3].

Entanglement and VE
In this subsection we study the connection between quantum entanglement and VE. First, we confine ourselves to the case of resonant oscillators ( = 0) in the USC regime, where coupling J 0 has the same order of magnitude with transition frequency ω 2 0 such that the following is the case.
Thus, classical dynamics is always unstable, S = 0. Note that the modulator (e.g., pumping source) has a frequency ω m = 2, while for = 0 the parametric oscillators have the same frequency, and this explains why the system is unstable for all J 0 > 0 at the resonance point ω m = 2ω 0 . Near the resonance point, for example, when ω 2 0 = 1.01, the instability appears for weak coupling and is maintained in the USC regime.
In Figure 6 we plot the dynamics of entanglement and VE for three values of ultrastrong coupling, J 0 = 0.1, 0.2, 0.3, and observe that entanglement and ground state excitations always coexist. Consequentely, one can use the classical instability as a signature of virtual excitations and entanglement in the USC regime. Furthermore, entanglement and VE exhibit similar dynamics. Indeed, before the quench time t q = π/2, the dynamics freezed, which is obvious since the Ermakov solutions (47a) and (47b) are constant and equal to one. Immediately after the quench, both dynamics revived exponentially. This similarity indicates that virtual excitations are necessary to maintain entanglement. Indeed, by increasing the coupling, the ground state populations of both parts increase, and the exchange of these VE, similarly to the Casimir effect, plays a vital role in mediating the coupling between oscillators [47] and, consequently, the maintenance of entanglement between them.
We also considered the non-resonant case where the quench parameter is nonzero, = 0, while we set ω 2 0 = 1.01. In Figure 7, we displayed the geometric average of excitations M 12 (t) = ( N 1 (t) N 2 (t) ) 1 2 as well as entanglement for two nonresonant cases, one with = 0.9, J 0 = 0.1 (red solid line) and another with = 0.99, J 0 = 0.01 (blue solid line). In both cases, we observed that when t ∈ [0, π/2], a tiny and constant number of excitations and entanglement are generated, while for t > π/2 the excitations dramatically increase compared to the entanglement. Before the quench, the VE are too tiny due to the weakness of the coupling and the Ermakov solutions, which are freezed to one, but the excitations drastically revived once the dynamics became quenched since the Ermakov solutions became unstable. The entanglement maximum occurs when M 12 is minimal, around the time t = 7π/8. Since entanglement is a genuine quantum feature, it is expected to be maximized when the occupation numbers of the modes are relatively small. Larger mean occupation numbers, such as those shown in Figure 7b, correspond to more classical behavior; thus, it is more difficult for the modes to be entangled. Note also that in cavity optomechanics, the standard procedure to calculate optomechanical entanglement is to linearize the equations around the classical mean values and then to find the correlations of zero mean quantum fluctuations around them (see, for example, Reference [17]). For the case with strong coupling J 0 = 0.1, the maximum value of entanglement is larger than the value obtained above in the resonant case = 0 with the same coupling (compare the red line in Figure 7a with the blue solid line in Figure 6a). On the other hand, for the case with weak coupling J 0 = 0.01, the presence of the time dependence results in the generation of some small amount of entanglement.

Quantifying Virtual Excitations and Entanglement via Instability
In this paragraph, we show how classical instabilities are deeply related to VE and entanglement of the ground state. To do so, we introduce the following instability quantifier: which quantifies how far the classical dynamics evolves from the stability edge (S = 0). In Figure 8, we plot entanglement E N , the geometric average of excitations M 12 = ( N 1 N 2 ) 1 2 and Ins versus the coupling strength J 0 , for = 0 ( Figure 8a) and = 0.1 (Figure 8b). We observe that the three quantities exhibit a similar dependence on the coupling J 0 . Specifically, they are monotonically increasing while attaining larger values in the USC regime. It is also worth it to note that Ins(J 0 = 0, = 0.1) ∼ 0.005 (oscillators are unstable), E N = 0 and M 12 = 0, while Ins(J 0 = 0, = 0) ∼ −1.2.10 −4 (oscillators are stable), E N = 0 and M 12 = 0. Finally, we point out that those similarities show that classical instabilities can be used as a signature to quantify and detect entanglement along with VE.

Hierarchy of Virtual Excitations
Here, we investigate the effects of coupling J 0 and quench parameter on the hierarchy of VE for both modes.
In Figure 9, we plotted the dynamics of excitations for the following two cases: 1.
In the figures of the first column, we set the quench parameter to 0.1 and use three values for the coupling, 0.1, 0.3 and 0.9 from top to bottom. For J 0 = 0.1, excitations in mode 1 exceed those of mode 2 during the entire interval. Then, by increasing the coupling to 0.3, mode 1 remains initially more populated than mode 2, but immediately after the quench t q = π 2 the hierarchy is inverted, and mode 2 becomes more populated up to some point where the hierarchy is inverted again. By further increasing the coupling to 0.9, the inversion of hierarchy becomes more pronounced until about t ∼ 15π 16 , where the excitations exhibit an inverse monotony. This phenomenon of hierarchy inversion can be observed as a redistribution of excitations, which becomes more important with the increase in USC coupling.

2.
In the figures of the second column, we fixed the ultra-strong coupling to 0.15 and varied the quench factor between the values 0.01, 0.05 and 0.08 from top to bottom. For small quench = 0.01, the excitations in both modes are close to each other, while by increasing the quench the excitations became more discernible but without inverted hierarchy.
To conclude, the hierarchy of VE becomes distinguishable by increasing the quench , i.e., when moving away from resonance. The increase in coupling inverts this hierarchy due to the exchange and redistribution of excitations.
(e) (f) Figure 9. Effects of coupling J 0 (first column) and quench parameter (second column) on the hierarchy of VE in both modes. The red dashed line corresponds to N 1 (t) and the blue solid line to N 2 (t) , while ω 2 0 = 1.01.

Conclusions
We have studied a system composed of two non-resonant coupled harmonic oscillators connected with periodically pumped coupling and frequencies with all of these time-dependent parameters modulated by a square wave. We have solved exactly the Schrödinger dynamics, which eventually results in the Ermakov differential equations. After suitable transformations, we have shown that these classical equations actually determine the evolution of the decoupled quantum Hamiltonian. Under the square wave time modulation, the stability study of these differential equations results in Meissner differential equations, and we have derived the condition of instability and the instability maps for them. We have shown with numerical simulations that, in the weak coupling regime and when the classical dynamics is stable, the two coupled quantum modes are empty of VE, but when increasing the coupling until the USC regime, the dynamics become unstable, and the ground state of the system becomes entangled and populated with VE. Additionally, in the USC regime, entanglement and VE exhibit a similar dynamics, which indicates that excitations are necessary to maintain entanglement between the modes. By defining an instability based quantifier, it is possible to identify entanglement and VE in the USC regime. Finally, we showed numerically that the hierarchy of VE becomes distinguishable as we move away from resonance, while it is inverted as the coupling strength is increased because of the exchange and redistribution of excitations. The present work, in addition to its theoretical merit, may be exploited to evaluate the entangling capacity of a quantum system composed by two coupled oscillators by simply examining the stability of a Meissner classical differential equation.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: