Metastable quantum entrainment

Quantum van der Pol oscillators are driven-dissipative systems displaying quantum synchronization phenomena. When forced by a squeezed drive, the frequency adjusts to half of the forcing displaying multiple preferred phases. Here we analyze the physical origin of this entrained response, establishing a connection with metastability in open quantum systems. We report a dynamical regime characterized by a huge separation of time scales, in which a dynamical mode displays a lifetime that can be orders of magnitude larger than the rest. In this regime, the long-time dynamics is captured by an incoherent process between two metastable states, which correspond to the preferred phases of the synchronized oscillator. In fact, we show that quantum entrainment is here characterized by fluctuations driving an incoherent process between two metastable phases, which ultimately limits its temporal coherence when moving into the quantum regime. Finally, we discuss connections with the phenomena of dissipative phase transitions and transient synchronization in open quantum systems.


I. INTRODUCTION
The investigation of the properties of the nonequilibrium dynamics of dissipative quantum systems has recently become a major research subject, with focus ranging from their dynamics and control, to phase transitions, collective phenomena or thermodynamic features [10,14,18,19,28,54,62,69]. One of the consequences of the interplay between driving, dissipation and interactions is the possibility of having the dynamics of an open quantum system dominated by a huge separation of time scales, a characteristic of metastability [44]. In this case, initially highly excited configurations quickly relax to a set of metastable states that act as attractors of the dynamics for a long intermediate timescale, until final relaxation to the true stationary state occurs. This multi-step dynamics has different physical origins, as population trapping in a driven three-level system [44], the coexistence of different phases in an open Ising model [56] or in a non-linear oscillator [32], a dissipative discrete-time crystal [21], the long-time response of a driven-dissipative Rabi model [39], or a quasistationary regime in a chain of interacting fermions subject to a localized particle loss [70]. A general approach towards metastability in open quantum systems has been recently introduced in Ref. [44], in which some of its crucial signatures in the spectral properties of the Liouvillian have been identified.
The characteristic timescale separation of metastability has also been found in between different dynamical regimes of driven-dissipative systems [11,64], or as a precursor of eventual spontaneous symmetry breaking in the thermodynamic or infinite-excitation limits of these systems [51]. Indeed, sudden changes can occur in the stationary state of these systems in these limits, what are known as dissipative phase transitions (DPTs) [37,51]. As in the case of quantum phase transitions, DPTs are associated with a closure of the spectral gap of the generator of the dynamics [37,51]. This asymptotic gap closure has profound consequences even far from the thermodynamic or classical limits, as it can lead to a very slow relaxation timescale and hence metastability in the quantum dynamics [44]. While for first order DPTs this usually occurs in a narrow region between the different regimes [11,12], as experimentally observed in Refs. [18,19], for symmetry-breaking DPTs this occurs in a whole dynamical regime characterized by the emergence of metastable symmetry-broken states [3,50,51].
A paradigmatic phenomenon of classical driven-dissipative systems is synchronization, which emerges in widely different contexts and forms [2,55]. Entrainment is a type of synchronization phenomenon in which a self-oscillating system [61] adjusts its frequency and phase to that of an external forcing or to a multiple/fraction of it. Otherwise, synchronization can also occur as a mutual or spontaneous phenomenon among interacting system components. In the past decade, intense research activity has addressed the fundamental question of whether synchronization can emerge also in quantum systems [20]. As for its classical counterpart, diverse approaches and systems have been considered, since the seminal works of Refs. [27,72]. Different signatures and indicators of quantum synchronization have been explored in the literature [20,25,30,31,38,41,42,47,65], while its emergence has been reported in disparate scenarios as systems of spins [23,30,57], networks of quantum harmonic oscillators [6,45], optomechanical systems [43,67], atomic lattices [7,63] and clouds [71,73], and micromasers [13]. Moreover, novel forms of synchronization have been identified [5,40,42,63], as the case of transient synchronization [24], in which spontaneous synchronization emerges due to a long-lived collective excitation of the quantum system, a phenomenon that has been linked to the presence of super/subradiance [4] and long-lived correlations [6,23,25,45,46,59].
A rich playground where synchronization can be extended from the classical to the quantum regime is offered by the paradigmatic van der Pol oscillator (vdP) [2,55], introduced in Refs. [41,65]. For this model, both entrainment by a harmonic drive and spontaneous synchronization between coupled oscillators have been found in the quantum regime [40,41,65,66]. Signatures of synchronization have been identified in phase-space representations of the stationary state as well as in the power spectrum. Indeed, the stationary Wigner distribution [9] of the quantum van der Pol (QvdP) displays a ring-like shape in the absence of entrainment, which turns to a localized lobe for large enough driving strength [41,65]. Moreover, the frequency of the maximum of the power spectrum is reported to shift towards that of the driving as the system becomes entrained [65]. An important caveat is that the more in the quantum regime the system operates, the harder it is to synchronize it due to the detrimental effect of quantum fluctuations [65]. As an alternative implementation strategy, a squeezed forcing has been reported to enhance entrainment in the quantum regime [60]. However, very deep into it, where the system behaves effectively as a few-level system, the enhancing effect of squeezing is limited [53].
The introduction of the squeezed forcing comes at the expense of changing qualitatively the dynamical scenario of entrainment [35,36,60]. This is an observation of utmost importance to understand the entrained quantum dynamics of this system, as we show here. While in the classical harmonically driven vdP oscillator entrainment corresponds to a fixed point attractor (in the rotating frame) [55,68], here it corresponds to bistability (in the rotating frame) [35,60], characterized by frequency locking to half of the frequency of the forcing and two possible locked-phases. As it turns out, the nature of the classical attractors has a deep influence in the quantum dynamics. In fact, for the driven QvdP, the entrained dynamics has been shown to be well described by a linearized model around the fixed point attractor [68]. This approach sheds light on the synchronized dynamics of this system and on the reported behavior of the power spectrum, while it unveils new dynamical features as phase-coherent dynamics [68]. In contrast, we show that the squeezed QvdP displays a distinct dynamical scenario in which the main features of the entrained response cannot be captured by a linearized model, a fact rooted in the emergence of bistability in the classical limit.
Motivated by the rich physics behind the squeezed vdP oscillator, we show in this work that such a model offers the opportunity to establish a connection between metastability and entrainment, that we argue to transcend this specific system and to be applicable to other synchronization scenarios in the presence of phase multistability. More specifically, we identify the huge timescale separation characteristic of a metastable dynamics as the distinctive feature of the entrained dynamics for this system. On the shortest time scale, any initial state rapidly relaxes to the manifold of states spanned by two metastable states. These are the two preferred phases of the oscillator when it is well entrained by the forcing. On the longest time scale, the dynamics is dominated by a slow fluctuation mode connecting both phases, which eventually drives the state of the system to an even incoherent mixture of both. In between, there is a long period of time in which the response of the system is a quasi-coherent subharmonic oscillation. This is precisely how subharmonic entrainment manifests in the quantum regime, and it dominates the behavior of the power spectrum. While these phases act as effective attractors of the dynamics (any initial condition quickly relaxes to them), quantum fluctuations turn phase bistability to phase metastability, thus limiting the temporal coherence of the synchronized response on the long-time.

