Phase synchronization in coupled bistable oscillators

TL;DR

This paper presents a theoretical study of synchronization in quantum bistable oscillators, revealing how coupling induces various phase preferences depending on the dynamical state, with potential implications for quantum synchronization control.

Contribution

It introduces a dissipative quantum oscillator model exhibiting bistability and analyzes how coupling leads to diverse phase synchronization behaviors.

Findings

In the limit-cycle regime, phase peaks at 0 and π.

In the low-occupation regime, phase peaks at π/2 and 3π/2.

Intermediate regimes show π/2-periodic phase distributions.

Abstract

We introduce a simple model system to study synchronization theoretically in quantum oscillators that are not just in limit-cycle states, but rather display a more complex bistable dynamics. Our oscillator model is purely dissipative, with a two-photon gain balanced by single- and three-photon loss processes. When the gain rate is low, loss processes dominate and the oscillator has a very low photon occupation number. In contrast, for large gain rates, the oscillator is driven into a limit-cycle state where photon numbers can become large. The bistability emerges between these limiting cases with a region of coexistence of limit-cycle and low-occupation states. Although an individual oscillator has no preferred phase, when two of them are coupled together a relative phase preference is generated which can indicate synchronization of the dynamics. We find that the form and strength of the relative phase preference varies widely depending on the dynamical states of the oscillators. In the limit-cycle regime, the phase distribution is $\pi$-periodic with peaks at $0$ and $\pi$, whilst in the low-occupation regime $\pi$-periodic phase distributions can be produced with peaks at $\pi/2$ and $3\pi/2$. Tuning the coupled system between these two regimes reveals a region where the relative phase distribution has $\pi/2$-periodicity.

Equations

$$\dot{\rho }=\mathcal{L}\rho ={\kappa }_1\mathcal{D}[\hat{a}]\left(\rho \right)+{\kappa }_2\mathcal{D}[({\hat{a}}^{†}{)}^2]\left(\rho \right)+{\kappa }_3\mathcal{D}[{\hat{a}}^3]\left(\rho \right),$$

$${\dot{\rho }}^{(k)}={\mathcal{M}}^{(k)}{\rho }^{(k)},$$

$${\dot{P}}_n=-G_n{P}_n+A_{n+1}{P}_{n+1}+B_{n-2}{P}_{n-2}+C_{n+3}{P}_{n+3},$$

$$G_n=$$

$$A_{n+1}=$$

$$B_{n-2}=$$

$$C_{n+3}=$$

$$0<\frac{3{\kappa }_3}{{\kappa }_2}\left(\frac{{\kappa }_1}{{\kappa }_2}-4\right)<1.$$

$$M=\frac{\mathrm{Re}({\lambda }_2^{(0)}-{\lambda }_1^{(0)})}{\mathrm{Re}({\lambda }_1^{(0)}-{\lambda }_0^{(0)})},$$

$${\hat{H}}_J=\hslash J\left({\hat{a}}_1^{†}{\hat{a}}_2+{\hat{a}}_1{\hat{a}}_2^{†}\right),$$

$$\dot{\rho }=-\frac{i}{\hslash }\left[{\hat{H}}_J,\rho \right]+\sum _{j=1,2}{\mathcal{L}}_j\rho ,$$

$$P(\varphi )$$

$$P(\varphi )=\frac{1}{2\pi }+J^{2}f({\kappa }_1,{\kappa }_2)[P_1^2-P_0{P}_2]\cos (2\varphi ),$$

$$\frac{d}{dt}⟨\hat{a}⟩=-\frac{{\kappa }_1}{2}⟨\hat{a}⟩+{\kappa }_2\left(2⟨\hat{a}⟩+⟨{\hat{a}}^{†}{\hat{a}}^2⟩\right)-\frac{3{\kappa }_3}{2}⟨({\hat{a}}^{†}{)}^2{\hat{a}}^3⟩.$$

$$\dot{r}+ir\dot{\phi }=-\frac{3{\kappa }_3}{2}{r}^5+{\kappa }_2{r}^3+\left(2{\kappa }_2-\frac{{\kappa }_1}{2}\right)r.$$

$$\dot{r}=\left[-\frac{3{\kappa }_3}{2}{r}^4+{\kappa }_2{r}^2+\left(2{\kappa }_2-\frac{{\kappa }_1}{2}\right)\right]r,$$

$$n_{\pm }=\frac{{\kappa }_2}{3{\kappa }_3}\left[1\pm \sqrt{1+\frac{3{\kappa }_3\left(4{\kappa }_2-{\kappa }_1\right)}{{\kappa }_2^2}}\right].$$

$$0<\frac{3{\kappa }_3}{{\kappa }_2}\left(\frac{{\kappa }_1}{{\kappa }_2}-4\right)<1.$$

$${\dot{\rho }}_{n,m}^{(p)}=$$

$$+A_{n+1}^{(p)}{\rho }_{n+1,m}^{(p)}+B_{n-2}^{(p)}{\rho }_{n-2,m}^{(p)}+C_{n+3}^{(p)}{\rho }_{n+3,m}^{(p)}$$

