Dynamics of entanglement between two harmonic modes in stable and unstable regimes

The exact dynamics of the entanglement between two harmonic modes generated by an angular momentum coupling is examined. Such system arises when considering a particle in a rotating anisotropic harmonic trap or a charged particle in a fixed harmonic potential in a magnetic field, and exhibits a rich dynamical structure, with stable, unstable and critical regimes according to the values of the rotational frequency or field and trap parameters. Consequently, it is shown that the entanglement generated from an initially separable gaussian state can exhibit quite distinct evolutions, ranging from quasiperiodic behavior in stable sectors to different types of unbounded increase in critical and unstable regions. The latter lead respectively to a logarithmic and linear growth of the entanglement entropy with time. It is also shown that entanglement can be controlled by tuning the frequency, such that it can be increased, kept constant or returned to a vanishing value just with stepwise frequency variations. Exact asymptotic expressions for the entanglement entropy in the different dynamical regimes are provided.


I. INTRODUCTION
The investigation of entanglement dynamics and growth in different physical systems is of great current interest [1][2][3]. Quantum entanglement is well known to be an essential resource for quantum teleportation [4] and pure state based quantum computation [5], where its increase with system size is necessary to achieve an exponential speedup over classical computation [6,7]. And a large entanglement growth with time after starting from a separable state indicates that the system dynamics cannot be simulated efficiently by classical means [8], turning it suitable for quantum simulations.
The aim of this work is to examine the dynamics of the entanglement between two harmonic modes generated by an angular momentum coupling, and its ability to reproduce typical regimes of entanglement growth in more complex many body systems, when starting from an initial separable gaussian state. The latter can be chosen, for instance, as the ground state of the non-interacting part of the Hamiltonian, thus reproducing the typical quantum quench scenario [1,2,8]. The present system can be physically realized by means of a charged particle in a uniform magnetic field within a harmonic potential or by a particle confined in a rotating harmonic trap [9][10][11][12], where the field or rotational frequency provides an easily controllable coupling strength. Accordingly, it has been widely used in quite different physical contexts, such as rotating nuclei [11,12], quantum dots in a magnetic field [13] and fast rotating Bose-Einstein condensates within the lowest Landau level approximation [14][15][16][17][18][19]. In spite of its simplicity, the model is able to exhibit a rich dynamical structure [20], with both stable and distinct types of unstable regimes, characterized by bounded as well as unbounded dynamics, when considering all possible values of the field or frequency in a general anisotropic potential. Nonetheless, being a quadratic Hamiltonian in the pertinent coordinates and momenta, the dynamics can be determined analytically in all regimes, and the entanglement between modes can be evaluated exactly through the gaussian state formalism [21][22][23][24][25]. For the same reason, the Hamiltonian is also suitable for simulation with optical techniques [26].
The main result we will show here is that due its nontrivial dynamical properties, the entanglement dynamics in the previous model can exhibit distinct regimes, including a quasiperiodic evolution in dynamically stable sectors, different types of logarithmic growth at the border between stable and unstable sectors (critical regime) and a linear increase in dynamically unstable sectors. The model is then able to mimic the three typical regimes for the entanglement growth with time after a quantum quench, arising in spin 1/2 chains with Ising type couplings, according to the results of refs. [1,8], which show a linear growth for short range couplings, a logarithmic growth for long range interactions and an oscillatory behavior for nearly infinite range interactions, when considering a half-chain bipartition. We also mention that the static ground state entanglement of the present model also exhibits critical behavior at the border of instability [27]. Mode entanglement dynamics in related harmonic models within stable regimes were previously studied in [28][29][30], while critical behavior and entanglement in ultrastrong-coupled oscillators (through a different interaction) were considered in [31]. Other relevant aspects of entanglement dynamics and generation in spin systems were discussed in [32][33][34][35].
In sec. II we discuss the exact dynamics of the system and describe the different regimes arising for strong coupling in anisotropic potentials. The entanglement evolution in gaussian states is then examined in detail in sec. III, including its exact evaluation through the covariance matrix formalism and the exact asymptotic behavior in the distinct dynamical regimes. Explicit results, including the possibility of entanglement control through a stepwise varying frequency, are also shown. Conclusions are finally provided in IV.

