Metastable quantum entrainment

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 [1–8]. 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 [9]. 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 [9], the coexistence of different phases in an open Ising model [10] or in a nonlinear oscillator [11]. Moreover, it has also been reported in a dissipative discrete time crystal with Rydberg atoms [12], in the long-time response of a driven-dissipative Rabi model [13], or in a chain of interacting fermions subject to a localized particle loss [14]. A general approach toward metastability in open quantum systems has been recently introduced in reference [9], 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 [15, 16], or as a precursor of eventual spontaneous symmetry breaking in the thermodynamic or infinite-excitation limits of these systems [17]. Indeed, sudden changes can occur in the stationary state of these systems in these limits, what are known as dissipative phase transitions (DPTs) [17, 18]. As in the case of quantum phase transitions, DPTs are associated with a closure of the spectral gap of the generator of the dynamics [17, 18]. 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 [9]. While for first order DPTs this usually occurs in a narrow region between the different regimes [16, 19], as experimentally observed in references [2, 3], for symmetry-breaking DPTs this occurs in a whole dynamical regime characterized by the emergence of metastable symmetry-broken states [17, 20, 21].

A paradigmatic phenomenon of classical driven-dissipative systems is synchronization, which emerges in widely different contexts and forms [22, 23]. Entrainment is a type of synchronization phenomenon in which a self-oscillating system [24] 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 [25]. As for its classical counterpart, diverse approaches and systems have been considered, since the seminal works of references [26, 27]. Different signatures and indicators of quantum synchronization have been explored in the literature [25, 28–35], while its emergence has been reported in disparate scenarios as systems of spins [31, 36, 37], networks of quantum harmonic oscillators [38, 39], optomechanical systems [40, 41], atomic lattices [42, 43] and clouds [44, 45], and micromasers [46]. Moreover, novel forms of synchronization have been identified [35, 43, 47, 48], as the case of transient synchronization [49], 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 [50] and long-lived correlations [28, 36, 38, 39, 51, 52].

A rich playground where synchronization can be extended from the classical to the quantum regime is offered by the paradigmatic van der Pol (vdP) oscillator [22, 23], introduced in references [29, 30]. For this model, both entrainment by a harmonic drive and spontaneous synchronization between coupled oscillators have been found in the quantum regime [29, 30, 47, 53]. 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 [54] 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 [29, 30]. Moreover, the frequency of the maximum of the power spectrum is reported to shift toward that of the driving as the system becomes entrained [30]. 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 [30]. As an alternative implementation strategy, a squeezed forcing has been reported to enhance entrainment in the quantum regime [55]. However, very deep into it, where the system behaves effectively as a few-level system, the enhancing effect of squeezing is limited [56].

The introduction of the squeezed forcing comes at the expense of changing qualitatively the dynamical scenario of entrainment [55, 57, 58]. 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) [22, 59], here it corresponds to bistability (in the rotating frame) [55, 58], 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 [59]. 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 [59]. 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.

In the laboratory frame, the squeezed vdP oscillator is described by the time-dependent master equation for the state of the system ${\hat{\rho }}_{L}$ (ℏ = 1) [55]:

$$\begin{equation}\begin{aligned}{\partial }_{t}{\hat{\rho }}_{L}& =-i[{\hat{H}}_{L}(t),{\hat{\rho }}_{L}]+\frac{{\gamma }_{1}}{2}\mathcal{D}[{\hat{a}}^{{\dagger}}]{\hat{\rho }}_{L}+\frac{{\gamma }_{2}}{2}\mathcal{D}[{\hat{a}}^{2}]{\hat{\rho }}_{L},\\ {\hat{H}}_{L}(t)& ={\omega }_{0}{\hat{a}}^{{\dagger}}\hat{a}+i\eta ({\hat{a}}^{2}{\mathrm{e}}^{\mathrm{i}2{\omega }_{s}t}-{\hat{a}}^{{\dagger}2}{\mathrm{e}}^{-\mathrm{i}2{\omega }_{s}t}),\end{aligned}\end{equation} \tag{ 1 }$$

where we have defined the dissipator $\mathcal{D}[\hat{L}]\hat{\rho }=2\hat{L}\hat{\rho }{\hat{L}}^{{\dagger}}-{\hat{L}}^{{\dagger}}\hat{L}\hat{\rho }-\hat{\rho }{\hat{L}}^{{\dagger}}\hat{L}$, $\hat{L}$ is the corresponding jump operator, $\hat{a}$ is the annihilation operator of the bosonic mode, γ₁ is the amplification rate, γ₂ the nonlinear damping rate, η is the squeezing strength of frequency 2ω_s , and ω₀ 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, ${\hat{U}}_{t}=\mathrm{exp}(-i{\omega }_{s}{\hat{a}}^{{\dagger}}\hat{a}t)$. This leads to a master equation with the same dissipative part but with a time-independent Hamiltonian ¹ : $\hat{H}={\Delta}{\hat{a}}^{{\dagger}}\hat{a}+i\eta ({\hat{a}}^{2}-{\hat{a}}^{{\dagger}2})$ with Δ = ω₀ − ω_s . Notice that we denote the state of the system in this frame as $\hat{\rho }$, 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 [29, 55] and optomechanical oscillators [30, 53, 55]. In both cases, the QvdP oscillator corresponds to the properly engineered effective dynamics of a mechanical degree of freedom. Then, amplification and nonlinear 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 equation (1) have been presented in [55] too, considering both optomechanical systems (e.g. modulating electrically the spring constant of the mechanical mode [60]) and trapped ions (e.g. by a combination of two Raman beams detuned by 2ω_s [61]).