$$+A_{m+1}^{(p)}{\rho }_{n,m+1}^{(p)}+B_{m-2}^{(p)}{\rho }_{n,m-2}^{(p)}+C_{m+3}^{(p)}{\rho }_{n,m+3}^{(p)},$$

$${\Delta }_{n,m}^{(p)}=$$

$$G_n^{(p)}=$$

$$+{\kappa }_3[(n+p)(n+p-1)(n+p-2)+n(n-1)(n-2)]\},$$

$$A_{n+1}^{(p)}=$$

$$B_{n-2}^{(p)}=$$

$$C_{n+3}^{(p)}=$$

$${\rho }_{0,0}^{(2)}=-iJ2\sqrt{2}{\left(2{\kappa }_1+14{\kappa }_2\right)}^{-1}{\rho }_{1,0}^{(1)},$$

$${\rho }_{1,0}^{(1)}=\frac{iJ\sqrt{2}\left(P_1^2-P_0{P}_2\right)}{\left(2{\kappa }_1+13{\kappa }_2\right)}.$$

$$P\left(\varphi \right)=\frac{1}{2\pi }+\frac{2J^2\left(P_1^2-P_0{P}_2\right)\cos \left(2\varphi \right)}{\pi \left({\kappa }_1+7{\kappa }_2\right)\left(2{\kappa }_1+13{\kappa }_2\right)}.$$

$$\dot{\alpha }$$

$$\dot{\varphi }=$$

$$\dot{r}=$$

$$-\frac{3{\kappa }_3}{32}r\left(5R^4+10R^2{r}^2+r^4\right)-JR\sin \varphi .$$

$$\dot{\varphi }\approx$$

$$\dot{r}\approx$$

Full text

I Introduction

The last few years has seen rapid progress in engineering and probing the properties of nonlinear oscillators in the quantum regime Dykman (2012); Ong et al. (2011, 2013); Wang et al. (2019). This has stimulated renewed theoretical interest in the properties of such systems, both at the level of individual oscillators and for more complex many-body realizations Zhang et al. (2017); Puri et al. (2017); Dykman et al. (2018); Rota et al. (2019). Of particular interest are phenomena, such as synchronization, which result from an interplay between nonlinearity and nonequilibrium features arising from a combination of gain and loss processes. In the classical regime, the standard paradigm for synchronization involves a nonlinear oscillator in a limit-cycle state which has a well-defined amplitude, but no preferred phase Pikovsky et al. (2001); Acebrón et al. (2005). Such oscillators have a tendency to adjust their rhythm to match either a weak external drive or that of other oscillators to which they may be coupled, typically leading to the emergence of a definite phase (or relative phase) for the oscillations. Although synchronization has been studied for a very long time in classical oscillators Pikovsky et al. (2001), the systematic study of this behavior in quantum oscillators outside the regime where semiclassical approximations work well Cresser et al. (1982) is quite recent Zhirov and Shepelyansky (2006); Manzano et al. (2013); Ludwig and Marquardt (2013); Mari et al. (2013); Lee and Sadeghpour (2013); Walter et al. (2014).

Studies of synchronization in quantum oscillators have explored issues such as the variety of ways in which the behavior differs from what is found in the semiclassical limit Lee and Sadeghpour (2013); Lörch et al. (2016); Davis-Tilley and Armour (2016); Lörch et al. (2017), how best to quantify the degree of synchronization in the quantum regime Mari et al. (2013); Lee and Sadeghpour (2013); Ameri et al. (2015); Hush et al. (2015), and the relationship between synchronization and entanglement Mari et al. (2013); Manzano et al. (2013); Lee et al. (2014); Ameri et al. (2015); Li et al. (2015); Davis-Tilley and Armour (2016); Galve et al. (2017); Roulet and Bruder (2018a). Much of the work on quantum synchronization has involved simple models such as the quantum van der Pol (QvdP) oscillator Lee and Sadeghpour (2013); Walter et al. (2014); Lörch et al. (2016, 2017) and spin-1 particles Roulet and Bruder (2018a, b); Koppenhöfer and Roulet (2019), although a range of other systems have also been considered including atomic ensembles Xu et al. (2014); Zhu et al. (2015) and optomechanical systems Mari et al. (2013); Ludwig and Marquardt (2013); Weiss et al. (2016). Significant efforts have also been devoted to proposing ways in which the behavior could be probed in experiment using systems such as trapped ions Lee and Sadeghpour (2013); Hush et al. (2015) or superconducting circuits Nigg (2018).

Although studies of synchronization in the quantum regime have employed a wide range of model systems, they have focussed (with occasional exceptions Zhirov and Shepelyansky (2008)) on systems whose dynamics is essentially just a limit-cycle. Here we instead explore how weak coupling generates synchronization, in the form of particular phase preferences, in a quantum oscillator which has a more complex bistable dynamics. We do this by proposing a minimal model for a quantum oscillator that displays a limit-cycle state as well as a low occupation-number state (in which the oscillator simply fluctuates about the origin) and can be tuned to a bistable regime in which both of these states coexist.