II. THE MODEL
In the laboratory frame, the squeezed van der Pol oscillator is described by the time-dependent master equation for the state of the systemρ l ( ̵ h = 1) [60]: where we have defined the dissipator D[L]ρ = 2LρL † −L †Lρ −ρL †L ,L is the corresponding jump operator,â is the annihilation operator of the bosonic mode, γ 1 is the amplification rate, γ 2 the nonlinear damping rate, η is the squeezing strength of frequency 2ω s , and ω 0 is the intrinsic frequency of the system. The explicit time dependence of this model enters only through the Hamiltonian term, which can be eliminated in a rotating frame, as can be achieved by a standard time-dependent unitary transformation,Û t = exp(−iω sâ †â t). This leads to a master equation with the same dissipative part but with a time-independent Hamiltonian 1 :Ĥ = ∆â †â + iη(â 2 −â †2 ) with ∆ = ω 0 − ω s . Notice that we denote the state of the system in this frame asρ, without any subscript, while the laboratory frame is indicated by the subscript 'L'. Possible experimental implementations of the QvdP oscillator have been described in the literature, considering platforms of trapped ions [41,60] and optomechanical oscillators [60,65,66]. In both cases, the QvdP oscillator corresponds to the properly engineered effective dynamics of a mechanical degree of freedom. Then, amplification and non-linear dissipation can be implemented by driving the fast optical degrees of freedom with lasers resonant with one-phonon absorption and two-phonon emission processes. Moreover, several possibilities for the implementation of the squeezed (two-boson) drive of Eq. (1) have been presented in [60] too, considering both optomechanical systems (e.g. modulating electrically the spring constant of the mechanical mode [58]) and trapped ions (e.g. by a combination of two Raman beams detuned by 2ω s [49]).
In this rotating frame, the master equation can be written in terms of the time-independent Liouvillian superoperator, defined by ∂ tρ = Lρ. Then, the dynamics of the system can be analyzed considering the eigenspectrum of this superoperator, which is composed by a set of eigenvalues λ j , and the corresponding set of right and left eigenmatrices defined by Lρ j = λ jρj and L †σ j = λ * jσ j , respectively [1,51] (more details in A). Generally, we can use the Liouvillian eigenmodes to decompose the state of the system at any timê Similarly, the eigenmode decomposition gives insight in the dynamics of expected values as well as multitime correlations, as we show below. In this expression we have defined the stationary state asρ ss =ρ 0 Tr[ρ 0 ], being the corresponding zero eigenvalue λ 0 = 0. Notice that the real part of these eigenvalues are non-positive, and we order them such that Re[λ 0 ] ≥ Re[λ 1 ] ≥ Re[λ 2 ] ≥ . . . . In the following, we will refer to the decay rates and frequencies of the eigenmodes respectively. While the formalism is fully quantum, we can identify the imbalance between nonlinear dissipation and linear amplification γ 2 γ 1 as the physical quantity controlling the excitation strength (i.e. boson number) in the family of QvdP systems [41,65], and hence the classical versus quantum system operation. In the limit γ 2 γ 1 ≫ 1 the QvdP oscillator is well approximated by a two-level system [40,53,65], while in the limit γ 2 γ 1 ≪ 1 the typical number of excitations becomes macroscopic. Indeed, in Ref. [8] we analyze in detail how the classical attractors emerge in this model as γ 2 γ 1 → 0. Here instead, we will focus on intermediate values of γ 2 γ 1 for which entrainment in the quantum regime has been reported [60].
A brief overview of the classical dynamics is presented in B, in which we show the system to display two different dynamical regimes depending solely on the relation between squeezing and detuning, η and ∆. In the rotating frame, for η < η c the stable attractor is a limit-cycle, while for η > η c the stable attractors are two fixed points. At the classical critical point η c = ∆ 2 the bifurcation between the two regimes occurs. In the laboratory frame the limit-cycle regime corresponds to the lack of entrainment, as generally its frequency is not commensurate with that of the forcing. On the other hand, the bistable fixed points oscillate harmonically at half of the frequency of the forcing, and thus this regime corresponds to subharmonic entrainment. Then, each of the fixed points corresponds to one of the possible bistable locked-phases, with a characteristic difference of phase of π.