In this rotating frame, the master equation can be written in terms of the time-independent Liouvillian superoperator, defined by ${\partial }_{t}\hat{\rho }=\mathcal{L}\hat{\rho }$. 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 $\mathcal{L}{\hat{\rho }}_{j}={\lambda }_{j}{\hat{\rho }}_{j}$ and ${\mathcal{L}}^{{\dagger}}{\hat{\sigma }}_{j}={\lambda }_{j}^{{\ast}}{\hat{\sigma }}_{j}$, respectively [17, 62] (more details in appendix A). Generally, we can use the Liouvillian eigenmodes to decompose the state of the system at any time

$$\begin{equation}\hat{\rho }(t)={\hat{\rho }}_{ss}+\sum\limits _{j\geqslant 1}\mathrm{Tr}[{\hat{\sigma }}_{j}^{{\dagger}}\hat{\rho }(0)]{\hat{\rho }}_{j}{e}^{{\lambda }_{j}t}.\end{equation} \tag{ 2 }$$

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 ${\hat{\rho }}_{ss}={\hat{\rho }}_{0}/\mathrm{Tr}[{\hat{\rho }}_{0}]$, being the corresponding zero eigenvalue λ₀ = 0. Notice that the real part of these eigenvalues are non-positive, and we order them such that Re[λ₀] ⩾ Re[λ₁] ⩾ Re[λ₂] ⩾ .... In the following, we will refer to the decay rates and frequencies of the eigenmodes

$$\begin{equation}{{\Gamma}}_{j}=\vert \text{Re}[{\lambda }_{j}]\vert \quad \;\text{and}\;\quad {\nu }_{j}=\text{Im}[{\lambda }_{j}],\end{equation} \tag{ 3 }$$

respectively.

While the formalism is fully quantum, we can identify the imbalance between nonlinear dissipation and linear amplification γ₂/γ₁ as the physical quantity controlling the excitation strength (i.e. boson number) in the family of QvdP systems [29, 30], and hence the classical versus quantum system operation. In the limit γ₂/γ₁ ≫ 1 the QvdP oscillator is well approximated by a two-level system [30, 47, 56], while in the limit γ₂/γ₁ ≪ 1 the typical number of excitations becomes macroscopic. Indeed, in reference [63] we analyze in detail how the classical attractors emerge in this model as γ₂/γ₁ → 0. Here instead, we will focus on intermediate values of γ₂/γ₁ for which entrainment in the quantum regime has been reported [55].

A brief overview of the classical dynamics is presented in appendix 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 π.

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.

3.1. 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 Γ₁ 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 figure 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.

Zoom In Zoom Out Reset image size