Our model involves only dissipative processes in which photons (quanta) are lost or gained, ensuring that an isolated oscillator never has a preferred phase. The key ingredient of the model responsible for generating bistability is a two-photon gain process. This is balanced by two channels of photon loss in which either a single photon or three photons are annihilated in the oscillator. Different dynamical states of the oscillator (low occupation-number regime, limit-cycle and bistability) are achieved by tuning the relative sizes of the gain and loss rates. In many ways the model is a logical extension of the much-studied QvdP oscillator which combines one-photon gain and two-photon loss Lee and Sadeghpour (2013). However, the QvdP model always displays a limit-cycle, albeit one whose size depends on the ratio of loss and gain rates.

We investigate in detail the phase synchronization that occurs when two of the bistable model oscillators are coupled via a weak photon exchange process. This leads to a rich range of behavior in the relative phase distribution, with a different pattern of phase preferences emerging depending on the underlying dynamical states of the oscillators. When the gain is strong and the oscillators are in limit-cycle states, maximal values of the relative phase distribution form at [math] and $\pi$, matching what is usually found for such states(as seen, e.g.  in the QvdP model Lee and Sadeghpour (2013)). For oscillators with low occupation numbers and no limit-cycle, however, relative phases of $\pi/2$ and $3\pi/2$ are preferred instead, a result which we argue can be understood as a result of the two-photon gain. All four peaks emerge simultaneously for a small, intermediate parameter regime. Furthermore, this $\pi/2$-periodic behavior is strongest when the bistability is most pronounced. Interestingly, and in strong contrast to the QvdP oscillator, we find that the strength of the synchronization in the limit-cycle state does not increase with increasing photon numbers.

The outline of the rest of this paper is as follows. In Sec. II, we start by introducing our bistable quantum oscillator model and exploring its steady-state and dynamical properties. Then we investigate the behavior of two coherently coupled oscillators in Sec. III, focussing in particular on the way in which the pattern and strength of features in the relative phase distribution of the system depend on the underlying dynamical properties of each oscillator. We conclude in Sec. IV and the Appendixes provide details about aspects of the calculations employed.

II Bistable Oscillator System

II.1 321-Oscillator Model

Our oscillator model involves three dissipative processes, as illustrated in Fig. 1(a). A two-photon gain process with rate $\kappa_{2}$ drives the oscillator to higher photon numbers, whilst a one-photon loss process damps it at a rate $\kappa_{1}$; an additional three-photon loss process at rate $\kappa_{3}$ is included to stabilise the system, ensuring that it has a steady state for any strength of the gain. The master equation for a single oscillator in the interaction picture is given by Lindblad (1976); Gardiner and Zoller (2000); Lee and Sadeghpour (2013),

[TABLE]

where $\hat{a}$ is the oscillator lowering operator and we have defined $\mathcal{D}[\hat{C}]\left(\rho\right)=\hat{C}\rho\hat{C}^{\dagger}-\tfrac{1}{2}\{\hat{C}^{\dagger}\hat{C},\rho\}$. Our model makes an interesting contrast with the QvdP oscillator Lee and Sadeghpour (2013); Walter et al. (2014), where one-photon gain is balanced by two-photon loss. The presence of a non-linear gain process in our model leads to important features, such as bistability, not seen for the QvdP.

The steady-state properties are readily found by exploiting the fact that the system is purely dissipative, so that the dynamics of the diagonal and off-diagonal matrix elements of the density operator in the number (Fock) basis are decoupled Scully and Zubairy (1997); Scarlatella et al. (2019). The master equation can therefore be rewritten as a set of $k$ equations

[TABLE]

where $\rho_{n}^{(k)}=\langle n|\rho|n+k\rangle$, with $|n\rangle$ the $n$-th number state, and $\mathcal{M}^{(k)}$ a matrix. For the diagonal elements, writing out Eq. (1) explicitly leads to the coupled set of equations

[TABLE]

for the probabilities, $P_{n}=\langle n|\rho|n\rangle$, with

[TABLE]

from which the form of $\mathcal{M}^{(0)}$ follows. In the steady state, the off-diagonal terms $\rho^{(k\neq 0)}$ all go to zero and the eigenvector of $\mathcal{M}^{(0)}$ with zero eigenvalue gives the $P_{n}$ distribution.

II.2 Phase diagram

The steady-state of the oscillator can be characterised by the behavior of the $P_{n}$ distribution along with the Wigner distribution Gardiner and Zoller (2000), $W\left(\alpha_{r},\alpha_{i}\right)$. Figures 1(b)-(e) show how the state of the system evolves as the ratio $\kappa_{1}/\kappa_{2}$ is changed for a small (fixed) value of $\kappa_{3}/\kappa_{2}$. When the non-linear gain dominates ($\kappa_{1}/\kappa_{2}\ll 1$) the oscillator is driven to large phonon numbers with an almost Gaussian $P_{n}$ distribution centered at a value $\langle n\rangle=\sum_{n}nP_{n}\gg 1$ [see Fig. 1(b)]. The corresponding Wigner distribution reveals a distribution with a ring of maxima [Fig. 1(c)], we classify this as a limit-cycle (LC) state, as it has a well-defined average amplitude, but no preferred phase Lee and Sadeghpour (2013). In the opposite limit of dominant loss ($\kappa_{1}/\kappa_{2}\gg 1$), the oscillator is damped to the lowest photon number states, leading to a sharp peak in the $P_{n}$ distribution at $n=0$. In this regime the Wigner distribution displays a single maximum at the origin [Fig. 1(e)] and we call this a fixed point (FP) state. In between these limits we find bistability (B) where features from both the LC and FP states can be found in the Wigner distribution [see Fig. 1(d)] and two peaks of similar size feature in the $P_{n}$ distribution Rodrigues et al. (2007).