III. TIMESCALE SEPARATION AND METASTABILITY
In this section we report on the signatures of metastability in the framework of the Liouvillian analysis: the presence of an eigenmode with a lifetime much longer than the rest, i.e. the opening of a spectral gap, and how this can be exploited to derive a simplified effective picture for the long-time dynamics of the system.

A. Opening of a spectral gap
The most important characteristic of the Liouvillian spectrum of the squeezed QvdP oscillator is that there is a whole parameter region where the minimum decay rate Γ 1 can be orders of magnitude smaller than Γ j≥2 . The inverse of such a small decay rate determines the longest relaxation timescale of our system. In Fig. 1 we demonstrate and analyze the formation of such a gap between the decay rates, while in the forthcoming sections we discuss its dynamical consequences and its relation to quantum entrainment.
The spectral gap is characterized through the ratio Γ 1 Γ 2 and can vary several orders of magnitude depending on the squeezing strength and the non-linear damping, Fig. 1 (a) and (b). As commented, we focus on parameter regimes for which the QvdP oscillator displays a moderate boson number, being neither strongly confined to the two lowest levels nor in the classical limit of large population. As the excitation number grows not only with the linear amplification but also with the squeezing strength, this corresponds to intermediate values of γ 2 γ 1 and small squeezing η γ 1 ( Fig. 1 (a)) or to larger non-linear damping while increasing the squeezing ( Fig. 1 (b)).
The ratio of decay rates displays the following characteristic dependence on the squeezing strength η: below a certain threshold value, this ratio remains constant and equal to one (bright region), while further increasing η the ratio diminishes steeply (dark region). This threshold value corresponds to an exceptional point (EP) of the Liouvillian [29,52], at which λ 1,2 become the same ( Fig. 1 (c-d)) while the corresponding eigenmatrices coalesce (not shown). The green-dotted line in Fig. 1 (a-b) indicates the squeezing strength at which the EP occurs (η EP ) as γ 2 γ 1 is varied. For η < η EP , λ 1,2 are complex conjugates, hence explaining the fact that Γ 1 Γ 2 = 1 for this region. While for η > η EP , λ 1,2 become real valued and Γ 1 decreases when increasing the squeezing strength, while Γ 2 increases, as exemplified in Fig. 1 (d).
Although the behavior of Γ 1 Γ 2 is qualitatively the same for all the considered range of γ 2 γ 1 , we find that there are important quantitative differences. In particular, the smaller is γ 2 γ 1 the wider is the gap for a given squeezing strength. This is intimately related to the behavior of Γ 1 as the classical limit is approached (i.e. γ 2 γ 1 → 0). Indeed, in this limit it is found that Γ 1 → 0 for η > η c while Γ 2 saturates to a non-zero value [8]. This means that the steady state becomes degenerate and the system undergoes a DPT [37,51]. This is intimately related to the emergence of the classical bistable attractors as we explain in detail in [8]. Notice that in this limit it is also found that η EP → η c . , eigenfrequencies and decay rates of the first two eigenmodes varying η γ1, respectively, and for γ2 γ1 = 0.1. In all panels ∆ γ1 = 0.1. Notice that in all this work the numerical results are obtained after truncation of the Fock space to a large enough Fock state such that the convergence of the results with this truncation can be assured. The needed truncation size generally increases as γ2 γ1 is diminished and as η γ1 is increased, and for the parameter regimes explored in this work this is bounded to the first 50 levels.
In broad terms, the spectral gap is enhanced in presence of a large number of excitations: the larger is the non-linear damping, the more squeezing is needed. The dynamics of the system displays a significant timescales separation for rates Γ 1 and Γ 2 differing by an order of magnitude or more, as in the (right) regions delimited by a red-dashed line (where Γ 1 Γ 2 ≤ 0.1) in Fig. 1 (a) and (b).