A. Hamiltonian
We consider two harmonic systems with coordinates and momenta Q µ , P µ , µ = x, y, coupled through their angular momentum L z = Q x P y − Q y P x . The Hamiltonian is Eq. (1) describes, for instance, the motion in the x, y plane of a particle of charge e and mass m within a harmonic trap of spring constantsK µ in a uniform field H along the z axis [11,12], if Ω = e|H| 2mc stands for half the cyclotron frequency and K µ =K µ + mΩ 2 .
It also determines the intrinsic motion of a particle in a harmonic trap with constants K µ which rotates around the z axis with frequency Ω. In this case [11,12], the actual Hamiltonian is H(t) = R(t)H 0 R † (t), with R(t) = e −iΩLzt/ the rotation operator, but averages of rotating observables O(t) = R(t)OR † (t) evolve like those of O under the time-independent "cranked" Hamiltonian (1).

B. Exact evolution
The Heisenberg equations of motion ido/dt = −[h, o] for the operators q µ , p µ (with t = Ω 0 T and T the actual time) become dqx dt = p x + ωq y , and can be written in matrix form as The system dynamics is then fully determined by the matrix H. We may write the general solution of (6) as where O ≡ O(0). In spite of their simplicity, Eqs.(5) can lead to quite distinct dynamical regimes according to the values of ω and k µ , as the eigenvalues of H, which is in general a non-hermitian matrix, can become imaginary or complex away from stable regions [20]. Moreover, H can also become non-diagonalizable at the boundaries between distinct regimes, exhibiting non-trivial Jordan canonical forms [20]. Nonetheless, as the eigenvalues of H are determined by 2 × 2 blocks, and given by λ ± and −λ ± , with where ε ± = kx±ky 2 and ∆ = ε 2 − + 4ω 2 ε + . We can then write the solution (8) explicitly as where uxx yy = (∆±ε−) cos λ+t+(∆∓ε−) cos λ−t