The bimodality of the $P_{n}$ distribution is captured by a sharp peak in the second moment Nielsen and Kittel (2006) $\mu_{(2)}=\langle n^{2}\rangle-\langle n\rangle^{2}$, as shown in Fig. 1(b). The distributions with the highest values of $\mu_{(2)}$ are found to be those with the most pronounced bistability, i.e. with two peaks of comparable area that are separated by a significant gap.

We use a standard mean-field (or semiclassical) approach to understand the origin of the oscillator’s dynamical states (the details are described in Appendix A). Two physically relevant mean-field solutions for the average photon number are found: a zero photon solution $n_{0}=0$ (corresponding to the fixed point), and a nonzero solution $n_{+}=\tfrac{\kappa_{2}}{3\kappa_{3}}\left[1+\sqrt{1+\tfrac{3\kappa_{3}}{\kappa_{2}^{2}}\left(4\kappa_{2}-\kappa_{1}\right)}\right]$ (corresponding to the limit-cycle). Linear stability analysis reveals that $n_{0}$ is stable for $\kappa_{1}>4\kappa_{2}$, whilst $n_{+}$ is both non-negative and stable for $3\kappa_{1}\kappa_{3}<12\kappa_{3}\kappa_{2}+\kappa_{2}^{2}$. Hence the mean-field approach predicts a region of bistability associated with the coexistence of these stable solutions in the parameter regime

[TABLE]

These predictions are compared with the behavior of the $P_{n}$ distribution in Fig. 1(b). The $n_{+}$ solution matches the average photon number $\langle n\rangle$ very well deep inside the limit-cycle regime, where photon numbers are large. However, the mean-field calculation doesn’t describe the extent of the bistable region very accurately, which is not surprising given the strong quantum non-linearity in this regime.

The behavior of the oscillator as a function of both $\kappa_{1}/\kappa_{2}$ and $\kappa_{3}/\kappa_{2}$ is summarised in Fig. 2. Classifying the oscillator state as either FP, B, or LC, based on the corresponding Wigner function leads to the ‘phase diagram’ shown in Fig. 2(a) whilst the behavior of the average occupation number is shown in Fig. 2(b). A sharp transition occurs between the FP and LC states via a bistable region which is narrow and very well defined (e.g. via a clear peak in the second moment) for $\kappa_{3}/\kappa_{2}\ll 1$ as $\kappa_{1}/\kappa_{2}$ is reduced. Average photon numbers tend to decrease as $\kappa_{3}/\kappa_{2}$ is increased, and although we can still use the Wigner distribution to distingusih the different states, the differences between them become much less distinct.

II.3 Dynamical Properties

We now turn to the dynamical properties of the system to understand whether the bistable states we have identified are also metastable in the sense that they display slow switching between the two states. To investigate this, we start by calculating how the eigenvalues with largest real parts (i.e. least negative, but non-zero) of the system behave Risken et al. (1987); Nirmal et al. (1996); Macieszczak et al. (2016).

The largest eigenvalue of the matrix $\mathcal{M}^{(k)}$ in Eq. (2) can be used to obtain the slowest timescale associated with the dynamics of $\rho^{(k)}$, $\tau_{k}=-1/{\rm{Re}}(\lambda^{(k)}_{0})$, examples of which are shown in Fig. 3(a). The $k=0$ case, $\tau_{0}$, describes the relaxation of the diagonal elements and it becomes very large for a range of $\kappa_{1}/\kappa_{2}$. The other times-scales $\tau_{1,2}$, describe the relaxation of phase preferences in the system. Although they never become as large as the peak values of $\tau_{0}$ and display no obvious signature of the bistability, they do change significantly as the system evolves from FP to LC states, becoming orders of magnitude larger in the latter case.

We can use the emergence of a single very slow timescale as a criterion to identify metastability in the system Macieszczak et al. (2016). We calculate the ratio of the differences between the three largest eigenvalues

[TABLE]

as shown in Fig. 3(b), and take (somewhat arbitrarily) a separation of $\lambda_{1}$ and $\lambda_{0}$ an order of magnitude larger than $\lambda_{2}$ and $\lambda_{1}$, i.e. $M\geq 10$, as a threshold for metastability. The parameter range identified via this criterion is, unsurprisingly, a small subset of that labelled on the basis of the steady-state as bistable in Fig. 2(a). Nevertheless, the peaks in $M$ and the second moment plotted in Fig. 2(a) match up well.