B. Effective long-time dynamics
The disparity between the weakest damping rate Γ 1 and the others allows to capture the system dynamics after a transient considering only the disparate timescales τ −1 1,2 = Γ 1,2 . This large separation of time scales makes an initially excited state to rapidly decay (∼ τ 2 ) to the manifold spanned byρ ss and the longest-lived eigenmodeρ 1 , where it displays a slow relaxation toρ ss at times of the order of τ 1 [44,56]. In the presence of a significant spectral gap, as reported in the previous section, there is a well-defined intermediate time scale, τ 2 ≪ t ≪ τ 1 , in which the contribution of the higher modes is negligible while the decay ofρ 1 is not yet appreciable, and thus the state of the system appears stationary in this time window. This is actually a signature of metastability [44], as we develop further in this section.
After times of the order of τ 2 , the dynamics is well-approximated just considering the contributions of the longestlived eigenmode and of the stationary state. As we have shown in the previous section, the parameter regime in which Γ 1 Γ 2 ≪ 1 corresponds to η > η EP , and thus λ 1,2 are real valued while the corresponding right and left eigenmatrices are Hermitian [1,51]. Then the long-time dynamics can be approximated aŝ which follows from Eq. (2) neglecting the contributions for the modes with j ≥ 2. By using the formalism introduced in Refs. [44,56], the approximate long-time dynamics of Eq. (4) can be recasted in terms of an effective stochastic process between two particular metastable states, which provides an insightful representation of this long-time response in a basis different from that provided by the stationary state and the longest-lived eigenmodeρ 1 , the latter not being a physical state since it is traceless [51]. In the following, we discuss the application of such formalism [44,56] to our system (further details are in C). From the long-time approximation of Eq. (4) it follows that the state of the system is restricted to the convex manifold of states spanned by the projection of the initial condition over the stationary state and the longest-lived eigenmode: We denote this manifold as the metastable manifold, as it captures the state of the system for the intermediate time scale τ 2 ≪ t ≪ τ 1 . In Refs. [44,56], it is shown that the metastable manifold can be parametrized in terms of the extreme metastable states (EMSs), defined as: where c max and c min are the maximum and minimum eigenvalues ofσ 1 . In our system, numerical analysis reveals that c min = −c max for the considered regime, i.e. when Γ 1 Γ 2 < 1. Moreover, we find that for a significant spectral gap (Γ 1 Γ 2 ≪ 1) these coefficients can be well approximated by c max = −c min ≈ 1, which yields the simpler expressions for the EMSs:μ 1(2) ≈ρ ss + (−)ρ 1 (see C). The projection of the state of the system at a given time onto the metastable manifold can be written in terms of the EMSs as: where the real coefficients p 1,2 (t) are defined in appendix C, while they are constraint to satisfy [44,56]: By means of the EMSs the long-time dynamics can be parametrized in terms of two physical states 2 , whose properties can be readily characterized. Crucially, the EMSs constitute the only basis for which p 1,2 (t) satisfy Eq. (8) in the whole metastable manifold. Indeed, in this case p 1,2 (t) can be interpreted as probabilities and the dynamics within the metastable manifold is given by [44,56]: whose solution makes Eq. (7) equivalent to Eq. (4). Thus, the approximate long-time dynamics for the state of the system has been recasted in the form of a two-state stochastic process between the EMSs, with a switching rate Γ 1 2 [22]. This characterizes the slow dynamics to which an initially excited state rapidly relaxes. From the solution of this dynamics (C) we find that p 1 (t → ∞) = p 2 (t → ∞) = 1 2, which could be also deduced from the fact that ρ ss = (μ 1 +μ 2 ) 2. This illustrates the notion of metastability in our system: any initial unbalanced mixture of the EMSs decays on a very long-time scale to the balanced one (ρ ss ), remaining apparently stable for the well-defined intermediate time scale τ 2 ≪ t ≪ τ 1 . In the following section we explore in detail the physical meaning of this statement by characterizing the properties ofμ 1,2 , through a visual representation based on the Wigner function, and of this effective incoherent dynamics between them.