2∆
, The matrix U(t) is real for any real values of ω, k µ and t, including unstable regimes where ∆ and/or λ ± can be imaginary or complex [20]. It represents always a linear canonical transformation of the q µ , p µ , satisfying (I denotes the 2 × 2 identity matrix) which ensures the preservation of commutation relations ( It corresponds to a proper Bogoliubov transformation of the associated boson operators. For ω = 0, we recover from Eqs. (11)-(12) the decoupled harmonic evolution q µ (t) = q µ cos On the other hand, in the isotropic case k x = k y = k (where ∆ = 2ω √ k and |λ ± | = | √ k ± ω|), [l z , h] = 0 and the evolution provided by Eqs. (11)-(12) is just the rotation of identical single mode evolutions: In particular, the Landau case (free particle in a magnetic field) corresponds to k x = k y = ω 2 , where λ + = 2ω and λ − = 0.

C. Dynamical regimes
The distinct dynamical regimes exhibited by this system for ω = 0 are summarized in Fig. 1. Let us first consider the standard stable case k x > 0, k y > 0 (first quadrant). The eigenvalues λ ± are here both real and non-zero in sectors A and B, defined by when k x > 0, k y > 0. A is the full stable sector where h is positive definite, whereas B is that where the system, though unstable, remains dynamically stable [20] (see also Appendix). If ω 2 lies between these values (sector D), λ − becomes imaginary (with λ + remaining real), leading to a frequency window where the system becomes dynamically unstable (unbounded motion), with sin(λ − t)/λ − = sinh(|λ − |t)/|λ − | in Eqs. (12). At the border between D and A or B (ω 2 = k x or ω 2 = k y ), λ − = 0 (with λ + > 0) and H becomes nondiagonalizable if k y = k x , although H 2 remains diagonalizable. The system becomes here equivalent to a stable oscillator plus a free particle [20] (see Appendix), and we should just replace sin(λ − t)/λ − by its limit t in Eqs. (12), which leads again to an unbounded motion.
Considering now the possibility of unstable potentials (k x < 0 and/or k y < 0, remaining quadrants), the dynamically stable sector B extends into this region provided k x > 0 > k y > −3k x (or viceversa) and where the upper bound applies only when ε + < 0 (i.e., −3k x < k y < −k x or viceversa). Eq. (17) defines a frequency window where the unstable system becomes Dynamical phase diagram of the system described by Hamiltonian (1). The evolution of the operators qµ, pµ is quasiperiodic in the dynamically stable sectors A, B, where the eigenfrequencies λ± are both real, but unbounded in the remaining sectors, with λ− imaginary in D, both λ± imaginary in C, and λ± complex conjugates in E. At the borders between these regions (except from the Landau point F ) the matrix H is non-diagonalizable and the evolution is also unbounded, with λ− = 0 at the borders between D and A or B, λ+ = 0 at the border between D and C, λ+ = λ− at the curve separating E from B and C and λ± = 0 at the critical points L. Dashed lines indicate the path described as ω is increased at fixed kµ, showing that the system dynamics will become unbounded (bounded) in a certain frequency window when starting at 1 (2). The black line indicates the isotropic case kx = ky, where entanglement will be periodic (see sec. III).
Finally, if both ∆ and λ = ε + + ω 2 vanish, which occurs when ε + = −ω 2 = −|ε − |/2, i.e., (or ω 2 = k y = −k x /3), the system exhibits a remarkable critical point (points L in Fig. 1), where λ ± = 0 and sectors B, C, D and E meet. Here both H and H 2 are non-diagonalizable, with H represented by a single 4 × 4 Jordan Block (inseparable pair [20]). By using this form or taking the λ → 0 limit in Eqs. (19), we obtain in this case a purely polynomial (and hence also unbounded) evolution, involving terms up to the third power of t: The elements of U(t) become Nonetheless, we remark that Eq. (13) remains satisfied (in both cases (19) and (21)).

A. Exact evaluation
Let us now consider the evolution of the entanglement between the x and y modes, starting from an initially separable pure gaussian state. Since the evolution is equivalent to the linear canonical transformation (8), the state will remain gaussian ∀ t, which entails that entanglement will be completely determined by the pertinent covariance matrix [22,23].
We may then assume that at t = 0, q µ = p µ = 0 for µ = x, y ( O = 0), such that these mean values will vanish ∀ t ( O t = O(t) = 0, as implied by Eq. (11)). We may then define the covariance matrix as x q x q y qxpx+pxqx 2 q x p y q x q y q 2 y q y p x qypy+py qy 2 qxpx+pxqx 2 q y p x p 2 x p x p y q x p y qypy+pyqy 2 which, according to Eqs. (8) and (13), will evolve as The entanglement between the two modes will now be determined by the symplectic eigenvaluef (t) = f (t)+1/2 of the single mode covariance matrix is a non-negative quantity representing the average boson occupation a † µ (t)a µ (t) of the mode (a µ (t) is the local boson operator satisfying a 2 µ (t) = 0), which is the same for both modes (f x (t) = f y (t)) when the global state is gaussian and pure. It is given by Eq. (24) is just the deviation from minimum uncertainty of the mode, and can be directly determined from the elements of (23). The von Neumann entanglement entropy between the two modes becomes where ρ µ (t) denotes the reduced state of the mode. Eq. (25) is an increasing concave function of f (t). For future use, we note that for large and small f (t), Other entanglement entropies, like the Renyi entropies S α (t) = ln Tr ρ α µ (t) 1−α , α > 0, and the linear entropy S 2 (t) = 1 − Trρ 2 µ (t) (of experimental interest as Trρ 2 and in general Trρ n can be measured without performing a full state tomography [3,36]), are obviously also determined by f (t), since Tr ρ α µ = [(1 + f µ ) α − f α µ ] −1 (α > 0). The initial covariance matrix C(0) will be here assumed of the form where α µ = 2 p 2 µ , such that α µ = k µ if the system is initially in the separable ground state of h 0 , as in the typical quantum quench scenario [1]. For fixed isotropic initial conditions we will just take α x = α y = 1.
For these initial conditions, we first notice that for small t, Eqs. (12) and (24) yield which indicates a quadratic initial increase of f (t) with time for any anisotropic initial covariance. Eq. (29) is 1 Sectors Average occupation Entanglement entropy TABLE I. The asymptotic evolution of the average occupation (24) and the entanglement entropy (25) in the different dynamical sectors indicated in Fig. 1. Entanglement is bounded in the stable sectors A, B, but increases linearly (with t) in the unstable sectors C, D, E, and logarithmically at the border between stable and unstable sectors, provided kx = ky.
In the isotropic case kx = ky = k it remains periodic for any value of k and anisotropic initial conditions. independent of the oscillator parameters k µ and proportional to ω 2 . However, for isotropic initial conditions α x = α y , quadratic terms vanish and we obtain instead a quartic initial increase, driven by the oscillator anisotropy ε − : Eq. (27) implies a similar initial behavior (except for a factor ln t) of the entanglement entropy. Next, in the isotropic case k x = k y = k, the exact expression for f (t) becomes quite simple, since the rotation is decoupled from the internal motion of the modes (Eq. (14)), and entanglement arises solely from rotation and initial anisotropy. We obtain Entanglement will then simply oscillate with frequency 4ω if α x = α y , being independent of the trap parameter k, since the latter affects just a local transformation decoupled from the rotation. Eq. (31) holds in fact even if k becomes negative (unstable potential) or vanishes.
In the general case, the previous decoupling no longer holds and the explicit expression for f (t) becomes quite long. The main point we want to show is that the different dynamical regimes lead to distinct behaviors of f (t), and hence of the generated entanglement entropy S(t), which are summarized in Table I. We now describe them in detail.

B. Evolution in stable sectors
In the dynamically stable sectors A and B of Fig. 1, both λ ± are real and non-zero, implying that the evolu- tion of f (t) and S(t) will be quasiperiodic, as seen in Fig.  2 (curves A 1 , A 2 and B). The initial state was chosen as the ground state of h 0 (α µ = k µ in (28)). Starting from point 1 in sector A (Fig. 1), the generated entanglement S(t) remains small when ω is well below the first critical value ω y = k y (curve A 1 ). As ω increases, S(t) will exhibit increasingly higher maxima, showing a typical resonant behavior for ω close to ω y (border with sector D), where λ − vanishes. Near this border, S(t) will essentially exhibit large amplitude low frequency oscillations determined by λ − , with maxima at t ≈ t m = mπ 2λ− (m odd), plus low amplitude high frequency oscillations determined by λ + , as seen in curve A 2 .
As ω increases, the system enters dynamically unstable sectors for ω y ≤ ω ≤ ω x = √ k x , and the evolution becomes unbounded (curves A-D, D 1 and D 2 , described in next subsection). For ω > ω x , the system reenters the dynamically stable regime and exhibits again the previous behaviors, with an oscillatory resonant type evolution for ω above but close to ω x (curve B in Fig. 2).
Close to instability but still within the stable regime, the maximum entanglement reached is of order ln |ω − ω µ |: For ω close to ω µ (µ = x, y) on the stable side, and for the initial conditions (28), f (t) will be maximum at t ≈ t m , with where λ + ≈ 2(ε + + ω 2 µ ) and Expression (32) (and hence S(t m )) will tend to decrease for decreasing anisotropy, i.e., increasing ratio k y /k x ≤ 1, as seen in Fig. 3 for m = 1, vanishing in the isotropic limit k y /k x → 1 (where f (t m ) = O(|k x − k y |) 1/2 ). On the other hand, the behavior for k y /k x → 0 will depend on the initial condition: If it is the ground state of H 0 (α µ = ω µ , curves a), f (t m ) will vanish at the , as obtained from Eq. (32). If the initial state is fixed, however, f (t m ) will approach a finite value for k y /k x → 0, and exhibit a monotonous decrease on average with increasing ratio k y /k x in both borders (curves b in Fig. 3), as also implied by (32). We also mention that the high frequency oscillations in f (t m ) and S(t m ) observed in Fig. 3 stem from the λ + /λ − ratio in the arguments of the trigonometric functions in Eq. (32). For ω close to ω µ , this ratio is minimum around k y /k x ≈ 1/5, which leads to the observed decrease in the oscillation frequency of S(t m ) in the vicinity of this ratio (top panel). frequency ω, for the oscillator parameters and initial state of Fig. 2. Entanglement is bounded for ω < ωy (sector A) and ω > ωx (sector B), but is proportional to t in the instability window ωy < ω < ωx (sector D).

C. Evolution in unstable sectors
Let us now examine in detail the evolution of S(t) in the dynamically unstable regimes. At the critical frequencies ω = ω µ , µ = y, x (borders A-D and B-D), λ − vanishes and Eqs. (12) and (24) lead, for large t and the initial conditions (28), to the critical evolution where λ + = 2(ε + + ω 2 µ ) > 0. This entails a linear increase, on average, of f (t) in this limit, and hence, a logarithmic growth of S(t), according to Eq. (26): where S 0 (t) = 1 + ln[f (t)/t] is a bounded function oscillating with frequency λ + . This behavior (curve A-D in Fig. 2) is the ω → ω µ limit of the previous resonant regime.
On the other hand, in the unstable sector D (ω y < ω < ω x ), λ − becomes imaginary. This leads to an exponential term in f (t) ( sin λ−t λ− → sinh |λ−|t |λ−| ), which will dominate the large t evolution: In this sector Eqs. (12), (24) and (27) imply, for large t, and hence, a linear growth (on average) of the entanglement entropy with time (curves D 1 , D 2 in Fig. 2). Therefore, in the unstable window ω y ≤ ω ≤ ω x , there is an unbounded growth with time of the entanglement entropy, which will originate a pronounced maximum in the generated entanglement at a given fixed time and anisotropy as a function of ω, as appreciated in Fig. 4. We now examine the behavior at the other sectors of Fig. 1. In the unstable sectors C and E, where one or both of the constants k µ are negative, λ ± are imaginary or complex (Fig. I). This implies an exponential increase of f (t), as indicated in table I, entailing again a linear asymptotic growth of the entanglement entropy with time: S(t) ≈ (|λ + | + |λ − |)t in C and S(t) ≈ 2|Im(λ ± )|t in E, neglecting constant or bounded terms.
On the other hand, at the border between sectors B and E, which corresponds to the critical curve ∆ = 0 between both points L in Fig. 1, we obtain, for large t and k x = k y (with the initial conditions (28)), the asymptotic behavior where λ = ε + + ω 2 > 0. This leads to √ αxαy ]. Hence, the unbounded growth of f (t) and S(t) is here more rapid than that at the previous borders A-D and B-D (ω = ω y or ω x ) (quadratic instead of linear increase of f (t)). At the border E-C the asymptotic behavior of f (t) is still exponential (i.e., linear growth of S(t)).
Finally, a further remarkable critical behavior arises at the special critical points L, obtained for condition (20), where all sectors B, C, D and E meet. We obtain here a purely polynomial evolution of (f (t) + 1/2) 2 , as implied by Eqs. (21). For large t, this leads to a quartic increase of f (t):  6. Evolution of the entanglement entropy for a stepwise varying frequency ω, starting from the separable ground state of H0 (with ky = 0.5kx > 0). In curve a we have set ω/ωx = 0.5, 0.7, 0 and 0.21 for successive time intervals of length ωx∆t = 30, such that the system is close to the first instability at the second interval (0.7ωx ≈ 0.99ωy, with ωµ = kµ), while in curve b the only change is ω = 0.75ωx ≈ 1.06ωy in the second interval, such that the system enters there the unstable regime leading to a linear entanglement growth. This plot shows that entanglement can be increased, kept constant and returned to a vanishing value just by tuning the frequency ω.
implying the following logarithmic increase of S(t): where S 2 ≈ 1+ln[ αxαy+ω 2 6 √ αxαy ω 3 ]. Hence, the increase is here still more rapid than at both previous borders. These critical behaviors are all depicted in Fig. 5.