We can get more insight into the dynamical properties of the system by looking at quantum jump trajectories Plenio and Knight (1998), obtained by unravelling the master equation. The system is evolved in time with the non-Hermitian Hamiltonian $H_{MC}=-\frac{i}{2}\left(\kappa_{1}\hat{a}^{\dagger}\hat{a}+\kappa_{2}\hat{a}^{2}(\hat{a}^{\dagger})^{2}+\kappa_{3}(\hat{a}^{\dagger})^{3}\hat{a}^{3}\right)$. For each short time step, quantum jumps of three different kinds (one-photon loss, two-photon gain, or three-photon loss process) can occur with a probability that depends on the state of the system (e.g. the probability of two photon gain occurring over the interval $\delta t$ is given by $\kappa_{2}\langle\psi(t)|\hat{a}^{2}(\hat{a}^{\dagger})^{2}|\psi(t)\rangle\delta t$ where $|\psi(t)\rangle$ is the state of the oscillator).

The frequency of the different jump processes is illustrated in Fig. 4 for sample trajectories obtained for parameters corresponding to the three different states of the system (LC, FP and B). Within the LC state the oscillator has a large number of photons and consequently displays a high level of activity (i.e. frequent jumps). In contrast, all jump processes are strongly suppressed (and the three-photon loss especially so) in the FP state because of the very low occupation numbers in this regime. Within the region which is both bistable and metastable (based on the behavior of the eigenvalues), the oscillator switches between periods in which it exhibits high and low levels of activity; the switching continues indefinitely, never settling into one state or the other. This intermittency in the dynamics of the trajectories is consistent with our interpretation of this regime as metastable Macieszczak et al. (2016).

III Synchronization of Coupled Oscillators

We now explore how phase ordering and synchronization occurs when two of these oscillators are coupled together weakly. For simplicity, we consider two identical oscillators and assume a coherent (photon-exchange) interaction of the form Lee and Sadeghpour (2013); Walter et al. (2014); Hush et al. (2015)

[TABLE]

where $\hat{a}_{j}$ is a lowering operator for oscillator $j$ and $J$ is the strength of the coupling. The master equation of the coupled system is given by

[TABLE]

where the dissipation terms follow from Eq. (1): $\mathcal{L}_{j}\rho=\kappa_{1}\mathcal{D}[\hat{a}_{j}](\rho)+\kappa_{2}\mathcal{D}[(\hat{a}^{\dagger}_{j})^{2}](\rho)+\kappa_{3}\mathcal{D}[\hat{a}^{3}_{j}](\rho)$.

Classical limit-cycle oscillators typically synchronize when they are coupled together weakly, developing a preference for one or more relative phase values Pikovsky et al. (2001); Kuznetsov et al. (2009). In the quantum regime, the phase states $|\varphi\rangle=\sum_{n}{\rm{e}}^{i\varphi}|n\rangle/2\pi$, can be used to construct a relative phase probability distribution Barnett and Pegg (1990); Luis and Sánchez-Soto (1996); Barak and Ben-Aryeh (2005); Lörch et al. (2017); Davis-Tilley and Armour (2016); Hush et al. (2015)

[TABLE]

where $\phi=\varphi_{1}-\varphi_{2}$ is the relative phase of the two oscillators. If the two oscillators are uncoupled, their phases are independent and the relative phase distribution is uniform, $P(\phi)=\tfrac{1}{2\pi}$. In the following we use the relative phase distribution to explore the impact of coupling on the behavior of our oscillator model in each of the dynamical states it displays.

Since we are interested in the effect of a very weak coupling between the oscillators, we use a perturbation method Lörch et al. (2016); Davis-Tilley and Armour (2016); Koppenhöfer and Roulet (2019) to calculate the relative phase distribution, details of which are described in Appendix B. Figure 5(a) shows how the relative phase distribution evolves as a function of $\kappa_{1}/\kappa_{2}$ and the system passes from FP to LC via the bistable region. Interestingly, the system displays $\pi$-periodic phase distributions with peaks which reach a similar size in both the LC and FP regimes. However, the location of the peaks are different in the two cases and in between (the bistable region) the phase distribution appears to flatten.

A more detailed picture of the phase behavior is obtained by looking at the Fourier coefficients of the relative phase distribution, $F_{k}=\text{Re}[\sum_{n,m=0}^{\infty}\langle n+k,m|\rho|n,m+k\rangle]$. For our coupled oscillator system, these coefficients are non-zero for even $k$ and typically get smaller very rapidly with increasing $k$. The two most important coefficients $F_{2}$ and $F_{4}$ are shown in Fig. 5(b) (calculated to second and fourth order in $J$ respectively). The $\pi$-periodicity of the distribution apparent in Fig. 5(a) stems from the fact that the magnitude of $F_{2}$ is almost everywhere much larger than that of $F_{4}$. However, as the system transitions from LC to FP (via the bistablity) $F_{2}$ changes sign to produce the shift in the locations of the peaks; as $F_{2}$ passes through zero $F_{4}$ dominates, giving rise to an unusual $\pi/2$ periodic distribution.

The different $\pi$-periodic patterns that arise in the phase distribution can be understood using simple arguments that exploit the specific characteristics of the system in the LC and FP states. Well-within the LC regime photon occupation numbers are large and semiclassical approaches work well as Fig. 1(b) illustrates for the single oscillator case. A straightforward calculation described in Appendix C recovers the preference for relative phases of [math] and $\pi$ which is generic for coherently coupled limit-cycle oscillators Lee and Sadeghpour (2013); Davis-Tilley and Armour (2016).