IV. METASTABLE PREFERRED PHASES
As we show here, the characterization of the EMSsμ 1,2 provides a first important insight on the relation between the metastable dynamics and entrainment. In fact, we can gain intuition from their visual rendering as provided by their Wigner representation. In Fig. 2 we plot the Wigner distribution forρ ss , andμ 1,2 . As illustrated in Fig. 2 (a) and (d) the Wigner distribution of the stationary state displays a manifest bimodal character. Here we consider two disparate parameter values γ 2 γ 1 = (0.1, 3) and η γ 1 = (0.2, 2), both satisfying Γ 1 Γ 2 ≪ 1 as can be checked from Fig.  1. In fact, we find this bimodality to be present in the whole metastable regime. This can be intuitively understood by considering the Wigner distribution ofμ 1,2 , as shown in panels (b), (c) and (e), (f). We can appreciate that each of the lobes corresponds indeed to one of the metastable states.
The bimodality of the stationary state is a reminiscence of classical bistability. However, in stark contrast to the classical case, according to the quantum formalism these states are metastable, and the effective dynamics of Eq. (9) tells us that quantum fluctuations eventually drive the system to an even mixture of both. In fact, if we consider an initial condition consisting only of one lobe, i.e. p j (0) = 1 and p k (0) = 0 with j ≠ k, we can see how on times of the order of τ 1 the population of both lobes becomes progressively the same, until the stationary state is finally reached.
In fact, the Wigner distribution of the stationary state is a well known indicator of the emergence of phase-locking between the oscillator and the external forcing [41,42,60,65], or between coupled oscillators [40,43,66]. Indeed, for the driven QvdP oscillator, a sufficient strong forcing transforms the stationary Wigner distribution from a ring-like shape, indicating lack of a preferred phase, to a localized lobe, indicating the emergence of a preferred phase [41,65]. The latter corresponds to phase-locking in the quantum regime, in which quantum fluctuations play a non-negligible role [68]. Moving to the case of the squeezed QvdP oscillator, the emergence of phase preference is accompanied by the bimodality of the Wigner distribution [60] as displayed in Fig. 2. Therefore, each of the EMSsμ 1,2 can be interpreted as one of the two possible preferred phases to which the QvdP oscillator settles for large enough squeezing strength. It follows that the incoherent process between these two metastable preferred phases limits phase-locking in the long-time limit, as the stationary state that is reached is an even incoherent mixture of both. The manifestation of this incoherent process limiting phase-locking in different dynamical quantities is explored in more detail in the following section.
Beyond their phase-space representation, the properties of the metastable states can be understood from the fact that they are a linear combination ofρ ss andρ 1 with c max = −c min , and from the symmetries of the Liouvillian. Before proceeding further, notice that our master equation is invariant under the transformationâ → −â,â † → −â † , i.e. it displays parity-symmetry [1,51]. This means that we can define the unitary transformation Z 2ρ = e iπâ †âρ e −iπâ †â , whose action commutes with that of the Liouvillian [Z 2 , L]ρ = 0. Hence, the eigenmodes are either parity symmetric or parity antisymmetric, which means Z 2ρj = z jρj with z j = ±1, respectively [1,51]. An important observation is that the stationary state and the longest-lived eigenmode belong to different symmetry sectors: While the stationary state must be parity symmetric [1,51], the symmetry ofρ 1 is assessed numerically by realizing that its contribution for parity symmetric observables vanishes identically, i.e. Tr[ρ 1 (â † ) mân ] = 0 for m + n = even, while its contribution for parity antisymmetric ones is generally non-zero, i.e. Tr[ρ 1 (â † ) mân ] ≠ 0 for m + n = odd, where m and n are integers. As a consequence of this and of c max = −c min , it follows that the EMSs do not have a well defined parity-symmetry, and actually they are parity-broken states satisfying: Thus, both EMSs yield the same value for the expected value of parity symmetric observables, while for parity antisymmetric ones their value solely differs by a phase e iπ . Then, parity-symmetric observables are insensitive to this metastable dynamics, as for any Pρ(t) [see Eq. (7)] bothμ 1,2 contribute exactly the same, making irrelevant the particular evolution of p 1,2 (t) as they sum is one for all t. In contrast, notice the particular case of the amplitude 3 in which since it is an antisymmetric observable we have ⟨â⟩ 1 = −⟨â⟩ 2 . This is precisely the same phase relation between the two possible classical amplitudes in the synchronized regime, as briefly commented in Sec. II. In this sense, it turns out that, as the classical limit is approached, observables computed over the EMSs approach the corresponding classical values for each bistable fixed point [8]. These observations together with the Wigner representations of Fig. 2, suggest an intimate connection between the emergence of preferred phases in the quantum regime, and thus synchronization, and metastability.

V. METASTABLE ENTRAINED DYNAMICS
In this section we consider how the characteristic metastable response of the system manifests in the amplitude dynamics and two-time correlations, in order to characterize entertainment. We find that the metastable regime corresponds to the regime in which the system is entrained by the external forcing. Still, the temporal coherence of this entrained response is ultimately limited by the incoherent process between the two possible preferred phases. Finally, we find that while the long-time amplitude dynamics can be understood from a classical stochastic process, the quantum nature of the fluctuations can manifest in the two-time correlations.

