Quantum state preparation for coupled period tripling oscillators

TL;DR

This paper explores the quantum transition to correlated states in coupled oscillators exhibiting period tripling, highlighting the role of geometric phases and phase correlations in the quantum regime.

Contribution

It introduces a quantum framework for understanding period tripling in coupled oscillators, emphasizing the role of geometric phases and phase correlations.

Findings

Correlations form between oscillation phases during the transition.

Transient symmetry breaking occurs due to geometric phases.

Wigner distribution shows a minimum at zero amplitude in the quantum regime.

Abstract

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.

Figures

Equations

$$H=H_s+H_d+H_i$$

$$H_s=\sum _n\Delta a_n^{†}{a}_n+K(a_n^{†}{)}^2{a}_n^2,$$

$$H_d=-r\sum _n\left[a_n^3+(a_n^{†}{)}^3\right]$$

$$H_0=$$

$$-r(X^3-3YXY)/\sqrt{2},$$

$$H_i=-\sum _{m\ne n}^N{V}_{mn}{a}_m^{†}{a}_n.$$

$$E(\theta )$$

$$\mathbb{T}\{j_1{j}_2...j_N\}=\{j_N{j}_1...j_{N-1}\},\mathbb{R}\{j_1{j}_2...j_N\}=\{j_N{j}_{N-1}...j_1\},$$

$${\mathbb{N}}_3\{j_1...j_N\}=\{j_1-\mathrm{1...}{j}_N-1\}.$$

$${\partial }_{t}W=\left[i\Delta {\partial }_{\alpha }\alpha +2iK{\partial }_{\alpha }\alpha (|\alpha {|}^2-1)-i\frac{K}{2}{\partial }_{\alpha }^2{\partial }_{{\alpha }^{\ast }}\alpha \right]W$$

$$+\left[-ir\left(3{\partial }_{\alpha }({\alpha }^{\ast }{)}^2+\frac{1}{4}{\partial }_{\alpha }^3\right)+\frac{1}{2}\kappa ({\partial }_{\alpha }\alpha +\frac{1}{2}{\partial }_{\alpha }{\partial }_{{\alpha }^{\ast }})\right]W+c.c..$$

$$H_d=-r\sum _n\left[a_n^2+(a_n^{†}{)}^2\right],$$

Full text

Introduction.–

The adiabatic theorem in quantum mechanics Born and Fock (1928) 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 Kaminsky et al. (2004); Boixo et al. (2014, 2016). 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 Kadowaki and Nishimori (1998); Farhi et al. (2001); Childs et al. (2001).

The interest in adiabatic many-body dynamics has now extended to systems of quantum oscillators Goto (2016); Nigg et al. (2017); Puri et al. (2017a); Savona (2017); Mamaev et al. (2018); Goto et al. (2018); Zhao et al. (2018); Dykman et al. (2018); Goto (2019); Rota and Savona (2019). 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 Goto (2016); Zhang and Dykman (2017); Puri et al. (2017b); Goto et al. (2018); Wang et al. (2019), and how this can be used for adiabatic quantum computing with oscillator arrays Fitzpatrick et al. (2017); Ma et al. (2019). Coupled coherent parametrically driven oscillators can go through a quantum phase transition into a correlated state (a “time-crystal” effect with no disorder)  Dykman et al. (2018).

Parametric oscillators can be mapped  Goto (2016) onto an “Ising machine”, which has recently been demonstrated in the classical regime with 100-2000 optical spins McMahon et al. (2016); Inagaki et al. (2016).

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 Zhang and Dykman (2017) between the quantum states at the minima of the effective oscillator Hamiltonian in Fig. 1 first noticed in Ref. Guo et al., 2013. 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 one-dimensional 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 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 weakly-damped 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 setup and Hamiltonian.–

We study arrays of $N$ coupled driven oscillators. The Hamiltonian

[TABLE]

consists of the Hamiltonian of the undriven oscillators ($H_{s}$), the driving term ($H_{d}$), and the interaction ($H_{i}$). We assume that all oscillators are identical and have inversion symmetry, and we keep in $H_{s}$ the lowest-order intrinsic nonlinearity (called Duffing or Kerr nonlinearity). In the frame that rotates at 1/3 the drive frequency $\omega_{F}$ and in the familiar rotating wave approximation (RWA) Walls and Milburn (2008)