The pattern of peaks in the relative phase distribution at $\pi/2$ and $3\pi/2$ seen for the FP can be understood by focussing on the limit where $\kappa_{1}/\kappa_{2}\gg 1$. In this regime photon occupation numbers are very small, the value of $\kappa_{3}$ becomes irrelevant, and we can simplify the perturbation theory calculation by assuming only the three lowest Fock states of the oscillators have non-negligible occupations (see Appendix B for details). In this limit we find

[TABLE]

where $f(\kappa_{1},\kappa_{2})$ is a (positive-valued) function of the two rates, whilst $P_{0}$, $P_{1}$ and $P_{2}$ are the occupation probabilities of the three lowest Fock states for the corresponding uncoupled oscillator. Clearly the behavior is always $\pi$-periodic, but the positions of the maxima depend on the occupation probabilities which in turn depend on the details of how the system is driven and damped in this regime. In our oscillator model, the two-photon driving gives a boost to $P_{2}$, generating a steady-state with a photon distribution where $P_{2}P_{0}/P_{1}^{2}>1$ which leads to the peaks at $\pi/2$ and $3\pi/2$ seen in Fig. 5. Thus the fact that our oscillator model displays a $\pi$-periodic relative phase distribution with peaks at $\pi/2$ and $3\pi/2$ can be seen as being due to the particular form of the lowest order nonlinearity in the system, two-photon gain, which shapes the steady-state number distribution.

Applying a very similar analysis to oscillators with different gain and loss processes casts light on the different ways relative phase preferences can develop (within the same low-occupation number limit). In particular, for an oscillator coupled to a thermal bath [with one-photon gain and loss related by a ratio of rates $\overline{n}/(\overline{n}+1)$, with thermal occupation number $\overline{n}$] and the QvdP oscillator Lee and Sadeghpour (2013) (one-photon gain and two-photon loss) one obtains expressions like that in Eq. (10). In each case $f$ takes a form that is different in detail, but the sign of $P_{1}^{2}-P_{0}P_{2}$ still determines the location of the peaks. For the QvdP oscillator in the strongly damped limit $P_{2}P_{0}<P_{1}^{2}$, hence one finds the same phase behavior as in the large photon number regime Lee and Sadeghpour (2013) (peaks at [math] and $\pi$). In contrast, for a thermal oscillator $P_{2}P_{0}=P_{1}^{2}$ and so no phase preference is expected.

Finally, we look at how the phase behavior of the system varies over the full phase diagram. Figure 6 shows how the components $F_{2}$ and $F_{4}$ behave for a range of $\kappa_{3}/\kappa_{2}$ as well as $\kappa_{1}/\kappa_{2}$. The transition between negative and positive values of $F_{2}$ is sharpest, and the corresponding peak in $F_{4}$ strongest, in the region where $\kappa_{3}/\kappa_{2}\ll 1$, which is where the system also displays metastability.

Figure 6 also reveals a somewhat surprising feature of the phase behavior. The strongest phase preferences are not associated with the limit-cycle regime. Furthermore, deep within the limit-cycle regime the magnitude of $F_{2}$ starts to saturate, even as the average photon number continues to grow, i.e. as $\kappa_{3}/\kappa_{2}$ is further reduced in the bottom left quadrant of Fig. 6(a). This is in contrast to other quantum limit-cycle oscillators, such as the QvdP oscillator, for which phase synchronization effects are enhanced by increasing the photon number Lee and Sadeghpour (2013). We can gain some insight into why this occurs from the phase dynamics of an uncoupled oscillator. At least in the semiclassical regime of large photon numbers, we expect to see stronger synchronization effects emerge when coupling is introduced in systems where the phase diffusion is weaker Pikovsky et al. (2001); Davis-Tilley and Armour (2016).

Using an approximate analytic approach (see Appendix D for details), we find that the phase diffusion rate is simply proportional to $\kappa_{2}$ in the regime where gain dominates over losses and photon numbers are large. This contrasts with the behavior of the QvdP oscillator where the phase diffusion is $\propto 1/\langle n\rangle$, which is why for this system synchronization effects get stronger as photon numbers are increased.

IV Conclusions

We have introduced a simple oscillator model with a two-photon gain process balanced by one- and three-photon losses that can be used to engineer a bistable oscillator state. The bistability occurs when the gain/loss rates are tuned between regimes where the oscillator displays a limit-cycle state characterised by a large amplitude (but no preferred phase) and a fixed-point state where occupation numbers are low. Quantum trajectory simulations show clear evidence of metastability in the bistable region, signalled by intermittency in the frequency of the quantum jumps.

When two such oscillators are coupled together weakly their relative phase distribution displays a rich pattern of behavior. Whilst the system has the usual predominantly $\pi$-periodic distribution with peaks at [math] and $\pi$ within the limit-cycle regime, in the limit of low occupation numbers the peaks appear at $\pi/2$ and $3\pi/2$. In between these two regimes, where the bistability arises, the distribution can instead be $\pi/2$-periodic. We identify the presence of two-photon gain as the key factor giving rise to this unusual behavior.