A. Amplitude dynamics
After a short transient of the order of τ 2 , the amplitude dynamics for an arbitrary initial state is restricted to the metastable manifold and is well approximated by where p 1,2 (t) are the solutions of Eq. (9). Notice that the expectation value of the amplitude is defined as ⟨â(t)⟩ = Tr[âρ(t)], which can be written in terms of the Liouvillian eigenmodes as detailed in A. Then, within the approximation of the long-time transient given in Eq. (4), and using the EMSs decomposition of the metastable manifold, we can obtain the expression given in Eq. (12). The accuracy of Eq. (12) is numerically confirmed in Fig.  3 (a) for three different values of γ 2 γ 1 . In this panel we plot in logarithmic scale the imaginary part of ⟨â(t)⟩ (the real part would be similar) comparing the exact results of the full model (color solid lines) with the ones following the reduced EMSs effective long-time dynamics (black broken lines). We see that the initial excited (coherent) state rapidly relaxes to the metastable manifold which is accurately described by Eq. (12). Notice that Eq. (12) is the same we would obtain from the two-state stochastic process using the rules of classical stochastics [22], which state that the average dynamics of a first moment is the sum of the contributions of each state weighted by p 1,2 (t) (see also D). The behavior observed in Fig. 3 (a) can be fully understood recalling the Liouvillian analysis of Sec. III, where we found that whenever Γ 1 Γ 2 ≪ 1, there is a well defined intermediate time scale τ 2 ≪ t ≪ τ 1 in which the dynamics is apparently stable, as the contributions of the modes with j ≥ 2 have already decayed out while the final decay is not yet appreciable as Γ 1 t ≪ 1. Thus ⟨â(t)⟩ ≈ ⟨â⟩ 1 p 1 (0) + ⟨â⟩ 2 p 2 (0) on this time scale. This corresponds to the plateaus displayed for intermediate times in Fig. 3 (a), in between the corresponding markers in the time axis (indicating τ 1,2 ), and which we now show to be intimately related to entrainment.
The tight connection between metastability and entrainment can be established going back to the laboratory frame 4 where for this intermediate timescale we find the approximate solution According to this, the QvdP oscillator displays an apparently stable perfectly entrained subharmonic response. Provided that ω s ≫ γ 1 , 5 the system will oscillate coherently for many cycles at exactly half of the frequency of the forcing until the effects of the incoherent process described by Eqs. (9) and (12) start to be evident at times of the order of τ 1 . This is exemplified in panel (b), in which we compare the approximate non-decaying subharmonic response of Eq. (13) (black-broken line) with the exact result (color solid line) and for intermediate times inside the plateau, finding a very good agreement. While in a long transient entrainment is then achieved in presence of metastability, Eq. (12) also illustrates precisely how the incoherent process between the two metastable phases eventually hinders the 1:2 entrained response (or period-doubled dynamics) in the long-time limit. In particular, as the stationary state is a symmetric mixture of both phases, and since they differ by a π-phase, the expected value of the amplitude eventually decays out on the long timescale given by τ 1 following the dynamics given by Eqs. (9). Again, this is to be contrasted with the classical case, in which the system settles in one of the bistable locked-phases oscillating subharmonically forever, due to the absence of fluctuations connecting both phases.