The spectral gap is characterized through the ratio Γ₁/Γ₂ and can vary several orders of magnitude depending on the squeezing strength and the non-linear damping, figures 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 γ₂/γ₁ and small squeezing η/γ₁ (figure 1(a)) or to larger non-linear damping while increasing the squeezing (figure 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 [64, 65], at which λ_1,2 become the same (figures 1(c) and (d)) while the corresponding eigenmatrices coalesce (not shown). The green-dotted line in figures 1(a) and (b) indicates the squeezing strength at which the EP occurs (η_EP) as γ₂/γ₁ is varied. For η < η_EP, λ_1,2 are complex conjugates, hence explaining the fact that Γ₁/Γ₂ = 1 for this region. While for η > η_EP, λ_1,2 become real valued and Γ₁ decreases when increasing the squeezing strength, while Γ₂ increases, as exemplified in figure 1(d).

Although the behavior of Γ₁/Γ₂ is qualitatively the same for all the considered range of γ₂/γ₁, we find that there are important quantitative differences. In particular, the smaller is γ₂/γ₁ the wider is the gap for a given squeezing strength. This is intimately related to the behavior of Γ₁ as the classical limit is approached (i.e. γ₂/γ₁ → 0). Indeed, in this limit it is found that Γ₁ → 0 for η > η_c while Γ₂ saturates to a non-zero value [63]. This means that the steady state becomes degenerate and the system undergoes a DPT [17, 18]. This is intimately related to the emergence of the classical bistable attractors as we explain in detail in [63]. Notice that in this limit it is also found that η_EP → η_c .

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 Γ₁ and Γ₂ differing by an order of magnitude or more, as in the (right) regions delimited by a red-dashed line (where Γ₁/Γ₂ ⩽ 0.1) in figures 1(a) and (b).

3.2. Effective long-time dynamics

The disparity between the weakest damping rate Γ₁ and the others allows to capture the system dynamics after a transient considering only the disparate timescales ${\tau }_{1,2}^{-1}={{\Gamma}}_{1,2}$. This large separation of time scales makes an initially excited state to rapidly decay (∼τ₂) to the manifold spanned by ${\hat{\rho }}_{ss}$ and the longest-lived eigenmode ${\hat{\rho }}_{1}$, where it displays a slow relaxation to ${\hat{\rho }}_{ss}$ at times of the order of τ₁ [9, 10]. In the presence of a significant spectral gap, as reported in the previous section, there is a well-defined intermediate time scale, τ₂ ≪ t ≪ τ₁, in which the contribution of the higher modes is negligible while the decay of ${\hat{\rho }}_{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 [9], as we develop further in this section.

After times of the order of τ₂, the dynamics is well-approximated just considering the contributions of the longest-lived eigenmode and of the stationary state. As we have shown in the previous section, the parameter regime in which Γ₁/Γ₂ ≪ 1 corresponds to η > η_EP, and thus λ_1,2 are real valued while the corresponding right and left eigenmatrices are Hermitian [17, 62]. Then the long-time dynamics can be approximated as

$$\begin{equation}\hat{\rho }(t)\approx {\hat{\rho }}_{ss}+\mathrm{Tr}[{\hat{\sigma }}_{1}\hat{\rho }(0)]{\hat{\rho }}_{1}{e}^{-{{\Gamma}}_{1}t},\end{equation} \tag{ 4 }$$

which follows from equation (2) neglecting the contributions for the modes with j ⩾ 2. By using the formalism introduced in references [9, 10], the approximate long-time dynamics of equation (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 ${\hat{\rho }}_{1}$, the latter not being a physical state since it is traceless [17]. In the following, we discuss the application of such formalism [9, 10] to our system (further details are in appendix C).

From the long-time approximation of equation (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:

$$\begin{equation}\mathcal{P}\hat{\rho }={\hat{\rho }}_{ss}+\mathrm{Tr}[{\hat{\sigma }}_{1}\hat{\rho }]{\hat{\rho }}_{1}.\end{equation} \tag{ 5 }$$

We denote this manifold as the metastable manifold, as it captures the state of the system for the intermediate time scale τ₂ ≪ t ≪ τ₁. In references [9, 10], it is shown that the metastable manifold can be parametrized in terms of the extreme metastable states (EMSs), defined as:

$$\begin{equation}{\hat{\mu }}_{1}={\hat{\rho }}_{ss}+{c}_{\mathrm{max}}{\hat{\rho }}_{1},\hspace{25.0pt}{\hat{\mu }}_{2}={\hat{\rho }}_{ss}+{c}_{\mathrm{min}}{\hat{\rho }}_{1},\end{equation} \tag{ 6 }$$

where c_max and c_min are the maximum and minimum eigenvalues of ${\hat{\sigma }}_{1}$. In our system, numerical analysis reveals that c_min = −c_max for the considered regime, i.e. when Γ₁/Γ₂ < 1. Moreover, we find that for a significant spectral gap (Γ₁/Γ₂ ≪ 1) these coefficients can be well approximated by c_max = −c_min ≈ 1, which yields the simpler expressions for the EMSs: ${\hat{\mu }}_{1(2)}\approx {\hat{\rho }}_{ss}+(-){\hat{\rho }}_{1}$ (see appendix 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:

$$\begin{equation}\mathcal{P}\hat{\rho }(t)={p}_{1}(t){\hat{\mu }}_{1}+{p}_{2}(t){\hat{\mu }}_{2},\end{equation} \tag{ 7 }$$

where the real coefficients p_1,2(t) are defined in appendix C, while they are constraint to satisfy [9, 10]:

$$\begin{equation}{p}_{1,2}(t)\geqslant 0,\quad {p}_{1}(t)+{p}_{2}(t)=1.\end{equation} \tag{ 8 }$$

By means of the EMSs the long-time dynamics can be parametrized in terms of two physical states, ² whose properties can be readily characterized. Crucially, the EMSs constitute the only basis for which p_1,2(t) satisfy equation (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 [9, 10]:

$$\begin{equation}{\dot{p}}_{1}(t)=-{\dot{p}}_{2}(t)=-\frac{{{\Gamma}}_{1}}{2}\left[{p}_{1}(t)-{p}_{2}(t)\right],\end{equation} \tag{ 9 }$$

whose solution makes equation (7) equivalent to equation (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 Γ₁/2 [66]. This characterizes the slow dynamics to which an initially excited state rapidly relaxes. From the solution of this dynamics (appendix C) we find that p₁(t → ∞) = p₂(t → ∞) = 1/2, which could be also deduced from the fact that ${\hat{\rho }}_{ss}=({\hat{\mu }}_{1}+{\hat{\mu }}_{2})/2$. This illustrates the notion of metastability in our system: any initial unbalanced mixture of the EMSs decays on a very long timescale to the balanced one (${\hat{\rho }}_{ss}$), remaining apparently stable for the well-defined intermediate timescale τ₂ ≪ t ≪ τ₁. In the following section we explore in detail the physical meaning of this statement by characterizing the properties of ${\hat{\mu }}_{1,2}$, through a visual representation based on the Wigner function, and of this effective incoherent dynamics between them.

As we show here, the characterization of the EMSs ${\hat{\mu }}_{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 figure 2 we plot the Wigner distribution for ${\hat{\rho }}_{ss}$, and ${\hat{\mu }}_{1,2}$. As illustrated in figures 2(a) and (d) the Wigner distribution of the stationary state displays a manifest bimodal character. Here we consider two disparate parameter values γ₂/γ₁ = (0.1, 3) and η/γ₁ = (0.2, 2), both satisfying Γ₁/Γ₂ ≪ 1 as can be checked from figure 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 ${\hat{\mu }}_{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.

Zoom In Zoom Out Reset image size

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 equation (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 τ₁ 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 [29, 30, 35, 55], or between coupled oscillators [40, 47, 53]. 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 [29, 30]. The latter corresponds to phase-locking in the quantum regime, in which quantum fluctuations play a non-negligible role [59]. Moving to the case of the squeezed QvdP oscillator, the emergence of phase preference is accompanied by the bimodality of the Wigner distribution [55] as displayed in figure 2. Therefore, each of the EMSs ${\hat{\mu }}_{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 ${\hat{\rho }}_{ss}$ and ${\hat{\rho }}_{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 $\hat{a}\to -\hat{a}$, ${\hat{a}}^{{\dagger}}\to -{\hat{a}}^{{\dagger}}$, i.e. it displays parity-symmetry [17, 62]. This means that we can define the unitary transformation ${\mathcal{Z}}_{2}\hat{\rho }={\text{e}}^{\text{i}\pi {\hat{a}}^{{\dagger}}\hat{a}}\hat{\rho }{\text{e}}^{-\text{i}\pi {\hat{a}}^{{\dagger}}\hat{a}}$, whose action commutes with that of the Liouvillian $[{\mathcal{Z}}_{2},\mathcal{L}]\hat{\rho }=0$. Hence, the eigenmodes are either parity symmetric or parity antisymmetric, which means ${\mathcal{Z}}_{2}{\hat{\rho }}_{j}={z}_{j}{\hat{\rho }}_{j}$ with z_j = ±1, respectively [17, 62]. An important observation is that the stationary state and the longest-lived eigenmode belong to different symmetry sectors:

$$\begin{equation}{\mathcal{Z}}_{2}{\hat{\rho }}_{ss}={\hat{\rho }}_{ss},\hspace{25.0pt}{\mathcal{Z}}_{2}{\hat{\rho }}_{1}=-{\hat{\rho }}_{1}.\end{equation} \tag{ 10 }$$

While the stationary state must be parity symmetric [17, 62], the symmetry of ${\hat{\rho }}_{1}$ is assessed numerically by realizing that its contribution for parity symmetric observables vanishes identically, i.e. $\mathrm{Tr}[{\hat{\rho }}_{1}{({\hat{a}}^{{\dagger}})}^{m}{\hat{a}}^{n}]=0$ for m + n = even, while its contribution for parity antisymmetric ones is generally non-zero, i.e. $\mathrm{Tr}[{\hat{\rho }}_{1}{({\hat{a}}^{{\dagger}})}^{m}{\hat{a}}^{n}]\ne 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:

$$\begin{equation}{\mathcal{Z}}_{2}{\hat{\mu }}_{1(2)}={\hat{\mu }}_{2(1)}.\end{equation} \tag{ 11 }$$

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 $\mathcal{P}\hat{\rho }(t)$ (see equation (7)) both ${\hat{\mu }}_{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 ³ in which since it is an antisymmetric observable we have ${\langle \hat{a}\rangle }_{1}=-{\langle \hat{a}\rangle }_{2}$. This is precisely the same phase relation between the two possible classical amplitudes in the synchronized regime, as briefly commented in section 2. 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 [63]. These observations together with the Wigner representations of figure 2, suggest an intimate connection between the emergence of preferred phases in the quantum regime, and thus synchronization, and metastability.

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.

5.1. Amplitude dynamics

After a short transient of the order of τ₂, the amplitude dynamics for an arbitrary initial state is restricted to the metastable manifold and is well approximated by

$$\begin{equation}\langle \hat{a}(t)\rangle \approx {\langle \hat{a}\rangle }_{1}{p}_{1}(t)+{\langle \hat{a}\rangle }_{2}{p}_{2}(t),\end{equation} \tag{ 12 }$$

where p_1,2(t) are the solutions of equation (9). Notice that the expectation value of the amplitude is defined as $\langle \hat{a}(t)\rangle =\mathrm{Tr}[\hat{a}\hat{\rho }(t)]$, which can be written in terms of the Liouvillian eigenmodes as detailed in appendix A. Then, within the approximation of the long-time transient given in equation (4), and using the EMSs decomposition of the metastable manifold, we can obtain the expression given in equation (12). The accuracy of equation (12) is numerically confirmed in figure 3(a) for three different values of γ₂/γ₁. In this panel we plot in logarithmic scale the imaginary part of $\langle \hat{a}(t)\rangle$ (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 equation (12). Notice that equation (12) is the same we would obtain from the two-state stochastic process using the rules of classical stochastics [66], 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 appendix D).

Zoom In Zoom Out Reset image size

The behavior observed in figure 3(a) can be fully understood recalling the Liouvillian analysis of section 3, where we found that whenever Γ₁/Γ₂ ≪ 1, there is a well defined intermediate time scale τ₂ ≪ t ≪ τ₁ 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 Γ₁ t ≪ 1. Thus $\langle \hat{a}(t)\rangle \approx {\langle \hat{a}\rangle }_{1}{p}_{1}(0)+{\langle \hat{a}\rangle }_{2}{p}_{2}(0)$ on this time scale. This corresponds to the plateaus displayed for intermediate times in figure 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 ⁴ where for this intermediate timescale we find the approximate solution

$$\begin{equation}\langle \hat{a}{(t)\rangle }_{L}\approx \langle {\hat{a}\rangle }_{1}{\text{e}}^{-\text{i}{\omega }_{s}t}({p}_{1}(0)-{p}_{2}(0)),\end{equation} \tag{ 13 }$$

According to this, the QvdP oscillator displays an apparently stable perfectly entrained subharmonic response. Provided that ω_s ≫ γ₁, ⁵ 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 equations (9) and (12) start to be evident at times of the order of τ₁. This is exemplified in panel (b), in which we compare the approximate non-decaying subharmonic response of equation (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, equation (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 τ₁ following the dynamics given by equation (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.

5.2. Two-time correlations and observed frequency

We now consider the dynamics of the amplitude two-time correlation in the stationary state, as from this two-time correlation a well-known indicator of frequency entrainment can be computed. This is the observed frequency [30, 53, 55, 57], ω_obs, defined as the maximum of the emission or power spectrum: ⁶

$$\begin{equation}{\omega }_{\text{obs}}=\text{argmax}[S(\omega )]\end{equation} \tag{ 14 }$$

with

$$\begin{equation}S(\omega )={\int }_{-\infty }^{\infty }\mathrm{d}\tau \enspace {\text{e}}^{-\text{i}\omega \tau }{\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}(0)\rangle }_{ss}=2\enspace \text{Re}\left[{\int }_{0}^{\infty }\mathrm{d}\tau \enspace {\text{e}}^{-\text{i}\omega \tau }{\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}(0)\rangle }_{ss}\right],\end{equation} \tag{ 15 }$$

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. ${\langle {\hat{a}}^{{\dagger}}(-\tau )\hat{a}(0)\rangle }_{ss}={\langle {\hat{a}}^{{\dagger}}(0)\hat{a}(\tau )\rangle }_{ss}$, while we have also used the fact that ${\langle {\hat{a}}^{{\dagger}}(0)\hat{a}(\tau )\rangle }_{ss}={[{\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}(0)\rangle }_{ss}]}^{{\ast}}$ (see appendix A).

Perfect (subharmonic) entrainment arises when the system oscillates at (half of) the driving frequency. As equation (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 [30, 55]. 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, equation (A.6) in appendix A. Since in the metastable regime we have that Γ₁/Γ_j⩾2 ≪ 1, the resonance associated with this mode should stand out in the spectrum. Indeed, we show in figure 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 Γ₁ diminishes when decreasing γ₂/γ₁: in the three cases amplification/damping rates and squeezing are varied but maintaining Γ₁/Γ₂ ≈ 0.05, and we can see how diminishing γ₂/γ₁ the resonances become significantly sharper, a signature of the vanishing of Γ₁ in the classical limit [63]. In fact, the decay rate of the dominant mode takes the values Γ₁/γ₁ ≈ (0.017, 0.087, 0.152) for the cases γ₂/γ₁ = (0.1, 1, 3), respectively. We further notice that in figure 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 Γ₁, which diminishes with increasing squeezing strength above the EP, and the rest of parameters fixed (e.g. see figure 1(d)); as for Γ₁, it increases with the nonlinear 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 equation (9) [9, 10], we find that it is more illustrative to write it in terms of the (fully equivalent) equation (4). In particular from equations (A.6) and (4) it follows that for Γ₁/Γ₂ ≪ 1:

$$\begin{equation}{\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}(0)\rangle }_{ss}\approx \mathrm{Tr}[{\hat{\sigma }}_{1}\hat{a}{\hat{\rho }}_{ss}]\mathrm{Tr}[{\hat{a}}^{{\dagger}}{\hat{\rho }}_{1}]{e}^{-{{\Gamma}}_{1}\tau }.\end{equation} \tag{ 16 }$$

We highlight the factor $\mathrm{Tr}[{\hat{\sigma }}_{1}\hat{a}{\hat{\rho }}_{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 $\hat{a}{\hat{\rho }}_{ss}$) [54], 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 [66], which is found to be $C(\tau )=\vert {\langle \hat{a}\rangle }_{1}{\vert }^{2}{e}^{-{{\Gamma}}_{1}\tau }$ (see appendix D). Since generally $\mathrm{Tr}[{\hat{\sigma }}_{1}\hat{a}{\hat{\rho }}_{ss}]\ne {\langle \hat{a}\rangle }_{1}$, we find that ${\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}(0)\rangle }_{ss}\ne C(\tau )$. 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 figures 3(d) and (e). Thus, in contrast to the case of the long-time amplitude dynamics equation (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 figures 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 γ₂/γ₁ (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. equation (16). Here, we can observe that the asymptotic decay of ω_obs toward 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.

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 [1], 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 [67].

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) [29, 30, 59]. There the dominant fluctuation modes describe the dynamics around the unique fixed point attractor of the entrained regime [59], which can indeed be analyzed by linearizing the master equation around this stable point [59]. 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) [55, 59]. In fact, in Ref. [59] 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 [59]. Moreover, in the overdamped regime the width of the Lorentzian-like peak is reported to increase with the driving strength [55], 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 Γ₁ in figure 1(d)) or in [55]. 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 [9]. A particularly relevant example for us is the case of DPTs with spontaneous parity-symmetry breaking [17]. This is indeed what happens in this system in the infinite-excitation limit where the classical attractors emerge, and the parity-broken metastable preferred phases acquire a divergent lifetime [63]. Other examples of metastable dynamics associated with parity-breaking transitions have been reported in references [17, 20, 21]. Moreover, such long timescales associated with fluctuations between multiple possible phases have also been reported for quantum systems of parametric oscillators [68], period-tripling oscillators [69], and optomechanical oscillators [41]. Thus, slow-relaxation timescales 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 [70, 71] emerging in a far-from-equilibrium scenario, in which nonlinear 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 reference [49]. 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 [28, 38, 39, 50, 51], as well as strong dynamical symmetries [43]. 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 [49]. 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 [22, 23].

We acknowledge interesting comments from Rosario Fazio and the use of the python package QuTip [72, 73] in order to obtain some of the numerical results. We acknowledge the Spanish State Research Agency, through the QUARESC Project QUARESC Project (PID2019-109094GB-C21/AEI/10.13039/501100011033) and through the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R&D (MDM-2017-0711); We also acknowledge funding by CAIB through the QUAREC project (PRD2018/47). AC is funded by the CAIB PhD program. GLG is funded by the Spanish Ministerio de Educación y Formación Profesional/Ministerio Universidades and co-funded by the University of the Balearic Islands through the Beatriz Galindo program (BG20/00085).

The data that support the findings of this study are available upon reasonable request from the authors.

As commented in the main text, the Liouvillian superoperator can be diagonalized obtaining a set of eigenvalues λ_j , and the corresponding right and left eigenmatrices ${\hat{\rho }}_{j}$ and ${\hat{\sigma }}_{j}^{{\dagger}}$, respectively. These eigenvalues are non-positive and either real or appear in complex conjugate pairs. In the former case, the corresponding eigenmatrices can be defined Hermitian [17]. Notice that except for spectral singularities [65], these eigenmatrices form a biorthogonal basis that can be normalized by means of the Hilbert–Schmidt product: $\mathrm{Tr}[{\hat{\sigma }}_{j}^{{\dagger}}{\hat{\rho }}_{k}]={\delta }_{jk}$. Moreover, in case the real part of the eigenvalues is negative, the corresponding eigenmodes are traceless, and hence they are not physical states [17]. In this case we can choose whether to rescale ${\hat{\sigma }}_{j}^{{\dagger}}$ or ${\hat{\rho }}_{j}$ freely, as long as they satisfy the normalization condition. The stationary state is obtained from the properly normalized right eigenmatrix with zero eigenvalue, i.e. ${\hat{\rho }}_{ss}={\hat{\rho }}_{0}/\mathrm{Tr}[{\hat{\rho }}_{0}]$ while the left one is the identity ${\hat{\sigma }}_{0}=I$.

Since except for spectral singular points these eigenmatrices form a basis of the system Hilbert space, they can be used to decompose the state of the system at any time as done in equation (2). By taking the trace over a given operator we can obtain its dynamics. In particular for the amplitude dynamics we obtain:

$$\begin{equation}\langle \hat{a}(t)\rangle =\sum\limits _{j\geqslant 1}\mathrm{Tr}[{\hat{\sigma }}_{j}^{{\dagger}}\hat{\rho }(0)]\mathrm{Tr}[\hat{a}{\hat{\rho }}_{j}]{e}^{{\lambda }_{j}t}.\end{equation} \tag{ A.1 }$$

Here we have used the fact that the Liouvillian displays a parity symmetry, and thus the eigenmodes are either parity symmetric or parity antisymmetric, and the stationary state is parity symmetric. In fact, it follows that for operators with a given symmetry, only eigenmodes with the same symmetry contribute to their dynamics. For the particular case of the amplitude only parity antisymmetric eigenmodes contribute to the sum, as $\mathrm{Tr}[\hat{a}{\hat{\rho }}_{j}]=0$ is identically zero if ${\mathcal{Z}}_{2}{\hat{\rho }}_{j}={\hat{\rho }}_{j}$. Notice that this is also the case for the amplitude two-time correlations calculated in the stationary state. But first, let us recall how to compute these two-time correlations from the master equation [54]:

$$\begin{equation}\langle {\hat{a}}^{{\dagger}}({t}_{1})\hat{a}({t}_{2})\rangle =\mathrm{Tr}[{\hat{a}}^{{\dagger}}{e}^{\mathcal{L}({t}_{1}-{t}_{2})}\left(\hat{a}\hat{\rho }({t}_{2})\right)],\quad \text{if}\enspace {t}_{1}\geqslant {t}_{2},\end{equation} \tag{ A.2 }$$

and

$$\begin{equation}\langle {\hat{a}}^{{\dagger}}({t}_{1})\hat{a}({t}_{2})\rangle =\mathrm{Tr}[\hat{a}{e}^{\mathcal{L}({t}_{2}-{t}_{1})}\left(\hat{\rho }({t}_{1}){\hat{a}}^{{\dagger}}\right)],\quad \text{if}\enspace {t}_{2}\geqslant {t}_{1},\end{equation} \tag{ A.3 }$$

where the superoperator ${e}^{\mathcal{L}\vert {t}_{1}-{t}_{2}\vert }$ applies to everything into the big rounded parenthesis, and it stands for time-evolving for a period of time |t₁ − t₂| and according to $\mathcal{L}$ the 'initial condition' inside the big rounded parenthesis. Following [54], these correlations can also be described using the variables τ ⩾ 0 and t = min(t₁, t₂), from which we obtain:

$$\begin{equation}\langle {\hat{a}}^{{\dagger}}(t+\tau )\hat{a}(t)\rangle =\mathrm{Tr}[{\hat{a}}^{{\dagger}}{e}^{\mathcal{L}\tau }\left(\hat{a}\hat{\rho }(t)\right)],\quad \text{for}\enspace \tau \geqslant 0,\end{equation} \tag{ A.4 }$$

and an analogous expression for the case of (A.3). The stationary two-time correlations studied in the main text are then defined as

$$\begin{equation}\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}{(0)\rangle }_{ss}=\underset{t\to \infty }{\mathrm{lim}}\enspace \langle {\hat{a}}^{{\dagger}}(t+\tau )\hat{a}(t)\rangle =\mathrm{Tr}[{\hat{a}}^{{\dagger}}{e}^{\mathcal{L}\tau }(\hat{a}{\hat{\rho }}_{ss})],\end{equation} \tag{ A.5 }$$

and analogously for the other correlation. In terms of the Liouvillian eigenmodes we can write:

$$\begin{equation}{\langle {\hat{a}}^{{\dagger}}(\tau )\hat{a}(0)\rangle }_{ss}=\sum\limits _{j\geqslant 1}\enspace \mathrm{Tr}[{\hat{\sigma }}_{j}^{{\dagger}}\hat{a}{\hat{\rho }}_{ss}]\mathrm{Tr}[{\hat{a}}^{{\dagger}}{\hat{\rho }}_{j}]{e}^{{\lambda }_{j}\tau },\end{equation} \tag{ A.6 }$$

in which we find again that the term j = 0 does not contribute as ${\langle \hat{a}\rangle }_{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 $\langle {\hat{a}}^{{\dagger}}\rangle {\langle \hat{a}\rangle }_{ss}$. Since λ₀ = 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 [54].

We briefly overview the dynamical regimes that this model displays in the classical limit. In this limit we define the complex amplitude $\alpha =\langle \hat{a}\rangle$, whose equation of motion is obtained from that of $\langle \hat{a}\rangle$ facotrizing higher-order moments:

$$\begin{equation}{\partial }_{t}\alpha =-i{\Delta}\alpha +\frac{{\gamma }_{1}}{2}\alpha -{\gamma }_{2}\vert \alpha {\vert }^{2}\alpha -2\eta {\alpha }^{{\ast}}.\end{equation} \tag{ B.1 }$$

This equation displays a single bifurcation at the classical critical point η_c = |Δ|/2. For η < η_c , the stable attractor is a limit-cycle of fundamental frequency ${\Omega}={\Delta}\sqrt{1-{(2\eta /{\Delta})}^{2}}$. While for η > η_c the stable attractors are two bistable fixed points that only differ by a complex phase e^iπ ; their amplitudes are given by α_± = ±Re^iϕ with $R=\sqrt{\frac{{\gamma }_{1}}{2{\gamma }_{2}}+\frac{1}{{\gamma }_{2}}\sqrt{4{\eta }^{2}-{{\Delta}}^{2}}}$, while the phase is 2ϕ = π − sin⁻¹[|Δ|/(2η)] [55, 58]. Back into the laboratory frame, the limit-cycle regime corresponds to the lack of entrainment since Ω and 2ω_s are not generally commensurate. On the other hand, the bistable solutions become ${\alpha }_{\pm }(t)=\pm {\mathrm{R}\mathrm{e}}^{\mathrm{i}\phi -\mathrm{i}{\omega }_{s}t}$, and thus this regime corresponds to subharmonic entrainment with two possible locked phases [22, 23].

As we have explained in the main text, the long-time dynamics of the system is restricted to the so-called metastable manifold [9, 10]. Notice that since in the whole region in which Γ₁/Γ₂ ≠ 1, λ₁ is real, we find that ${\hat{\sigma }}_{1}$ and ${\hat{\rho }}_{1}$ are Hermitian. This property enables to parametrize the metastable manifold in terms of the EMSs defined in equation (6), as well as the projectors on these states, defined by:

$$\begin{equation}{\hat{P}}_{1}=\frac{1}{{\Delta}c}({\hat{\sigma }}_{1}-{c}_{\mathrm{min}}I),\hspace{25.0pt}{\hat{P}}_{2}=\frac{1}{{\Delta}c}(-{\hat{\sigma }}_{1}+{c}_{\mathrm{max}}I),\end{equation} \tag{ C.1 }$$

with Δc = c_max − c_min, while c_max and c_min are the maximum and minimum eigenvalues of ${\hat{\sigma }}_{1}$. ⁷ As commented, numerical analysis reveals that c_min = −c_max for the considered regime, i.e. when Γ₁/Γ₂ < 1.

The projection of the initial state onto the metastable manifold can be decomposed as the projection onto the EMSs: $\mathcal{P}\hat{\rho }(0)=\mathrm{Tr}[{\hat{P}}_{1}\hat{\rho }(0)]{\hat{\mu }}_{1}+\mathrm{Tr}[{\hat{P}}_{2}\hat{\rho }(0)]{\hat{\mu }}_{2}$. Importantly, we have that by construction ${\hat{P}}_{1,2}\geqslant 0$, ${\hat{P}}_{1}+{\hat{P}}_{2}=I$ and $\mathrm{Tr}[{\hat{P}}_{i}{\hat{\mu }}_{j}]={\delta }_{ij}$. Moreover, we can define the expectation values ${p}_{1(2)}(t)=\mathrm{Tr}[{\hat{P}}_{1(2)}\hat{\rho }(t)]$. It follows that p_1,2(t) are positive and they sum to one. Hence p_1,2(t) can be interpreted as probabilities as stated in the main text. Notice that this is only the case for the EMSs, as they form the only basis in which ${\hat{P}}_{1,2}$ are positive in all the metastable manifold [9]. The long-time dynamics can be recasted in this basis leading to equation (9), which in virtude of the properties of p_1,2(t) can be interpreted as a stochastic process [9]. Then the solution to this equation is:

$$\begin{equation}{p}_{1(2)}(t)=\frac{{p}_{1(2)}^{0}}{2}(1+{e}^{-{{\Gamma}}_{1}t})+\frac{{p}_{2(1)}^{0}}{2}(1-{e}^{-{{\Gamma}}_{1}t}),\end{equation} \tag{ C.2 }$$

where the initial conditions are given by the projection of the initial state onto the metastable manifold in terms of the EMSs: ${p}_{1(2)}^{0}=\mathrm{Tr}[{\hat{P}}_{1(2)}\hat{\rho }(0)]$.

Finally we compare the exact expressions for the EMSs, i.e. ${\hat{\mu }}_{1}={\hat{\rho }}_{ss}+{c}_{\mathrm{max}}\enspace {\hat{\rho }}_{1}$ with the approximate ones, that we distinguish in this appendix with a prime, i.e. ${\hat{\mu }}_{1}^{\prime }={\hat{\rho }}_{ss}+{\hat{\rho }}_{1}$. We do so by computing the trace distance between both, defined as $D({\hat{\mu }}_{1},{\hat{\mu }}_{1}^{\prime })=\mathrm{Tr}[\sqrt{{({\hat{\mu }}_{1}-{\hat{\mu }}_{1}^{\prime })}^{{\dagger}}({\hat{\mu }}_{1}-{\hat{\mu }}_{1}^{\prime })}]/2$, see figure C1. Notice that in this figure we have assigned the value one to the region below the EP in which Γ₁ = Γ₂, 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 Γ₁/Γ₂ = 0.1. Comparing figures C1(a) and (b) with figures 1(a) and (b), we can see how the trace distance becomes vanishingly small as Γ₁/Γ₂ 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 ${\hat{\mu }}_{1}\approx {\hat{\mu }}_{1}^{\prime }$ to a very good approximation. The same result is found for ${\hat{\mu }}_{2}$ and ${\hat{\mu }}_{2}^{\prime }$ (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: ${\hat{\mu }}_{1,2}^{\prime }$ and ${\hat{P}}_{1(2)}^{\prime }\approx [I+(-){\hat{\sigma }}_{1}]/2$.

Zoom In Zoom Out Reset image size

In this appendix we consider equation (9) as a fully classical two-state stochastic process and we compute the amplitude dynamics and two-time correlations according to the classical rules [66]. In particular, we consider that state 1 is characterized by the complex amplitude ${\alpha }_{1}={\langle \hat{a}\rangle }_{1}$, while state 2 by the complex amplitude ${\alpha }_{2}=-{\langle \hat{a}\rangle }_{1}$. The solution for the classical two-state process is again given by equation (C.2), while the difference with the effective long-time quantum dynamics resides, in principle, in the recipe to compute statistical quantities. Firstly, we have that the classical averaged amplitude is given by [66]:

$$\begin{equation}\bar{\alpha (t)}=\sum\limits _{j=1}^{2}{\alpha }_{j}{p}_{j}(t),\end{equation} \tag{ D.1 }$$

where the overbar indicates a classical average over the stochastic process. As commented in the main text this expression coincides with equation (12). Notice that the stationary value for the amplitude is zero, i.e. $\bar{{\alpha }_{ss}}=0$, as expected. Secondly, the classical recipe for the two-time amplitude correlations in the stationary state is given by [66]:

$$\begin{equation}C(\tau )=\underset{t\to \infty }{\mathrm{lim}}\enspace \bar{{\alpha }^{{\ast}}(t+\tau )\alpha (t)}=\sum\limits _{i=1}^{2}\sum\limits _{j=1}^{2}{\alpha }_{i}^{{\ast}}{\alpha }_{j}{\tilde{p}}_{ij}(\tau ){p}_{j}(t\to \infty ),\end{equation} \tag{ D.2 }$$

where we have defined the tilded probabilities ${\tilde{p}}_{ij}(\tau )$ as the solutions p_i (τ) given in equation (C.2) with the special initial condition ${p}_{j}^{(0)}=1$. From this expression we recover the result quoted in the main text, that is $C(\tau )=\vert {\alpha }_{1}{\vert }^{2}{e}^{-{{\Gamma}}_{1}\tau }$.