Our goal in this work has been to explore the properties of the relative phase distribution and its connection to the underlying oscillator dynamics in the simplest possible model system displaying limit-cycles and bistability. We plan to address the question of how the complex behavior we uncovered could be realised, and detected, in future work. It would also be interesting to investigate how the dynamics is affected by strong couplings Ishibashi and Kanamoto (2017) and the patterns of phase behavior that arise in systems where more than two oscillators are coupled together Ludwig and Marquardt (2013); Lee and Sadeghpour (2013); Davis-Tilley et al. (2018).

Acknowledgements

This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/N50970X/1] through a studentship held by MJ. W.L. acknowledges support from the UKIERI-UGC The- matic Partnership No. IND/CONT/G/16-17/73, EP- SRC Grant No. EP/M014266/1 and EP/R04340X/1.

Appendix A Mean-Field Analysis

An equation of motion for the expectation value of the annihilation operator, $\frac{d}{dt}\left\langle\hat{a}\right\rangle=\text{Tr}\left[\hat{a}\dot{\rho}\right]$, is found from the master equation [Eq. (1)]

[TABLE]

In a mean-field (or semiclassical) approach we break the correlations between operators such that $\langle\hat{A}\hat{B}\rangle=\langle\hat{A}\rangle\langle\hat{B}\rangle$. We choose to carry out this approximation after first normal ordering the operators.

Making the substitution $\langle a\rangle=r{\rm{e}}^{i\varphi}$ where $r$ and $\varphi$ are classical amplitude and phase variables, allows us to rewrite Eq. (11) as

[TABLE]

Evaluating the real and imaginary parts separately, leads to $\dot{\varphi}=0$ and

[TABLE]

with steady state solutions $n_{0}=0$ and

[TABLE]

These solutions correspond to the average photon number and so are subject to the constraints of being both positive and real. The stability of the solutions can be determined, e.g. through the properties of the relevant Jacobian.

The negative branch, $n_{-}$, is never both physical and stable, the zero photon solution, $n_{0}$, is stable for $\frac{\kappa_{1}}{\kappa_{2}}>4$, and the positive branch, $n_{+}$, is stable whenever it is physical, i.e. for $\frac{\kappa_{1}}{\kappa_{2}}<4+\frac{\kappa_{2}}{3\kappa_{3}}$. This mean-field calculation results in two stable solutions for the average photon number and therefore a predicted bistability in the photon number for the parameter regime

[TABLE]

Appendix B Perturbation theory

B.1 General Method

Perturbation theory provides a convenient way of calculating the way in which the relative phase distribution behaves for weak coupling. The steady-state of the uncoupled ($J=0$), two-oscillator system [Eq. (7)] only has diagonal terms. Treating the coupling as a perturbation Lörch et al. (2017); Koppenhöfer and Roulet (2019) allows us to calculate the terms in the first off-diagonal as a function of the uncoupled oscillator terms. Each subsequent off-diagonal can, in turn, be calculated from previous ones Davis-Tilley and Armour (2016).

Writing Eq. (7) in the number state basis, with $\rho_{n,m}^{(p)}=\langle n+p,m|\rho|n,m+p\rangle$, leads to a set of simultaneous equations

[TABLE]

with

[TABLE]

In the steady-state, this reduces to sets of simultaneous equations with the coupling term, $\Delta^{(p)}_{n,m}$ coupling together terms with different $p$-values. The zeroth-order terms are the diagonal ($p=0$) elements, the uncoupled probabilities $\rho_{n,m}^{(0)}=P_{n}P_{m}$, which necessarily sum to unity. The first-order terms are obtained by substituting the zeroth-order terms into the expression for $\Delta_{n,m}^{(p)}$, leading to non-zero contributions for $p=1$ and the process is continued to higher order in $J$ recursively.

The first-order terms obey the relation $\rho_{m,n}^{(1)}=-\rho_{n,m}^{(1)}$ and hence sum to zero Davis-Tilley and Armour (2016), which means that they make no contribution to the relative phase distribution [Eq. (LABEL:eq:relph)] since it depends on sums of the off-diagonal elements. The sum of the $p=2$ terms, however, is real and finite and so does contribute resulting in a $\pi$-periodic relative phase distribution. Continuing to higher orders, we find that all of the odd-$p$ terms sum to zero, and so only the even $p$ sums contribute to the relative phase distribution. In particular, the $p=4$ terms, lead to a ${\pi}/{2}$-periodic contribution which can dominate the phase distribution when the $\pi$-periodic terms vanish.

B.2 Low Occupation-Number Regime

This calculation can be simplified and solved analytically in the limit of very low photon numbers. We proceed by assuming only the lowest three photon states are appreciably occupied, i.e. $P_{n>2}=0$, and hence truncate the state space to include only $|0\rangle$, $|1\rangle$, and $|2\rangle$. Due to the size of the Hilbert space, only a single term contributes to the relative phase distribution, $P(\phi)=\frac{1}{2\pi}+\frac{1}{\pi}\text{Re}\left[e^{2i\phi}\rho_{0,0}^{(2)}\right]$. In the steady-state, Eq. (16) with $p=2$ leads to

[TABLE]

using the relation $\rho_{m,n}^{(1)}=-\rho_{n,m}^{(1)}$. Equation (16) with $p=1$ gives

[TABLE]

This results in the relative phase distribution

[TABLE]