B. Two-time correlations and observed frequency
We now consider the dynamics of the amplitude two-time correlation in the stationary state, as from this twotime correlation a well-known indicator of frequency entrainment can be computed. This is the observed frequency [36,60,65,66], ω obs , defined as the maximum of the emission or power spectrum 6 : with where the subscript 'ss' denote that this correlation is calculated in the stationary state in which two-time correlations only depend on the difference between the time arguments, i.e. ⟨â † (−τ )â(0)⟩ ss = ⟨â † (0)â(τ )⟩ ss , while we have also used the fact that ⟨â † (0)â(τ )⟩ ss = [⟨â † (τ )â(0)⟩ ss ] * (see A). Perfect (subharmonic) entrainment arises when the system oscillates at (half of) the driving frequency. As Eq. (15) is in the rotating frame, the observed frequency will be zero when the system is well entrained, while it is known to be close to the intrinsic detuning, ∆, when there is no synchronization [60,65]. In fact, in the metastable regime the emission spectrum is dominated by the contribution of the longest-lived eigenmode. Intuition about this can be gained from the eigendecomposition of the relevant two-time correlation, Eq. (A6) in A. Since in the metastable regime we have that Γ 1 Γ j≥2 ≪ 1, the resonance associated with this mode should stand out in the spectrum. Indeed, we show in Fig. 3 (c) how the contribution of this mode fits accurately the main peak of the spectrum in this regime. Furthermore this panel also illustrates how Γ 1 diminishes when decreasing γ 2 γ 1 : in the three cases amplification/damping rates and squeezing are varied but maintaining Γ 1 Γ 2 ≈ 0.05, and we can see how diminishing γ 2 γ 1 the resonances become significantly sharper, a signature of the vanishing of Γ 1 in the classical limit [8]. In fact, the decay rate of the dominant mode takes the values Γ 1 γ 1 ≈ (0.017, 0.087, 0.152) for the cases γ 2 γ 1 = (0.1, 1, 3), respectively. We further notice that in Fig. 3 (c) we are varying the squeezing strength and the non-linear damping at the same time, both parameters affecting the width of the peaks. In particular, as we have seen, the width of the main peak is accurately captured by Γ 1 , which diminishes with increasing squeezing strength above the EP, and the rest of parameters fixed [e.g. see Fig.  1 (d)]; as for Γ 1 , it increases with the non-linear damping strength and the rest of parameters fixed (not shown here).
While we recall that the large-τ dynamics of two-time correlations can be written in terms of Eq. (9) [44,56], we find that it is more illustrative to write it in terms of the (fully equivalent) Eq. (4). In particular from Eqs. (A6) and (4) it follows that for Γ 1 Γ 2 ≪ 1: We highlight the factor Tr[σ 1âρss ] in this equation stemming from the quantumness of the model. Indeed a two-time correlation (and thus the emission spectrum) follows from an initial perturbation of the stationary state (hereâρ ss ) [9], and can be interpreted as the unavoidable disturbance of a measurement process. Generally, this factor makes this two-time correlation different from the corresponding one of a classical two-state stochastic process [22], which is found to be C(τ ) = ⟨â⟩ 1 2 e −Γ1τ (see D). Since generally Tr[σ 1âρss ] ≠ ⟨â⟩ 1 , we find that ⟨â † (τ )â(0)⟩ ss ≠ C(τ ). This manifests in the fact that the Fourier transform of C(τ ) is centered at the origin for η ≥ η EP while this is not the case for the quantum case, as we shall see in Fig. 3 (d) and (e). Thus, in contrast to the case of the long-time amplitude dynamics Eq. (12), multi-time correlations do not follow straightforwardly from the corresponding classical law of a two-state stochastic process. In these two-time correlations, the quantum nature of the fluctuations and of the degrees of freedom becomes manifest.
We now proceed to analyze the behavior of the observed frequency. In Figs. 3 (d) and (e) we exemplify for different parameter values how as the squeezing strength η increases the system becomes entrained, i.e. ω obs ∆ goes from one to zero. Notice how the transition to entrainment is sharper the smaller is γ 2 γ 1 (large number of excitations). For η > η EP we have plotted in black-broken lines the observed frequency calculated from the effective long-time dynamics, i.e. Eq. (16). Here, we can observe that the asymptotic decay of ω obs towards zero is very well captured by this approximate calculation. This can be further appreciated in the insets of these two panels, in which the well-entrained region, that of ω obs close to zero, is shown in more detail. This constitutes a further confirmation of what we have been discussing previously: the fact that the entrained dynamics is characterized by the metastable response of the state of the system. Thus the signatures of entrainment in the power spectrum are captured by the effective long-time dynamics between the metastable phases, that follows from the opening of a spectral gap.

