Quantum state preparation for coupled period tripling oscillators

We investigate the quantum transition to a correlated state of coupled oscillators in the regime where they display period tripling in response to a drive at triple the eigenfrequency. Correlations are formed between the discrete oscillation phases of individual oscillators. The evolution toward the ordered state is accompanied by the transient breaking of the symmetry between seemingly equivalent configurations. We attribute this to the nontrivial geometric phase that characterizes period tripling. We also show that the Wigner distribution of a single damped quantum oscillator can display a minimum at the classically stable zero-amplitude state.

Introduction.-The adiabatic theorem in quantum mechanics [1] states that a quantum system in the instantaneous ground state of a time-dependent Hamiltonian will approximately remain there if the Hamiltonian changes slowly compared to the gap to the first excited state. Recently the adiabatic dynamics in many-body systems has been extensively studied with arrays of qubits [2][3][4]. One promising application is adiabatic quantum computing, where the initial Hamiltonian is well-understood, so that initialization of its ground state is straightforward, and the final Hamiltonian encodes the cost function of an optimization problem that is hard to solve on a classical computer [5][6][7].
The interest in adiabatic many-body dynamics has now extended to systems of quantum oscillators [8][9][10][11][12][13][14][15][16][17]. This was triggered by the observation how, with turning on parametric driving close to twice the oscillator eigenfrequency, the ground state of a single oscillator adiabatically connects to the cat state [8,13,[18][19][20], and how this can be used for adiabatic quantum computing with oscillator arrays [21,22]. Coupled coherent parametrically driven oscillators can go through a quantum phase transition into a correlated state (a "timecrystal" effect with no disorder) [15].
Parametric oscillators can be mapped [8] onto an "Ising machine", which has recently been demonstrated in the classical regime with 100-2000 optical spins [23,24].
The many-body dynamics of driven coupled oscillators can be radically different if the driving frequency is close to triple the oscillator eigenfrequency. An isolated oscillator can display period tripling in this case. A particular feature of the effect is the geometric phase [18] between the quantum states at the minima of the effective oscillator Hamiltonian in Fig. 1 first noticed in Ref. 25. It can be thought of as resulting from a "magnetic field" that pierces the oscillator phase space.
In this paper we study how the geometric phase of the quantum period tripling and the high degeneracy of the period-3 states affect the dynamics of coupled quantum oscillators. Specifically, we study how the system goes into a coherent many-body state as the driving field is slowly turned on and tuned close to resonance. The results refer to a onedimensional oscillator array with either attractive or repulsive couplings. Such couplings favor, respectively, the same or different phases of the period-3 oscillations and are analogous to ferro-or antiferromagnetic coupling in the case of spins. The The plot corresponds to H 0 , Eq. (4), in units of the Kerr parameter K; X and Y are the scaled coordinate and momentum, r = 1.4K, and ∆ = 0. In the main text, the minima are enumerated counterclockwise as j = 0, 1, 2 starting with the minimum on the axis Y = 0. (b) Wigner distribution in the lowest fully symmetric eigenstate of H 0 for r = 1.4K and ∆ = 0. (c) Eigenvalues of H 0 as functions of r/K where ∆ = ∆ ini (1 − r/r max ), ∆ ini = 6K. For r = 0 the spectrum is that of a weakly anharmonic oscillator and the levels are color-coded as n = 3k (red), 3k + 1 (blue), 3k + 2 (green), with k = 0, 1, 2, .... With the increasing r/K the levels with different k merge into triples of tunnel-split intrawell levels of H 0 . case of antiferromagnetic coupling is particularly interesting because multiple configurations can lead to neighboring oscillators having different phases. We note that, because of the geometric phase, the oscillator chain cannot be simply mapped on a chain of spin-1 particles.
We also study the stationary distribution of a single weaklydamped oscillator in the ultra-quantum regime to explore whether period tripling can qualitatively change this distribution compared to what would be expected in the semiclassical limit. The very possibility of such a change is a consequence of the peculiar semiclassical dynamics where the unstable period-3 states approach the stable state with the increasing drive, but do not merge with this state. Physical where a n and a † n are the ladder operators of the nth oscillator. In Eq. (2) we introduced the detuning ∆ = ω 0 − ω F /3 of the drive with respect to the oscillator eigenfrequency ω 0 ; K is the nonlinearity parameter, and we set = 1.
The Hamiltonian that describes the driving corresponds to the energy of an oscillator in the driving field, which is proportional to the field multiplying the cube of the oscillator coordinate, with r being the scaled field amplitude. The term (3) can arise also from a coupling linear in the coordinate or momentum taking into account the oscillator nonlinearity, cf. [18]. From Eqs.
(2) and (3), we can write the RWA Hamiltonian H 0 = H s + H d of an individual oscillator as where X and Y correspond to the scaled coordinate and momentum, X = (a † + a)/ √ 2 and Y = i(a † − a)/ √ 2. The classical phase-space energy surface corresponding to H 0 is shown in Fig. 1 along with an example of the Wigner distribution. The Hamiltonian has a three-fold symmetry in the oscillator phase space, a feature of period tripling. The three minima away from X = Y = 0 emerge for r 2 > 8K(∆ − 2K)/9. Classically, they correspond to different phases θ = 0, 2π/3, and 4π/3 of the period-3 oscillations.
We assume linear coupling between the oscillators. After an RWA it is described by the interaction Hamiltonian To reveal the novel features of the many-body dynamics coming from period tripling, we consider the simplest model of the oscillator array: a nearest-neighbor coupling, V mn = Vδ m,n±1 , and periodic boundary conditions. For the "ferromagnetic" and "antiferromagnetic" cases, V > 0 and V < 0, respectively. Below we loosely use the term "energy" for the eigenvalues of the Hamiltonian H. The Hamiltonian is invariant under simultaneous rotation of all oscillators by −2π/3, which is realized by the unitary operator N 3 = exp −(2πi/3) n a † n a n . The other symmetry operations are translation T † a n T = a n+1 and reversing R † a n R = a N+1−n the order of the oscillators.
Measurement of states.-For the classification and measurement of the states we use the resolution of unity with coherent states,Î = 1 π ∞ 0 |α| d|α| 2π 0 dθ |α α|, with α = |α| exp(iθ), to define the measurement operators In terms of the oscillator Fock states |k in the absence of driv- is the Gamma-function, and we use the convention θ > 0. The approximate effect of E(θ) is a projection on the sector of phase space bounded by the polar angles −θ and θ. As the coherent states do not form an orthogonal basis, E(θ) is not a projector, but corresponds to a more general form of measurement that can be described in the framework of Positive Operator Valued Probability Measures (POVMs) [27].
We define P 0 = E(π/3), corresponding to the third of phase space limited by the polar angles ±π/3. For one oscillator, where the phase-space rotation operator is , we define the rotated operators P 1 = N † 3 P 0 N 3 and P 2 = N † 3 P 1 N 3 corresponding to, respectively, the sectors rotated by 2π/3 and 4π/3. As P 0 + P 1 + P 2 =Î, the P-operators form a POVM, and we define the corresponding probabilities as p j = ψ| P j |ψ , where |ψ is the oscillator wave function. These definitions naturally generalize to arrays of oscillators. For two oscillators the probability of the first oscillator to be in sector j and the second oscillator to be in sector k is p jk = ψ| P j ⊗ P k |ψ . In the general case p j 1 ... j N = ψ| Π N n=1 P j n |ψ . Quasi-adiabatic state preparation.-We will assume that each oscillator is initialized in the vacuum state, |ψ ini = N n=1 |vac n , which is the ground state for r = 0 if the initial detuning ∆ ini is positive and large compared to the coupling strength. We then ramp up the scaled driving amplitude r linearly to its maximal value, r(t) = (t/t f )r max , where t f is the ramp time. Simultaneously the detuning is linearly decreased to 0, i.e. ∆(t) = (1 − t/t f ) ∆ ini . All other parameters are kept constant. In the numerical plots, all energies and frequencies are in units of K.
We are interested in the state of the system at the end of the sweep. If the oscillators are uncoupled and the sweep is fully adiabatic, the state of each of them for not too small r will be a symmetric superposition of states |ψ j ( j = 0, 1, 2) localized on the phase plane (X, Y) in Fig. 1 at the minima of the Hamiltonian function H 0 (X, Y) [18]. The states |ψ j correspond to classical period-3 oscillations with the phases that differ by 2π/3 for different j. We associate j = 0, 1 and 2 with the directions 0, 2π/3 and 4π/3 on the phase plane toward the wells of H 0 , respectively, or equivalently, with the number of the well. If the oscillator is in the state j, the POVM measurement will give the probability p j ≈ δ j j .
Coupling the oscillators leads to correlations between their oscillation phases to minimize the coupling energy. Without the drive (r = 0), the energy of an individual oscillator is independent of its phase, whereas the multi-oscillator state is invariant only with respect to the continuous global phase, the rotation operator exp(−iθ n a † n a n ) commutes with H s + H i . For r = 0 and |V| ∆ ini the ground-state multi-oscillator wave function is the product of the ground-state wave functions of the individual oscillators, and then p j 1 ...
Not only does the drive break the continuous phase symmetry of an individual oscillator, but it also reduces the level spacing within the triples of its neighboring energy levels, see Fig. 1. Therefore the oscillator coupling becomes effectively stronger with increasing r and its effect becomes more pronounced. For large r, the low-energy multi-oscillator states are combinations of the products |ψ j 1 ... |ψ j N of intrawell states |ψ j n of individual oscillators. Our measurement directly reveals such combinations.
Symmetry arguments.-The multi-oscillator initial (r = 0) state |ψ ini = N n=1 |vac n provides the totally symmetric representation of the group generated by the operators N 3 , T, and R. Since the full Hamiltonian (1) is invariant under these symmetry operations, the state |ψ(t) obtained by evolving |ψ ini will remain totally symmetric. Such a state is not necessarily the ground state of the full Hamiltonian. However, it is the lowest-energy totally symmetric state. If the evolution is slow on the scale determined by the gaps between the totally symmetric states, the final state |ψ(t f ) will be the lowest-energy totally symmetric state.
In Figs. 2 and 3 we show, using our POVM-based measurement for a system with three and four oscillators, that |ψ(t f ) can be indeed close to the adiabatic state [28].
In our simulations the driving parameter r was ramped up to r max = 1.4K. As seen from Fig. 1 (c), for these values of r and ∆ = 0 the three lowest energies of a single oscillator are close to each other and turn into the tunnel-split energies of the linear combinations of the intrawell states of H 0 .
The products of weakly perturbed intrawell states of individual oscillators |ψ j 1 ... |ψ j N can be denoted as { j 1 ... j N }, where j n refers to the nth oscillator. To first order, the cou- The totally symmetric state of the coupled oscillators is found in a standard way by summing the wave functions obtained by repeatedly applying the operators T, N, and R to k |ψ j k . Configuration symmetry breaking in the transient regime.
-For the case of ferromagnetic interaction, the probability to find all oscillators aligned along one direction in the ground state, i.e., to be in the configuration { j j...} with j = 0, 1, or 2 for large r, is close to 1/3, independent of the number of oscillators. This probability is indeed approached in the sweep, as seen from the black lines in Figs. 2 (a) and 3 (a).
For anti-ferromagnetic interaction the situation is more interesting, as seen from Figs. 2 (b) and 3 (b). For three oscillators the configuration that minimizes the antiferromagnetic coupling energy for large r is { j 1 j 2 j 3 } with all j 1,2,3 being different from each other. There are six such configurations. The totally symmetric state can be obtained by applying successively the symmetry operators (7) to the configuration {012}. Respectively, for the adiabatic state preparation, the probability to find the system in one of the configurations will be 1/6. This is indeed seen in Fig. 2 (b).
For four oscillators, there are two configurations that both minimize the coupling energy for large r, to leading order in H i and neglecting tunneling. They are {0102} and {0101}, and the respective totally symmetric states built out of them. The only difference between the configurations {0102} and {0101} is that, in the first of them, oscillator 4 is in the well rotated clockwise with respect to the neighboring oscillators, whereas in the second, this oscillator is in the well rotated counterclockwise. The equivalence of the configurations is broken by the geometric phase between the intrawell states. The energy splitting between the corresponding totally symmetric states is small, leading to strong nonadiabaticity with varying r and ∆ and to a similar population of the states. In turn, this leads to the strong slow oscillations of the configuration populations in Fig. 3(b). The oscillation period increases as the energy splitting falls off. Which of the totally symmetric states has a lower energy depends on the values of r and ∆, similar to the case of a single oscillator [18]. There are no reasons to expect that the symmetric combination of these states has the lowest energy.
The effect of the geometric phase is seen also in Fig. 2 (a,b).  evolution only the lowest energy one (the black dot) will be occupied. To first order in the coupling, its energy is shifted down from the energy of noninteracting oscillators by 3|V|X 2 0 and 3|V|X 2 0 /2 for the ferro-and antiferromagnetic coupling respectively; here, X 0 is the distance of the phase-space minima of H 0 from the origin. These expressions are an overestimate by ∼ 30% for r = 1.4. The excited fully symmetric states (the red dots) are also partly occupied. In leading order, they have the same energy for ferro-and antiferromagnetic coupling.
As a test, we studied a frustrated triangle of oscillators, where the first and the third oscillators are coupled antiferromagnetically, but the second oscillator is coupled ferromagnetically to the other two. In the absence of the geometric phase, the configurations {000}, {011}, and {022} would be expected to have the same energy. However, we found that, for the same parameters r = 1.4K, ∆ = 0, and |V| = 0.4K, the symmetrized configuration {000} has the lowest energy.
Open period-3 system.-The peculiar features of the quantum-coherent dynamics of period-3 oscillations is expected to have a counterpart in the dissipative dynamics. Some aspects of this dynamics can be revealed by studying the stationary distribution of a dissipative oscillator in the ultra-quantum regime. We assume that the dissipation comes from a term linear in a, a † that couples the oscillator to a thermal reservoir. The dissipation-induced change of the density matrix ∂ t ρ is described by the standard operator Lρ = 1 2 κ(2aρa † − a † aρ − ρa † a); here κ is the energy decay rate, and we have set the oscillator Planck numbern = 0.
The difference between the classical and quantum dynamics is most easily seen from the equation for the Wigner distribution W(α, α * ). It can be derived in a standard way [26], Here, the terms with the first-order derivatives describe classical dynamics in the rotating frame in the absence of quan-tum fluctuations. For (3r/2K) 2 > (2 − ∆/K) 2 + (κ/2K) 2 − 2 + ∆/K, the classical oscillator has three stable states with nonzero |α|; they correspond to period-3 oscillations in the lab frame; the state at α = 0 is also stable. Within the classical theory one expects the stationary probability density to display peaks at the stable states. A classical description refers to the case where |α| varies on the scale 1 and corresponds, in particular, to disregarding the terms with the third derivatives in Eq. (8). However, we found that in the ultra-quantum regime, where |α| ∼ 1 in the classical period-3 states, the third derivatives change the distribution qualitatively. The maximum at α = 0 can turn into a minimum, see Fig. 4 (a) and (b). The minimum emerges once the drive becomes sufficiently strong and is most pronounced for r/K ∼ 1. As for all parameters in this article, for α = 0 the eigenvalues of the Hessian had the same sign, we can use the sign of the Laplacian ∂ α ∂ α * W to distinguish whether W is maximal or minimal at α = 0.
The local minimum of W at the origin disappears for larger frequency detuning, higher decay rate, or higher temperature, where quantum effects are less pronounced, see Appendix B, The small curvature ∂ α ∂ α * W(0) for large r results from the saddle points of H 0 approaching α = 0. Therefore quantum fluctuations become strong and wash away the classical stability of the state α = 0.
Conclusions.-As seen from the above analysis, for period tripling, driven coupled oscillators exhibit a quantum transition to a correlated state that is qualitatively different from the classical transition. For ferromagnetic coupling, with a slowly increasing drive, a quantum system adiabatically goes into a correlated state of period-3 oscillations. In contrast, a classical system will stay in the zero-amplitude state. Interestingly, the probabilities of different seemingly equivalent transient quantum configurations are different, hinting at an effect of the intrinsic geometric phase of the oscillators. These unusual features of quantum oscillator arrays can be studied with coupled nanomechanical resonators and optical cavities. A particularly promising platform is provided by coupled circuit-QED microwave cavities [21,22], as they combine strong enough nonlinearities and long coherence times. In a single cavity, period tripling has already been observed [29].
Another unexpected feature of period tripling is that, in the presence of dissipation, the stationary distributions of the quantum and classical oscillators are qualitatively different. In a certain parameter range, the quantum Wigner probability distribution displays a local minimum, rather than a maximum at the classically stable zero-amplitude state.
Our results show that period tripling in quantum oscillators allows studying new many-body phenomena far from thermal equilibrium, which have no analog in classical systems and in equilibrium quantum systems.
In this section we show more detailed results on the region where the Wigner distribution has a minimum at the classically stable state of zero vibration amplitude. As explained in the main text, in the parameter range we have explored, the difference between the maximum (classical regime) and the minimum (quantum regime) is given by the sign of the Laplacian of the steady state Wigner distribution at the origin in the oscillator phase space. Figure 6 shows scans of the Laplacian for variable detuning ∆ and variable Planck numbern. On increasingn, the region of exhibiting quantum behavior (green area where ∂ α ∂ α * W < 0) shrinks, as expected, since the oscillator becomes more "classical". On increasing the frequency detuning, this area shifts toward larger field amplitudes. Figure 7 illustrates the Wigner density for the parameters marked at the left dashed line of Fig. 4 (b).   Figure 7. Steady-state Wigner density of a dissipative oscillator in the period-tripling regime for the parameters used in Fig. 4 (b), i.e., κ = 0.5K,n = 0, ∆ = 0 and r = 0.75K. (a) Solution of the full master equation corresponding to Eq. (8). (b) Results of the semiclassical approximation obtained by solving a Fokker-Planck equation. As can be seen in Fig. 4 (b), the Laplacian in the center is zero for the quantum case, while the classical result show three maxima at the period-3 states and a maximum at the origin where the oscillator amplitude is zero. by Eq. (3) in the main text, is replaced by stemming from a parametric modulation at frequency ω F close to twice the oscillator eigenfrequency ω 0 ; in this case ∆ = ω 0 − ω F /2 in Eq.
(2) of the main text. We map the states of the oscillator to a bit using the measurement operators introduced in Eq. (6) with P 1 = E(π/2) on the right half plane and P 0 = exp(iπa † a)P 1 on the left halfplane. The probability p j for an oscillator in state |ψ to be in bit j is then p j = ψ| P j |ψ , where p 0 + p 1 = 1 as expected for the P-operators that form a POVM.
For period doubling the parameters are in a regime close to the adiabatic limit, so that the maximal probabilities are almost reached. This maximum occurs at 1/2, except for antiferromagnetic coupling and three oscillators, where it is at 1/6. In the four-oscillator case, the ferromagnetic and antiferromagnetic coupling are equivalent up to a basis transformation, therefore the curves in panels (c) and (d) of Fig. 8 agree.
The probability for the system to remain in the lowest fully symmetric state state or to switch to higher-lying fully symmetric states crucially depends on the rate of change of the system parameters and the energy gap to the excited states. In a simplified picture the system dynamics can be understood as a series of Landau-Zener transitions occurring at each avoided crossing the system passes through. For the Landau-Zener Hamiltonian H = β 2 tσ z + Ωσ x , where β parameterizes the sweep rate and 2Ω is the minimal energy gap, the transition probability to the higher-lying state at each of these crossings is then approximately given by P LZ = 1 − exp(−πΩ 2 /β 2 ). In our setup, increasing ∆ and V leads to larger energy gaps.
Together with the sweep time t f they fully characterize a sweep. Conveniently, we can fix r max , as the important transitions only occur during the phase transition, but not far above threshold, where the correlations between the oscillators are already effectively frozen. Also, we already have fixed K = 1 as a reference, as all other units are given in units of K. It is therefore sufficient to scan the parameters (∆, V, t f ). Due to the numerical cost we focus on the case of three oscillators. As a function of these parameters, Fig. 9 shows the probabilities of the configuration minimizing H i for three attractively coupled oscillators. For reference, the corresponding probabilities for period doubling are also shown. The left panels refer to period doubling and the right panels refer to period tripling. Note that for a perfectly adiabatic sweep, due to, respectively, the double and triple degeneracy, the maximal probability to be reached are 1/2 and 1/3. Figures 9 (a) and (b) show the threshold of coupling V and sweep time t f for the state evolution to remain adiabatic, panels (b) and (d) show it as a function of V and ∆ ini . While it is always better to sweep more slowly, there is a trade-off for the choice of the initial detuning: The optimal value for ∆ ini increases for larger coupling V. For both plots, we conclude that the requirements are more demanding for period tripling as compared to doubling.