This is a $\pi$-periodic distribution and the position of the peaks is determined by the steady state of the uncoupled oscillators. The two-photon driving in our model ensures $P_{2}P_{0}>P_{1}^{2}$, which leads to peaks at $\pi/2$ and $3\pi/2$.

Appendix C Mean-Field Analysis of Coupled Oscillators

The same mean-field procedure discussed above in Appendix A can be applied to the case of two coupled oscillators with $\alpha=\left\langle\hat{a}_{1}\right\rangle$ and $\beta=\left\langle\hat{a}_{2}\right\rangle$, leading to the equations of motion

[TABLE]

Changing to polar co-ordinates, with the definitions $\alpha=r_{1}e^{i\phi_{1}}$ and $\beta=r_{2}e^{i\phi_{2}}$, and then rewriting the equations in terms of the sum-and-difference variables $r=r_{1}-r_{2}$, $R=r_{1}+r_{2}$, and $\phi=\phi_{1}-\phi_{2}$ we find

[TABLE]

For large photon occupation numbers (as is the case in the limit-cycle regime) and for weak couplings ($J/\kappa_{2}\ll 1$), we have $r\ll R$ and can approximate these equations as

[TABLE]

The weakness of the coupling and the small size of the ratio $r/R$ leads to a separation of timescales with $\phi$ relaxing much more slowly than $r$. Adiabatic elimination of $r$ leads to the simple relation $\dot{\phi}=-{\partial U(\phi)}/{\partial\phi}$ with the pseudo-potential Kuznetsov et al. (2009)

[TABLE]

This potential predicts relative phase preferences for [math] and $\pi$, its stable minima.

Appendix D Phase Diffusion

In this Appendix we return to the case of a single oscillator and obtain an estimate for the phase diffusion rate in the semiclassical limit where photon numbers are large. In the semiclassical regime at least, the strength of phase synchronization in coupled oscillators is determined by competition between the coupling and the rate of phase diffusion in the individual oscillators, with slower phase diffusion leading to stronger phase preferences Davis-Tilley and Armour (2016); Pikovsky et al. (2001)

The phase distribution for a single oscillator takes the form Knight and Gerry (2004); Barnett and Pegg (1990); Davis-Tilley and Armour (2016)

[TABLE]

with $\Phi^{(k)}=\sum_{n=0}^{\infty}\rho^{(k)}_{n}$, where $\rho^{(k)}_{n}=\langle n|\rho|n+k\rangle$. Although the behavior is in general quite complex, we can obtain a simple approximate description in the semiclassical limit where the density matrix is tightly peaked around a large average photon occupation number Scully and Zubairy (1997); Davis-Tilley and Armour (2016).

Using (1) [and the notation introduced in (16)], we can obtain the equation of motion for $\Phi^{(k)}$

[TABLE]

In the semiclassical limit, i.e. the strong gain regime where $\gamma=\kappa_{1}/\kappa_{2}\ll 1$ and $\Gamma=\kappa_{3}/\kappa_{2}\ll 1$, the photon number saturates to a large value $\langle n\rangle\simeq 2\kappa_{2}/(3\kappa_{3})$ [see Eq. (14)]. We proceed by assuming we can replace $n$ by $\langle n\rangle$, and expand the square-roots appearing in $A_{n}^{(k)}$, $B_{n}^{(k)}$, etc, treating $1/\langle n\rangle$ together with $\gamma$ and $\Gamma$ as small quantities Scully and Zubairy (1997); Davis-Tilley and Armour (2016); Schieve and McGowan (1993). This leads to the simplified equation

[TABLE]

Hence to leading order the relaxation timescale for the $k$-th component, $\Phi^{(k)}$, is simply proportional to $1/\kappa_{2}$, $\tau_{k}^{LC}\kappa_{2}\simeq 4/(5k^{2})$.

Notice that the semiclassical approximation here assumes all the off-diagonal elements ($\rho_{n}^{(k)}$) decay at the same rate, hence the decay time for $k=1$ is also an approximate linewidth for the oscillator Scully and Zubairy (1997). The slowest timescales $\tau_{k}$ associated with the matrices $\mathcal{M}^{(k)}$ in Eq. (2) are found to plateau in the limit of $\kappa_{1}/\kappa_{2}\ll 1$ and $\kappa_{3}/\kappa_{2}\ll 1$ [illustrated on the left-hand side of Fig. 3(a)] and numerically we find $\tau_{1}\simeq 0.8/\kappa_{2}$ and $\tau_{2}\simeq 0.2/\kappa_{2}$ in this regime, matching up very well with $\tau_{1}^{LC}$ and $\tau_{2}^{LC}$, respectively.

Finally, using the definition Eq. 33 and Eq. 35, we see that the phase distribution obeys a diffusion equation Davis-Tilley and Armour (2016); Schieve and McGowan (1993)

[TABLE]

This is very different to what is found obtained in a similar calculation for the QvdP oscillator (or indeed the laser Scully and Zubairy (1997)) where the diffusion constant is $\propto 1/\langle n\rangle$. Hence for such systems phase diffusion gets weaker (and the linewidth narrower), so that synchronization effects get stronger, as the photon number increases Lee and Sadeghpour (2013).