[TABLE]

where $a_{n}$ and $a_{n}^{\dagger}$ are the ladder operators of the $n$th oscillator. In Eq. (2) we introduced the detuning $\Delta=\omega_{0}-\omega_{F}/3$ of the drive with respect to the oscillator eigenfrequency $\omega_{0}$; $K$ is the nonlinearity parameter, and we set $\hbar=1$.

The Hamiltonian that describes the driving

[TABLE]

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. Zhang and Dykman (2017).

From Eqs. (2) and (3), we can write the RWA Hamiltonian $H_{0}=H_{s}+H_{d}$ of an individual oscillator as

[TABLE]

where $X$ and $Y$ correspond to the scaled coordinate and momentum, $X=(a^{\dagger}+a)/\sqrt{2}$ and $Y=i(a^{\dagger}-a)/\sqrt{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(\Delta-2K)/9$. Classically, they correspond to different phases $\theta=0,\,2\pi/3$, and $4\pi/3$ of the period-3 oscillations.

We assume linear coupling between the oscillators. After an RWA it is described by the interaction Hamiltonian

[TABLE]

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\delta_{m,n\pm 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\pi/3$, which is realized by the unitary operator ${\mathbb{N}}_{3}=\exp\left[-(2\pi i/3)\sum\nolimits_{n}a_{n}^{\dagger}a_{n}\right]$. The other symmetry operations are translation ${\mathbb{T}}^{\dagger}a_{n}{\mathbb{T}}=a_{n+1}$ and reversing ${\mathbb{R}}^{\dagger}a_{n}{\mathbb{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, $\hat{I}=\frac{1}{\pi}\int_{0}^{\infty}|\alpha|\,\mathrm{d}|\alpha|\int_{0}^{2\pi}\mathrm{d}\theta\ket{\alpha}\bra{\alpha}$, with $\alpha=|\alpha|\exp(i\theta)$, to define the measurement operators

[TABLE]

In terms of the oscillator Fock states $\ket{k}$ in the absence of driving, E(\theta)=\frac{1}{\pi}\sum_{k,k^{\prime}=0}^{\infty}\frac{\Gamma\bigl{(}(k+k^{\prime}+2)/2\bigr{)}}{\sqrt{k!k^{\prime}!}}\cdot\frac{\sin\left[(k-k^{\prime})\theta\right]}{k-k^{\prime}}\cdot\ket{k}\bra{k}^{\prime}, where $\Gamma(x)$ is the Gamma-function, and we use the convention $\theta>0$. The approximate effect of $E(\theta)$ is a projection on the sector of phase space bounded by the polar angles $-\theta$ and $\theta$. As the coherent states do not form an orthogonal basis, $E(\theta)$ 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) Busch et al. (1996).

We define $P_{0}=E(\pi/3)$, corresponding to the third of phase space limited by the polar angles $\pm\pi/3$. For one oscillator, where the phase-space rotation operator is $N_{3}=\exp[-(2\pi i/3)a^{\dagger}a]$, we define the rotated operators $P_{1}=N_{3}^{\dagger}P_{0}N_{3}$ and $P_{2}=N_{3}^{\dagger}P_{1}N_{3}$ corresponding to, respectively, the sectors rotated by $2\pi/3$ and $4\pi/3$. As $P_{0}+P_{1}+P_{2}=\hat{I}$, the $P$-operators form a POVM, and we define the corresponding probabilities as $p_{j}=\bra{\psi}P_{j}\ket{\psi}$, where $\ket{\psi}$ 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}=\bra{\psi}P_{j}\otimes P_{k}\ket{\psi}$. In the general case $p_{j_{1}...j_{N}}=\bra{\psi}\Pi_{n=1}^{N}P_{j_{n}}\ket{\psi}$.

Quasi-adiabatic state preparation.–

We will assume that each oscillator is initialized in the vacuum state, $\ket{\psi}_{\rm ini}=\prod\nolimits_{n=1}^{N}\ket{\rm vac}_{n}$, which is the ground state for $r=0$ if the initial detuning $\Delta_{\mathrm{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_{\rm f})r_{\mathrm{max}}$, where $t_{\rm f}$ is the ramp time. Simultaneously the detuning is linearly decreased to [math], i.e. $\Delta(t)=\left(1-t/t_{\rm f}\right)\Delta_{\mathrm{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 $\ket{\psi}_{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)$ Zhang and Dykman (2017). The states $\ket{\psi}_{j}$ correspond to classical period-3 oscillations with the phases that differ by $2\pi/3$ for different $j$. We associate $j=0,1$ and $2$ with the directions $0,2\pi/3$ and $4\pi/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^{\prime}}\approx\delta_{jj^{\prime}}$.

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\theta\sum_{n}a^{\dagger}_{n}a_{n})$ commutes with $H_{s}+H_{i}$. For $r=0$ and $|V|\ll\Delta_{\rm 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}...j_{N}}=(1/3)^{N}$.

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 $\ket{\psi}_{j_{1}}...\ket{\psi}_{j_{N}}$ of intrawell states $\ket{\psi}_{j_{n}}$ of individual oscillators. Our measurement directly reveals such combinations.

Symmetry arguments.–

The multi-oscillator initial ($r=0$) state $\ket{\psi}_{\rm ini}=\prod\nolimits_{n=1}^{N}\ket{\mathrm{vac}}_{n}$ provides the totally symmetric representation of the group generated by the operators ${\mathbb{N}}_{3}$, ${\mathbb{T}}$, and ${\mathbb{R}}$. Since the full Hamiltonian (1) is invariant under these symmetry operations, the state $\ket{\psi(t)}$ obtained by evolving $\ket{\psi}_{\rm 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 $\ket{\psi(t_{\rm 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 $\ket{\psi(t_{\rm f})}$ can be indeed close to the adiabatic state 111The results on the evolution of the symmetric state of two oscillators are given in Appendix A.

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 $\Delta=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 $\ket{\psi}_{j_{1}}...\ket{\psi}_{j_{N}}$ can be denoted as $\{j_{1}...j_{N}\}$, where $j_{n}$ refers to the $n$th oscillator. To first order, the coupling energy in such a state is $-V\sum_{m}(X_{j_{m}}X_{j_{m+1}}+Y_{j_{m}}Y_{j_{m+1}})$ where $(X_{j},Y_{j})$ is the position of the $j$th minimum of $H_{0}$ on the phase plane. The operators ${\mathbb{N}}_{3},{\mathbb{T}}$ and $\mathbb{R}$ can be thought of as shift operators in the space of $\{j_{1}...j_{N}\}$,

[TABLE]

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 $\mathbb{T}$, $\mathbb{N}$, and $\mathbb{R}$ to $\prod_{k}\ket{\psi}_{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 $\{jj...\}$ 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 (Symmetry arguments.–) 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 $\Delta$ 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 $\Delta$, similar to the case of a single oscillator Zhang and Dykman (2017). 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). Here, the transient populations of the would-be equivalent orientations $\{001\}$ and $\{002\}$ are different. The probability oscillations are more pronounced for antiferromagnetic coupling, where nonadiabatic effects are stronger.

Further insight into the features of the multi-oscillator states is provided by their energy spectra. The lowest eigenvalues of the Hamiltonian of a three-oscillator array for $t=t_{\rm f}$, are shown in Figs. 2 (c) and (d). Out of 27 states (combinations of three intrawell states of three oscillators) one can form four fully symmetric states. Three of them are occupied both for ferro- and antiferromagnetic coupling. For the fully adiabatic 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_{0}^{2}$ and $3|V|X_{0}^{2}/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 $\sim 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,\Delta=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^{\dagger}$ that couples the oscillator to a thermal reservoir. The dissipation-induced change of the density matrix $\partial_{t}\rho$ is described by the standard operator $\hat{L}\rho=\frac{1}{2}\kappa(2a\rho a^{\dagger}-a^{\dagger}a\rho-\rho a^{\dagger}a)$; here $\kappa$ is the energy decay rate, and we have set the oscillator Planck number $\bar{n}=0$.

The difference between the classical and quantum dynamics is most easily seen from the equation for the Wigner distribution $W(\alpha,\alpha^{*})$. It can be derived in a standard way Walls and Milburn (2008),

[TABLE]

Here, the terms with the first-order derivatives describe classical dynamics in the rotating frame in the absence of quantum fluctuations. For $(3r/2K)^{2}>\sqrt{(2-\Delta/K)^{2}+(\kappa/2K)^{2}}-2+\Delta/K$, the classical oscillator has three stable states with nonzero $|\alpha|$; they correspond to period-3 oscillations in the lab frame; the state at $\alpha=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 $|\alpha|$ varies on the scale $\gg 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 $|\alpha|\sim 1$ in the classical period-3 states, the third derivatives change the distribution qualitatively. The maximum at $\alpha=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\sim 1$. As for all parameters in this article, for $\alpha=0$ the eigenvalues of the Hessian had the same sign, we can use the sign of the Laplacian $\partial_{\alpha}\partial_{\alpha^{*}}W$ to distinguish whether $W$ is maximal or minimal at $\alpha=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 $\partial_{\alpha}\partial_{\alpha^{*}}W(0)$ for large $r$ results from the saddle points of $H_{0}$ approaching $\alpha=0$. Therefore quantum fluctuations become strong and wash away the classical stability of the state $\alpha=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 Fitzpatrick et al. (2017); Ma et al. (2019), as they combine strong enough nonlinearities and long coherence times. In a single cavity, period tripling has already been observed Svensson et al. (2017).

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.

Acknowledgments.–

We are grateful to Gil Refael and Mark Rudner for a valuable discussion. CB and NL acknowledge financial support by the Swiss SNSF and the NCCR Quantum Science and Technology. MID’s research was supported in part by the National Science Foundation (Grant No. DMR-1806473) and the Moore Scholarship from Caltech. Y.Z. was supported by the National Science Foundation (DMR-1609326). All the quantum simulations have been performed using QuTiP Johansson et al. (2013), the semiclassical partial differential equations were solved with numpy Oliphant (06) and scipy Jones et al. (01).

Appendix A A two-oscillator system.–

In this section we provide the results that complement the results on period tripling oscillators presented in the main text. Figure 5 shows sweeps for two oscillators that are analogous to the sweeps of three and four oscillators shown in Figs. 2 and  3 of the main article, with the same parameter values. For ferromagnetic coupling, as in the case of 3 and 4 oscillators, the system approaches one of the 3 equivalent configurations $\{jj\}$, which leads to the probability $p_{00}\approx 1/3$. For antiferromagnetic coupling, the most probable configurations are $\{j,j+1\}$, giving $p_{01}\approx 1/6$ at the end of the sweep.

Appendix B A strongly non-classical Wigner distribution of a dissipative oscillator.–

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 $\Delta$ and variable Planck number $\bar{n}$. On increasing $\bar{n}$, the region of exhibiting quantum behavior (green area where $\partial_{\alpha}\partial_{\alpha^{*}}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).

Appendix C Comparison to period doubling.–

For reference, we briefly discuss the case of period doubling, where the drive Hamiltonian in the rotating frame, given by Eq. (3) in the main text, is replaced by

[TABLE]

stemming from a parametric modulation at frequency $\omega_{F}$ close to twice the oscillator eigenfrequency $\omega_{0}$; in this case $\Delta=\omega_{0}-\omega_{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(\pi/2)$ on the right half plane and $P_{0}=\exp(i\pi a^{\dagger}a)P_{1}$ on the left half-plane. The probability $p_{j}$ for an oscillator in state $\ket{\psi}$ to be in bit $j$ is then $p_{j}=\bra{\psi}P_{j}\ket{\psi}$, 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 anti-ferromagnetic coupling and three oscillators, where it is at $1/6$. In the four-oscillator case, the ferromagnetic and anti-ferromagnetic 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=\beta^{2}t\sigma_{z}+\Omega\sigma_{x}$, where $\beta$ parameterizes the sweep rate and $2\Omega$ is the minimal energy gap, the transition probability to the higher-lying state at each of these crossings is then approximately given by $P_{\mathrm{LZ}}=1-\exp(-\pi\Omega^{2}/\beta^{2})$. In our setup, increasing $\Delta$ and $V$ leads to larger energy gaps. Together with the sweep time $t_{\rm f}$ they fully characterize a sweep. Conveniently, we can fix $r_{\mathrm{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 $(\Delta,V,t_{\rm 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_{\rm f}$ for the state evolution to remain adiabatic, panels (b) and (d) show it as a function of $V$ and $\Delta_{\mathrm{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 $\Delta_{\mathrm{ini}}$ increases for larger coupling $V$. For both plots, we conclude that the requirements are more demanding for period tripling as compared to doubling.