VI. DISCUSSION AND CONCLUSIONS
In this work we have established the connection between the quantum entrainment in the squeezed QvdP oscillator and quantum metastability. We have reported that squeezing enables the opening of a spectral gap in the Liouvillian, which leads to a huge separation of timescales in the dynamics: after a short transient of time the system settles into the so-called metastable manifold in a mixture of two metastable states that depends on the initial condition. It then follows an incoherent process between the two parity-broken metastable preferred phases of the entrained oscillator. Indeed, the oscillation frequency settles to the value of half of the forcing. Still, quantum entrainment is ultimately limited by this incoherent process, as on a long time scale the temporal coherence of the subharmonic response eventually decays out. In this sense, future investigations could address the question of how the reported metastable entrained dynamics manifests at the quantum trajectory level, analyzing several possible stochastic unravellings of our master equation [69], while considering the dynamics of different observables and comparing with the fixed point attractor scenario of the driven QvdP oscillator, or that of coupled QvdP oscillators [17].
The distinctive features of the analyzed quantum entrained dynamics stem from the distinct classical attractors introduced by the squeezed forcing. Hence, classical phase bistability becomes phase metastability in the presence of quantum fluctuations, as reflected by the dominant fluctuation mode in the system. This more complex scenario is to be contrasted with quantum entrainment in the driven QvdP oscillator (in the absence of squeezing) [41,65,68]. There the dominant fluctuation modes describe the dynamics around the unique fixed point attractor of the entrained regime [68], which can indeed be analyzed by linearizing the master equation around this stable point [68]. This qualitative difference is behind the reported specific features of the power or incoherent spectrum (i.e. the part of the emission spectrum corresponding to the fluctuations) [60,68]. In fact, in Ref. [68] this fluctuations dynamics is analyzed in detail as well as its signatures in the power spectrum. Fluctuations are shown to display an overdamped regime and an underdamped one, which translate in the power spectrum displaying either a broad peak centered at the driving frequency or displaying sidebands around this frequency, respectively [68]. Moreover, in the overdamped regime the width of the Lorentzian-like peak is reported to increase with the driving strength [60], which means that the fluctuations or perturbations around the fixed point attractor decay faster as the system becomes better entrained. This can be interpreted as this attractor becoming more stable and less susceptible to the effects of fluctuations as the forcing strength is increased, i.e. to the coherent part of the dynamics becoming enhanced as the system becomes entrained, while the incoherent spectrum flattens. On the other hand, in the case of the squeezed QvdP oscillator, the width of the main Lorentzian peak decreases significantly as the squeezing strength is increased with the other parameters fixed as reported here [see the behavior of Γ 1 in Fig. 1 (d)] or in [60]. As our analysis reveals, this contrasting behavior of the two cases is due to the completely different physical origin of this peak, rooted in the fluctuation dynamics between the two fixed point attractors of the entrained regime. Thus, the decreasing of this linewidth is to be interpreted as the jumps between the two preferred phases becoming suppressed (or less frequent), which also leads to an enhanced temporal coherence for the entrained subharmonic response. Therefore, we conclude that a careful assessment of the fluctuation dynamics is crucial to understand the synchronized response of quantum systems and its key signatures in the power spectrum.
Precisely, the characteristic fluctuation dynamics reported here transcends the particular context of the squeezed QvdP oscillator. As discussed, metastability appears in many different driven-dissipative quantum systems [44]. A particularly relevant example for us is the case of DPTs with spontaneous parity-symmetry breaking [51]. This is indeed what happens in this system in the infinite-excitation limit where the classical attractors emerge, and the paritybroken metastable preferred phases acquire a divergent lifetime [8]. Other examples of metastable dynamics associated with parity-breaking transitions have been reported in Refs. [3,50,51]. Moreover, such long-time scales associated with fluctuations between multiple possible phases have also been reported for quantum systems of parametric oscillators [16], period-tripling oscillators [26], and optomechanical oscillators [67]. Thus, slow-relaxation time scales seem a characteristic feature of multistable dynamical systems in the quantum regime, and a seemingly metastable manifold and dynamics could be expected for these examples as well as in further synchronization scenarios with multiple preferred phases. In this sense, it would also be interesting to investigate whether the metastable entrained dynamics reported here can be seen as a form of quantum activation process [15,48] emerging in a far-from-equilibrium scenario, in which non-linear dissipation plays a fundamental role in shaping the properties of the metastable states.
Finally, the presence of a long-lived collective excitation has been found to be behind the emergence of spontaneous synchronization in a variety of open quantum systems. This is the case of transient synchronization recently reviewed in Ref. [24]. In fact, such an excitation can display a lifetime much longer than the rest, with the extreme possibility of attaining no-decay. The latter corresponds to stationary synchronization, being associated with highly symmetric situations in the presence of noiseless subsystems and decoherence free subspaces [4,6,25,45,46], as well as strong dynamical symmetries [63]. Transient synchronization and the reported metastable entrained dynamics have in common that the emergence of synchronization is associated with a significant timescale separation in the dynamics, i.e. the opening of a spectral gap, that leads to a long-lived synchronized response practically independent of the initial conditions [24]. In a bigger picture, these examples also illustrate the fact that, while synchronization in quantum systems is possible, the temporal coherence of such a dynamical response is often limited by the presence of quantum fluctuations inherent in open quantum systems, generally leading to a finite time-window for the observation of the synchronized dynamics, in stark contrast to the case of classical noiseless dynamical systems [2,55]. and analogously for the other correlation. In terms of the Liouvillian eigenmodes we can write: Tr[σ † jâρ ss ]Tr[â †ρ j ]e λj τ , in which we find again that the term j = 0 does not contribute as ⟨â⟩ ss = 0 due to parity symmetry. Otherwise, in the absence of parity symmetry there can be a non-zero contribution from j = 0, i.e. an additional term that reads ⟨â † ⟩⟨â⟩ ss . Since λ 0 = 0, this term is time-independent and yields a Dirac delta when Fourier transformed. This contribution is known in some contexts as the coherent part of the emission spectrum [9]. defined as D(μ 1 ,μ ′ 1 ) = Tr[ (μ 1 −μ ′ 1 ) † (μ 1 −μ ′ 1 )] 2, see Fig. 4. Notice that in this figure we have assigned the value one to the region below the EP in which Γ 1 = Γ 2 , which we have delimited by a green-dotted line in order to highlight that it is not part of the analysis. Moreover, we have used a red-dashed line to indicate the contour Γ 1 Γ 2 = 0.1. Comparing Fig. 4 (a), (b) with Fig. 1 (a), (b), we can see how the trace distance becomes vanishingly small as Γ 1 Γ 2 becomes small. Indeed both quantities seem to be very well correlated, as can be checked from these figures. This means that in the regions in which there is a huge separation of timescales, and hence it is meaningful to define the effective dynamics (9), we have thatμ 1 ≈μ ′ 1 to a very good approximation. The same result is found forμ 2 and µ ′ 2 (not shown). Then, as commented in the main text, in the metastable regime we can use the more convenient expressions for the EMSs and also for the projectors:μ ′ 1,2 andP ′ 1(2) ≈ [I + (−)σ 1 ] 2.