D. Entanglement control
We finally show in Fig. 6 the possibilities offered by this model for controlling the entanglement growth through a stepwise time dependent frequency, starting from the separable ground state of H 0 . After applying a "low" initial frequency ω = 0.5ω x for ω x t < 30, which leads to a weak quasiperiodic entanglement, by tuning ω to a value close to the first instability ω y = k y for a finite time (30 < ω x t < 60), it is possible to achieve a large entanglement increase (curve a). Then, by setting ω = 0 (i.e., switching off the field or rotation), entanglement is kept high and constant, since the evolution operator becomes a product of local mode evolutions. Finally, by turning the frequency on again up to a low value, entanglement can be made to exhibit strong oscillations, practically vanishing at the minimum if ω is appropriately tuned. Thus, disentanglement at specific times can be achieved if desired. The entanglement increase at the second interval can be enhanced by allowing the system to enter the instability region for a short time, as shown in curve b. and β + are positive but α − β − < 0, so that the potential is stable in one direction but unstable in the other (i.e., α − > 0, β − < 0). Here λ + is real but λ − is imaginary.
The decomposition (41) no longer holds at the critical curve ∆ = 0, which separates sector E from sectors C and D. At this curve (including points L), the system is inseparable, in the sense that it cannot be written as a sum of two independent quadratic systems, even if allowing for complex coordinates and momenta as in sector E. While the matrix H 2 (Eq. (9)) is always diagonalizable for ∆ = 0, i.e., whenever the decomposition (41) is feasible, both H and H 2 are non-diagonalizable when ∆ = 0. Here λ ± = λ, with λ real at the border between B and E, imaginary between C and E and zero at the points L. The optimum decomposition of h in these cases is discussed in [20].