Theory of coupled parametric oscillators beyond coupled Ising spins

TL;DR

This paper investigates the complex dynamics of two coupled parametric oscillators beyond the simple Ising spin analogy, revealing new phases and universal behavior through analytical and numerical methods.

Contribution

It provides a comprehensive theoretical analysis of coupled parametric oscillators with energy-conserving interactions, extending understanding beyond traditional Ising spin models.

Findings

Identification of phase boundaries and critical exponents.

Discovery of new phases with multiple attractors.

Universal character of the phase diagram independent of nonlinearity type.

Abstract

Periodically driven parametric oscillators offer a convenient way to simulate classical Ising spins. When many parametric oscillators are coupled dissipatively, they can be analogous to networks of Ising spins, forming an effective coherent Ising machine (CIM) that efficiently solves computationally hard optimization problems. In the companion paper, we studied experimentally the minimal realization of a CIM, i.e. two coupled parametric oscillators [L. Bello, M. Calvanese Strinati, E. G. Dalla Torre, and A. Pe'er, Phys. Rev. Lett. 123, 083901 (2019)]. We found that the presence of an energy-conserving coupling between the oscillators can dramatically change the dynamics, leading to everlasting beats, which transcend the Ising description. Here, we analyze this effect theoretically by solving numerically and, when possible, analytically the equations of motion of two parametric oscillators. Our main tools include: (i) a Floquet analysis of the linear equations, (ii) a multi-scale analysis based on a separation of time scales between the parametric oscillations and the beats, and (iii) the numerical identification of limit cycles and attractors. Using these tools, we fully determine the phase boundaries and critical exponents of the model, as a function of the intensity and the phase of the coupling and of the pump. Our study highlights the universal character of the phase diagram and its independence on the specific type of nonlinearity present in the system. Furthermore, we identify new phases of the model with more than two attractors, possibly describing a larger spin algebra.

Figures

Tables

Table 1. Table 1: Types of nonlinearity (pump depletion or Kerr/saturation) that we consider in the equations of motion in Eq. ( 32 ).

Nonlinearity

Pump depletion

$F_{\mathrm{NL},1}=-h\beta {\omega }_0^2\sin \left(2{\omega }_{0}t\right)x_1^3$

$F_{\mathrm{NL},2}=-h\beta {\omega }_0^2\sin \left(2{\omega }_{0}t+\varphi \right)x_2^3$

Kerr, saturation

$F_{\mathrm{NL};1}=-\beta {\omega }_0^2{x}_1^3$

$F_{\mathrm{NL};2}=-\beta {\omega }_0^2{x}_2^3$

Equations

$$\begin{array}[]{l}\ddot{x}_{1}+\omega^{2}_{0}[1+h\,\sin(\gamma t)]x_{1}+\omega_{0}g\,\dot{x}_{1}-\omega_{0}r\,\dot{x}_{2}=0\vspace{0.1cm}\\ \ddot{x}_{2}+\omega^{2}_{0}[1+h\,\sin(\gamma t)]x_{2}+\omega_{0}g\,\dot{x}_{2}+\omega_{0}r\,\dot{x}_{1}=0\end{array}\,\,.$$

$${\ddot{x}}_{+}+{\omega }_0^2[1+h\sin (\gamma t)]x_{+}+{\omega }_0(g+ir){\dot{x}}_{+}=0$$

$${\ddot{x}}_{-}+{\omega }_0^2[1+h\sin (\gamma t)]x_{-}+{\omega }_0(g-ir){\dot{x}}_{-}=0.$$

$${\ddot{y}}_{\pm }+{\omega }_0^2\left[1-\frac{{(g\pm ir)}^2}{4}+h\sin (\gamma t)\right]y_{\pm }=0.$$

$${\ddot{y}}_{\pm }+{\omega }_0^2\left[1+h\sin (\gamma t)\right]y_{\pm }=0.\text{vspace}0.2cm$$

$$y_{\pm }(t)=e^{-i\mu t}\sum _{n\in \mathbb{Z}}{A}_n{e}^{in\gamma t},$$

$$D_n(\mu )A_n+\frac{i}{2}{\omega }_0^{2}h\left(A_{n+1}-A_{n-1}\right)=0,$$

$$\left(\begin{array}[]{ccccccc}\ddots&&&\vdots&&&\reflectbox{$\ddots$}\\ &-i\,\cfrac{\omega^{2}_{0}h}{2}&D_{-1}(\mu)&i\,\cfrac{\omega^{2}_{0}h}{2}&0&0&\\ \\ \cdots&0&-i\,\cfrac{\omega^{2}_{0}h}{2}&D_{0}(\mu)&i\,\cfrac{\omega^{2}_{0}h}{2}&0&\cdots\\ \\ &0&0&-i\,\cfrac{\omega^{2}_{0}h}{2}&D_{+1}(\mu)&i\,\cfrac{\omega^{2}_{0}h}{2}&\\ \reflectbox{$\ddots$}&&&\vdots&&&\ddots\end{array}\right)\left(\begin{array}[]{c}\vdots\\ A_{-1}\\ \\ \\ A_{0}\\ \\ \\ A_{+1}\\ \vdots\end{array}\right)=0\,\,,$$

$${\rm det}\left(\begin{array}[]{cc}\omega^{2}_{0}-{(\mu+\gamma)}^{2}&i\,\cfrac{\omega^{2}_{0}h}{2}\\ -i\,\cfrac{\omega^{2}_{0}h}{2}&\omega^{2}_{0}-\mu^{2}\end{array}\right)\equiv P_{4}(\mu)=0\,\,,$$

$${\mu }_{\pm }(h,ϵ)=-{\omega }_0-\frac{ϵ}{2}$$

$$\pm \frac{1}{2}\sqrt{{ϵ}^2+4ϵ{\omega }_0+8{\omega }_0^2-2{\omega }_0\sqrt{4{(ϵ+2{\omega }_0)}^2+{\omega }_0^2{h}^2}}.$$

$${\mu }_{\pm }(h,ϵ)\simeq -{\omega }_0-\frac{ϵ}{2}\pm \frac{i}{2}\sqrt{{\left(\frac{{\omega }_{0}h}{2}\right)}^2-{ϵ}^2}.$$

$$x_{\pm }(t)\sim e^{\mp i\frac{{\omega }_{0}rt}{2}}{e}^{\frac{{\omega }_{0}t}{2}\left(\sqrt{{\left(\frac{h}{2}\right)}^2-{\left(\frac{ϵ}{{\omega }_0}\right)}^2}-g\right)}\cos (\gamma t/2).$$

$$\begin{array}[]{l}\ddot{x}_{1}+\omega^{2}_{0}\left[1+h\,\sin(\gamma t-\phi/2)\right]x_{1}+\omega_{0}g\,\dot{x}_{1}-\omega_{0}r\,\dot{x}_{2}=0\vspace{0.1cm}\\ \ddot{x}_{2}+\omega^{2}_{0}\left[1+h\,\sin(\gamma t+\phi/2)\right]x_{2}+\omega_{0}g\,\dot{x}_{2}+\omega_{0}r\,\dot{x}_{1}=0\,\,,\end{array}$$

$$\begin{array}[]{l}\ddot{x}_{1}+\omega^{2}_{0}\left[1+h\,\sin(\gamma t)\cos(\phi/2)\right]x_{1}\vspace{0.1cm}\\ -h\omega_{0}^{2}\sin(\phi/2)\,\cos(\gamma t)\,x_{1}+\omega_{0}g\,\dot{x}_{1}-\omega_{0}r\,\dot{x}_{2}=0\vspace{0.1cm}\\ \ddot{x}_{2}+\omega^{2}_{0}\left[1+h\,\sin(\gamma t)\cos(\phi/2)\right]x_{2}\vspace{0.1cm}\\ +h\omega_{0}^{2}\sin(\phi/2)\,\cos(\gamma t)\,x_{2}+\omega_{0}g\,\dot{x}_{2}+\omega_{0}r\,\dot{x}_{1}=0\,\,.\end{array}$$

$$\begin{array}[]{l}\ddot{x}_{+}+\omega^{2}_{0}\left[1+h\,\sin(\gamma t)\cos(\phi/2)\right]x_{+}\vspace{0.1cm}\\ +\omega_{0}(g+i\,r)\dot{x}_{+}-h\omega_{0}^{2}\sin(\phi/2)\,\cos(\gamma t)\,x_{-}=0\vspace{0.1cm}\\ \ddot{x}_{-}+\omega^{2}_{0}\left[1+h\,\sin(\gamma t)\cos(\phi/2)\right]x_{-}\vspace{0.1cm}\\ +\omega_{0}(g-i\,r)\dot{x}_{-}-h\omega_{0}^{2}\sin(\phi/2)\,\cos(\gamma t)\,x_{+}=0\,\,.\end{array}$$

$$\left(\begin{array}[]{c}x_{+}\\ x_{-}\end{array}\right)=e^{-i\mu t}\sum_{n\,\in\,\mathbb{Z}}\,e^{in\gamma t}\,\left(\begin{array}[]{c}H^{(+)}_{n}\\ H^{(-)}_{n}\end{array}\right)\,\,,$$

$$D_{\pm ,n}(\mu )H_n^{(\pm )}+i\frac{{\omega }_0{h}^2}{2}\cos \left(\frac{\varphi }{2}\right)\left(H_{n+1}^{(\pm )}-H_{n-1}^{(\pm )}\right)$$

$$-\frac{{\omega }_0{h}^2}{2}\sin \left(\frac{\varphi }{2}\right)\left(H_{n-1}^{(\mp )}+H_{n+1}^{(\mp )}\right)=0,$$

$$\left(\begin{array}[]{ccccccc}\ddots&&&\vdots&&&\reflectbox{$\ddots$}\\ &\mathbf{M}^{*}(h,\phi)&\mathbf{D}_{-1}(\mu)&\mathbf{M}(h,\phi)&0&0&\\ \\ \cdots&0&\mathbf{M}^{*}(h,\phi)&\mathbf{D}_{0}(\mu)&\mathbf{M}(h,\phi)&0&\cdots\\ \\ &0&0&\mathbf{M}^{*}(h,\phi)&\mathbf{D}_{+1}(\mu)&\mathbf{M}(h,\phi)&\\ \reflectbox{$\ddots$}&&&\vdots&&&\ddots\end{array}\right)\left(\begin{array}[]{c}\vdots\\ \mathbf{H}_{-1}\\ \\ \\ \mathbf{H}_{0}\\ \\ \\ \mathbf{H}_{+1}\\ \vdots\end{array}\right)=0\,\,,$$

$$\displaystyle\mathbf{D}_{n}(\mu)=\left(\begin{array}[]{cc}D_{+,n}&0\\ 0&D_{-,n}\end{array}\right)$$

$$\displaystyle\mathbf{M}(h,\phi)=\left(\begin{array}[]{cc}i\,\cfrac{h\omega^{2}_{0}}{2}\,\cos(\phi/2)&-\cfrac{h\omega^{2}_{0}}{2}\,\sin(\phi/2)\\ \\ -\cfrac{h\omega^{2}_{0}}{2}\,\sin(\phi/2)&i\,\cfrac{h\omega^{2}_{0}}{2}\,\cos(\phi/2)\end{array}\right)\,.$$

$$D_{\pm ,n}\simeq {\Omega }_r^2-{\left[n\gamma -\mu -\frac{{\omega }_0}{2}(ig\pm r)\right]}^2.$$

$${\mu }_{\pm ,\pm }^{(0)}(h,\gamma )\simeq \pm \frac{{\omega }_{0}r}{2}-\frac{\gamma }{2}-i\frac{{\omega }_{0}g}{2}\pm i{\Gamma }^{(0)}(h,\gamma ).$$

$${\mu }_{\pm }^{(\pm {\omega }_{0}r)}(h,\gamma )\simeq -\frac{\gamma }{2}-i\frac{{\omega }_{0}g}{2}\pm i{\Gamma }^{(\pm {\omega }_{0}r)}(h,\gamma ,r).$$

$$\begin{array}[]{l}\ddot{x}_{1}+\omega^{2}_{0}\left[1+h\left(1-\beta\,x^{2}_{1}\right)\,\sin(\gamma t)\right]x_{1}\\ \\ \hskip 99.58464pt+\omega_{0}g\,\dot{x}_{1}-\omega_{0}r\,\dot{x}_{2}=0\\ \\ \vspace{-0.2cm}\ddot{x}_{2}+\omega^{2}_{0}\left[1+h\left(1-\beta\,x^{2}_{2}\right)\,\sin(\gamma t+\phi)\right]x_{2}\\ \\ \hskip 99.58464pt+\omega_{0}g\,\dot{x}_{2}+\omega_{0}r\,\dot{x}_{1}=0\end{array}\,\,.$$

$$\frac{\partial A_R}{\partial \tilde{\tau }}=\left[\frac{\tilde{h}}{4}-\frac{1}{2}-\frac{\beta \tilde{h}}{2}\left(A_R^2+3A_I^2\right)\right]A_R+\frac{\tilde{r}}{2}{B}_R\frac{\partial A_I}{\partial \tilde{\tau }}=\left[-\frac{\tilde{h}}{4}-\frac{1}{2}+\frac{\beta \tilde{h}}{2}\left(A_I^2+3A_R^2\right)\right]A_I+\frac{\tilde{r}}{2}{B}_I$$

$$\frac{\partial B_R}{\partial \tilde{\tau }}=\frac{\tilde{h}}{4}\left[B_R\cos (\varphi )+B_I\sin (\varphi )\right]-\frac{1}{2}{B}_R-\frac{\beta \tilde{h}}{2}\left[B_R^3\cos (\varphi )+2B_I^3\sin (\varphi )+3B_R{B}_I^2\cos (\varphi )\right]-\frac{\tilde{r}}{2}{A}_R$$

$$\frac{\partial B_I}{\partial \tilde{\tau }}=\frac{\tilde{h}}{4}\left[B_R\sin (\varphi )-B_I\cos (\varphi )\right]-\frac{1}{2}{B}_I-\frac{\beta \tilde{h}}{2}\left[2B_R^3\sin (\varphi )-B_I^3\cos (\varphi )-3B_R^2{B}_I\cos (\varphi )\right]-\frac{\tilde{r}}{2}{A}_I.$$

$$\begin{array}[]{l}\ddot{x}_{+}+\omega^{2}_{0}[1+g\tilde{h}\sin(2\omega_{0}t)]x_{+}+\omega_{0}g\,\dot{x}_{+}+i\,\omega_{0}r_{0}\tilde{r}\,\dot{x}_{+}\vspace{0.1cm}\\ \hskip 28.45274pt-\cfrac{\beta g\tilde{h}\omega_{0}^{2}}{4}\,\sin(2\omega_{0}t)\,\left(x^{3}_{-}+3\,x^{2}_{+}x_{-}\right)=0\vspace{0.1cm}\\ \ddot{x}_{-}+\omega^{2}_{0}[1+g\tilde{h}\sin(2\omega_{0}t)]x_{-}+\omega_{0}g\,\dot{x}_{-}-i\,\omega_{0}r_{0}\tilde{r}\,\dot{x}_{-}\vspace{0.1cm}\\ \hskip 28.45274pt-\cfrac{\beta g\tilde{h}\omega_{0}^{2}}{4}\,\sin(2\omega_{0}t)\left(x^{3}_{+}+3\,x_{+}x_{-}^{2}\right)=0\end{array}\,\,.$$

$$\frac{{\partial }^2}{\partial t^2}{x}_{\pm }^{(0)}+{\omega }_0^2{x}_{\pm }^{(0)}=0$$

$$\frac{{\partial }^2}{\partial t^2}{x}_{\pm }^{(B)}+2\frac{\partial }{\partial \sigma }\frac{\partial }{\partial t}{x}_{\pm }^{(0)}+{\omega }_0^2{x}_{\pm }^{(B)}\pm i{\omega }_0\tilde{r}\frac{\partial }{\partial t}{x}_{\pm }^{(0)}=0,$$

$$\frac{{\partial }^2}{\partial t^2}{x}_{\pm }^{(1)}+2\frac{\partial }{\partial \tau }\frac{\partial }{\partial t}{x}_{\pm }^{(0)}+{\omega }_0^2{x}_{\pm }^{(1)}$$

$$+\tilde{h}{\omega }_0^2\sin (2{\omega }_{0}t)x_{\pm }^{(0)}+{\omega }_0\frac{\partial }{\partial t}{x}_{\pm }^{(0)}$$

Full text

1. I Introduction

2. II Linear parametric oscillators - Solution by Floquet theorem

3. II.1 In-phase pumps

4. II.2 Varying the pump phase - General derivation

5. II.2.1 Parametric resonance at $\gamma=2\Omega_{r}$

6. II.2.2 Parametric resonances at $\gamma=2\Omega_{r}\pm\omega_{0}r$

7. II.3 Linear instability regions

8. III Nonlinear case - Perturbative multiple-scale analysis

9. III.1 Equations for the long-time dynamics

10. III.2 Phase diagram of two coupled oscillators

11. III.3 Critical scaling by a three-scale analysis

12. IV Different types of nonlinearity

13. V Dissipative coupling and CIM

14. VI Conclusions

15. .1 Details on the derivation of Eq. ()

16. .2 Stability analysis of the single parametric oscillator

17. .2.1 Pump-depletion nonlinearity

18. .2.2 Kerr nonlinearity

19. .3 Oscillation threshold and beats

20. .4 Alternative form of the coupling - Multiple-scale analysis

21. .5 Expression of the Jacobian matrix of the system in Eq. ()

I Introduction

Parametric oscillations are one of the best known examples of nontrivial effect induced by a periodic drive. Over the past decades, parametric oscillators have attracted a significant attention thanks to their wide range of applications for electronic low-noise amplification Suhl (1957); Weiss (1957); Uhlir (1958); Wade and Heffner (1958); Danielson (1959). In recent years, parametric oscillators are used as generators of squeezed light Yurke (1984); Collett and Gardiner (1984); Wu et al. (1987); Lvovsky (2015), with applications in high-accuracy sensing Caves (1981); Harry (2010); Aasi and LIGO collaborators (2013); Steinlechner et al. (2013), quantum information and communication Furusawa et al. (1998); Ralph (1999a, b); Braunstein et al. (2000); Ciattoni et al. (2018); Shaked et al. (2018). They have also been studied in the context of nano- or microelectromechanical systems (NEMS or MEMS) both for practical applications and because they represent a suitable platform to analyse fundamental aspects of nonlinear dynamics Lifshitz and Cross (2003, 2009); Kenig et al. (2009a, b, 2011); Karabalin et al. (2011); Kenig et al. (2012); Salgado Sánchez et al. (2016).

A degenerate parametric oscillator is the canonical example of a period doubling instability. Due to the external periodic pump, the system can display two regimes: a stable regime, in which the oscillator is not excited, and a parametrically amplified regime, in which the oscillator oscillates at half the frequency of the pump. In the latter case, the equation of motion of the parametric oscillator admits two solutions, characterised by a relative shift of one period of the pump. A parametric oscillator is therefore the simplest example of a discrete time crystal that explicitly breaks time-translational symmetry Wilczek (2012); Khemani et al. (2016); Else et al. (2016); von Keyserlingk et al. (2016); Yao et al. (2017); Choi et al. (2017); Sacha and Zakrzewski (2017); Yao et al. (2018); O’Sullivan et al. (2018); Yao and Nayak (2018); Gambetta et al. (2019). As such, it is suitable for the simulation of a classical spin-1/2 (Ising spin) where the spin states (‘up’ or ‘down’) are represented by the two inequivalent solutions.

Recently, networks of many coupled degenerate optical parametric oscillators have been proposed as suitable platforms to simulate networks of classical spin-1/2 (Ising) systems on a large scale. This kind of simulator was referred to as coherent Ising machine (CIM) Wang et al. (2013), and was proposed as a new platform to efficiently solve complex combinatorial and minimization problems. Because of its potential applications in computation, its recent experimental realization Inagaki et al. (2016); Yamamoto et al. (2017); Böhm et al. (2018) has triggered a significant amount of work, both on the theoretical and computational side King et al. (2018); Hamerly et al. (2018, 2019); Cervera-Lierta (2018); Tiunov et al. (2019); Wang and Roychowdhury (2019).

In light of such potential applications, in this work, we focus on the minimal realization of such network, namely, two coupled degenerate parametric oscillators. Such system has been analyzed in previous studies in the context of MEMS Kenig et al. (2011); Karabalin et al. (2011); Kenig et al. (2012). In the companion letter Bello et al. (2019), we analyzed the system of two coupled parametric oscillators experimentally and theoretically, in view of its application as the building block for a CIM. The two parametric oscillators were experimentally implemented by two radio-frequency cavities, in the presence of a power-splitter non-dissipative coupling. Our main finding was that the system of two coupled oscillators presents a much richer phenomenology than previously analyzed, depending on the values of the system parameters. In addition to a phase where the two oscillators display the expected behaviour of two Ising spins Wang et al. (2013), it was found that, depending on the dissipative or non-dissipative nature of the coupling, the system gives rise to a phase where the dynamics is characterized by limit cycles. In this case, the amplitudes of the two oscillators exhibit periodic beats on top of the fast oscillations at half the frequency of the drive. Such phase lies beyond the dynamics of coupled Ising spins, and its presence may be either a useful resource for CIMs or an additional source of error, which requires further investigation. The goal of this paper is to provide a theoretical background for the experiment in Ref. Bello et al. (2019), as well as a complete characterization of the possible phases in the system of two coupled parametric oscillators.

Specifically, we first consider in Sec. II the linearized equations of motion of the model and solve them using Floquet theory. This approach allows us to determine the stability diagram of the model. Next, in Sec. III, we apply a multiple-scale analysis to determine the phase diagram in the presence of non-linearities. In our study we focus on the universal properties of the model, such as the critical exponents of the different instabilities and the role of the different types of non-linearities (Sec. IV). While most of the article is dedicated to energy preserving couplings, towards the end (Sec. V) we consider a dissipative coupling and connect our findings to the results of Ref. Wang et al. (2013) in the context of CIMs.

II Linear parametric oscillators - Solution by Floquet theorem

We open our analysis by introducing a model of coupled parametric oscillators with purely energy-preserving coupling and showing the explicit solution in the absence of nonlinear terms. In our study, we do not follow the canonical approach to this problem, used for instance in Ref. Landau and Lifshitz (1982) to solve a single parametric oscillator. Here, we instead rely on the Floquet theorem (see e.g. Refs. Magnus and Winkler (1979); Chicone (2006); Eckardt and Anisimovas (2015)), which can be easily generalized to more complicated situations.

II.1 In-phase pumps

In the absence of nonlinearities, a pair of coupled parametric oscillators is described by a set of two generalized linear Mathieu’s equations:

[TABLE]

In Eq. (1), $\omega_{0}$ denotes the proper frequency of the oscillators, $h$ and $\gamma$ represent the intensity and frequency of the pumps, respectively, $g$ is the intrinsic loss term, which we take equal for both oscillators. The coupling $r$ describes an energy preserving coupling between the oscillators: this coupling corresponds to rotations in the $(x_{1},x_{2})$ plane and preserves the total energy, which is proportional to $x^{2}_{1,2}$ for the $x_{1}$ or $x_{2}$ modes, respectively. In the experiment of Ref. Bello et al. (2019) this coupling was implemented by a power splitter coupler. In Eq. (1), the two oscillators are coupled such that the exchange of energy from $x_{1}$ to $x_{2}$, and vice versa, occurs with the same rate, which is determined by $r$. This assumption will be relaxed in Sec. V. In the limit of $r\rightarrow 0$. Eq. (1) becomes equivalent to two decoupled parametric oscillators described by two Mathieu’s equations.

The equations in Eq. (1) can be separated by performing the change of basis $x_{\pm}(t)=x_{1}(t)\pm i\,x_{2}(t)$. In such basis, we have two decoupled parametric oscillators with real and imaginary loss terms:

[TABLE]

The $\dot{x}_{\pm}$ terms in Eq. (2) can be reabsorbed into the definitions of the fields by introducing $x_{\pm}(t)=e^{-(g\pm i\,r)\omega_{0}t/2}\,y_{\pm}(t)$ and then Eq. (2) becomes

[TABLE]

For simplicity, we first focus on the case of $g,r\ll 1$, where one can neglect terms proportional to $(g\pm ir)^{2}$:

[TABLE]

In this limit, $y_{+}(t)$ and $y_{-}(t)$ obey the same equation. In the next subsection we will show how to release this constraint. Since the equations of motion are periodic with period $\mathcal{T}=2\pi/\gamma$, we can look for solutions of the form $y_{\pm}(t)=e^{-i\mu t}\,f(t)$, where $\mu$ is a complex frequency and $f(t)$ is a periodic function with period $\mathcal{T}$. In order to determine $\mu$, we can proceed as follows: using the periodicity of $f(t)$, we can express $y_{\pm}(t)$ in terms of its Fourier components:

[TABLE]

where $A_{n}$ identifies the amplitude of the $n$-th Fourier component, where $n$ is an integer number. If we plug Eq. (5) into Eq. (4) and equate to zero terms multiplying the same oscillating factor, we obtain a recursive equation for the coefficients $A_{n}$:

[TABLE]

where we define $D_{n}(\mu)=\omega^{2}_{0}-{(n\gamma-\mu)}^{2}$. We can more conveniently write Eq. (6) in the matrix form

[TABLE]

The expression in Eq. (7) can be rewritten as $\mathbf{M}_{\mu}\cdot\mathbf{A}=0$, where $\mathbf{A}$ identifies the column vector containing the Fourier components of $f(t)$, and $\mathbf{M}_{\mu}$ is the infinite-by-infinite matrix in Eq. (7). In general, from Eq. (7), the requirement for the existence of a nontrivial solution requires that the determinant of the matrix $\mathbf{M}_{\mu}$ vanishes, ${\rm det}(\mathbf{M}_{\mu})=0$.

We now make the following observation: as we see from Eq. (6), the parametric drive directly couples each Fourier component $A_{n}$ only to its nearest-neighbour ones $A_{n\pm 1}$. This means that $A_{n}$ is coupled to $A_{n\pm m}$, with $m>0$, via a coupling that is of the order of $h^{m}$. Within a perturbative fashion, at first order in $h$, we therefore see that the strongest effect of the parametric drive is coupling $A_{n}$ with $A_{n-1}$ or with $A_{n+1}$. By inspection of Eq. (6), we see that, for $h\rightarrow 0^{+}$, such coupling occurs when $D_{n}(\mu)=D_{n\pm 1}(\mu)=0$. This condition is satisfied for $\mu=\omega_{0}(2n\pm 1)$, which corresponds to the parametric resonance condition $\gamma=2\omega_{0}$.

Since the function in Eq. (5) is periodic with period $\gamma=2\omega_{0}$, we can without loss of generality focus on the case of $n=0$ and therefore consider the situation in which only $A_{0}$ and $A_{-1}$ are coupled by the parametric drive. Therefore, from Eq. (6) and Eq. (7) the requirement for the existence of a nontrivial solution reduces to the condition

[TABLE]

with $\gamma$ sufficiently close to $2\omega_{0}$. We can parametrize the deviation from the parametric resonance condition by introducing a small detuning $\epsilon$ and by rewriting $\gamma=2\omega_{0}+\epsilon$. The polynomial $P_{4}(\mu)$ defined in Eq. (8) has four complex roots. Since we are looking for the parametric resonance between $A_{0}$ and $A_{-1}$, we consider only the degenerate roots that converge to $\mu=-\omega_{0}$ when $h=0$ and $\epsilon=0$, which are found to be

[TABLE]

As evident from Eq. (5), if $\mu_{1,2}$ are complex, then their nonzero imaginary part quantifies the rate of the parametric exponential damping [${\rm Im}(\mu)<0$] or exponential amplification [${\rm Im}(\mu)>0$] of the solution $y_{\pm}(t)$. By expanding $\mu_{1,2}$ for small $h$ and $\epsilon$ and by discarding terms of the order of $\epsilon^{2}h$, the two solutions in Eq. (9) can be written as

[TABLE]

The region such that ${\rm Im}(\mu_{\pm})>0$ (i.e., $\omega_{0}h/2>|\epsilon|$) identifies the linear instability region, where parametric amplification occurs. Within this region, as evident from Eq. (5) and Eq. (10), the oscillator frequency is always exactly equal to $\omega_{0}+\epsilon/2\equiv\gamma/2$, i.e., half of the pump frequency, therefore manifesting its time-crystal nature. Thus, the growing solution of Eq. (4) is found to be $y_{+}(t)\sim e^{(\omega_{0}t/2)\sqrt{(h/2)^{2}-(\epsilon/\omega_{0})^{2}}}\,\cos(\gamma t/2)$, from which it follows that

[TABLE]

Equation (11) describes two solutions $x_{1}(t)={\rm Re}[x_{\pm}(t)]$ and $x_{2}(t)=\pm{\rm Im}[x_{\pm}(t)]$. Within the instability region, these solutions represent two parametrically-driven solutions with in-quadrature beats (in the limit $r\ll 1$). Notice that, from Eq. (11), it is evident that the parametric amplification occurs when the pump strength $h$ is above a threshold value: in this case, parametric amplification occurs when $\sqrt{(h\omega_{0}/2)^{2}-\epsilon^{2}}>\omega_{0}g$. At resonance $\epsilon=0$ (i.e., $\gamma=2\omega_{0}$), the threshold condition reads $h>2g$.

II.2 Varying the pump phase - General derivation

We now use the scheme introduced in Sec. II.1 in order to solve the system of linear Mathieu’s equations in the case when the two oscillators are pumped with a different phase, $\gamma t\to\gamma t\pm\phi/2$:

[TABLE]

where $\phi$ represents the phase difference between the two pumps. Since the oscillators normally lock to the phase of the pump, varying the phase difference between the pumps is equivalent to varying the phase between the two oscillators. We can equivalently rewrite Eq. (12) as

[TABLE]

In the basis $x_{\pm}(t)=x_{1}(t)\pm i\,x_{2}(t)$, Eq. (13) becomes

[TABLE]

By comparing Eq. (14) with Eq. (2), one finds that a finite dephasing $\phi\neq 0$ has two effects: (i) it reduces the strength of the parametric drive from $h$ to $h\cos(\phi/2)$ and (ii) gives birth to an effective coupling between the $x_{+}(t)$ and $x_{-}(t)$ oscillators, whose strength is proportional to $\sin(\phi/2)$. Since now the $x_{\pm}(t)$ oscillators obey different equations of motion, we write the Floquet form of $x_{\pm}(t)$ [see Eq. (5)] in the vector form:

[TABLE]

where $H^{(+)}_{n}$ and $H^{(-)}_{n}$ represent the $n$-th Fourier component of $x_{+}$ and $x_{-}$, respectively. By proceeding as done for Eq. (6), we obtain the recursion relations for $H^{(\pm)}_{n}$:

[TABLE]

where we define $D_{\pm,n}=\omega^{2}_{0}-{(n\gamma-\mu)}^{2}+i\,\omega_{0}(g\pm i\,r)(n\gamma-\mu)$. As done for Eq. (7), Eq. (16) can be also written in the matrix form

[TABLE]

where one defines

[TABLE]

We now define $\Omega_{r}=\omega_{0}\sqrt{1+(r^{2}-g^{2})/4}$ and focus on the limit of small $g$, where

[TABLE]

As in Sec. II.1, we focus on the four-by-four minor of the matrix in Eq. (17) which contains the blocks for $n=0$ and $n=-1$. As for Eq. (8), the requirement for the existence of a nontrivial solution results into the computation of the roots of an (eight-order) polynomial, $P_{8}(\mu)$. Among its eight roots, some of them (which we call the relevant roots) can acquire a nonzero imaginary part depending on which parametric resonance is met. As we will explain below, in this case, there are three distinct resonances, where parametric amplification can occur, which are given by $\gamma=2\Omega_{r},2\Omega_{r}\pm\omega_{0}r$.

II.2.1 Parametric resonance at $\gamma=2\Omega_{r}$

We first discuss the parametric resonance at $\gamma=2\Omega_{r}$. In the limit of $h=0$, one has four relevant roots that are doubly degenerate. They are $\mu^{(0)}_{\pm,\pm}\simeq\pm\omega_{0}r/2-\gamma/2-i\,\omega_{0}g/2$, for which $D_{\pm,0}=D_{\pm,-1}=0$. For finite $h$, the root degeneracy is removed and one finds

[TABLE]

As for the case in Eq. (10), the additional imaginary part $\Gamma^{(0)}(h,\gamma)$, determines the rate of parametric amplification for the resonance around $\gamma=2\Omega_{r}$. The rate $\Gamma^{(0)}$ is identical for $x_{+}$ and $x_{-}$, indicating that both quadratures are amplified by the same amount. Importantly, the oscillation frequencies of the two quadratures differ by $\omega_{0}r$, indicating that the system’s energy oscillates between $x_{1}$ and $x_{2}$. Also, $x_{\pm}(t)$ display fast oscillations at half the frequency of the pump but with in-quadrature beats at frequency $\omega_{0}r/2$ on top of such oscillations. As we will discuss in Sec. III (see also Appendix .3), in the presence of nonlinearities, this behaviour will evolve into a limit cycle that will eventually stabilize the amplitude of the beats.

II.2.2 Parametric resonances at $\gamma=2\Omega_{r}\pm\omega_{0}r$

We now move to the resonances at $\gamma=2\Omega_{r}\pm\omega_{0}r$. In this case, we have only two relevant roots, which are degenerate in the limit $h\rightarrow 0^{+}$. We denote the two roots by $\mu^{(\pm\omega_{0}r)}_{\pm}\simeq-\gamma/2-i\,\omega_{0}g/2$. A nonzero $h$ removes the degeneracy in the two relevant roots, adding an equal and opposite imaginary part, which now we call $\Gamma^{(\pm\omega_{0}r)}(h,\gamma,r)$. Therefore, for finite $h$, one has

[TABLE]

From Eq. (21), we see that both oscillators oscillate with a frequency equal to $\gamma/2$ within all the instability region (for small $h$). Therefore, the system behaves as a time crystal, in which fast oscillations of $x_{\pm}$ are locked in phase with the even or odd cycles of the pump.

The functions $\Gamma^{(0)}(h,\gamma)$ and $\Gamma^{(\pm\omega_{0}r)}(h,\gamma,r)$ in Eqs. (20) and (21) identify the three different instability regions around the three parametric resonances $\gamma=2\Omega_{r},2\Omega_{r}\pm\omega_{0}r$, respectively. In contrast to the case discussed in Sec. II.1, we could not find an analytical expression for $\Gamma^{(0)}(h,\gamma)$ and $\Gamma^{(\pm\omega_{0}r)}(h,\gamma,r)$. Therefore, in order to determine the stability phase diagram, we resort to numerics.

II.3 Linear instability regions

The regions of linear instability can be numerically determined from the imaginary part of the roots of the polynomial $P_{8}(\mu)$, for different values of $h$ and $\gamma$. This determines {\color[rgb]{0,0,0}-\omega_{0}g/2}+\Gamma^{(0)}(h,\gamma) or {\color[rgb]{0,0,0}-\omega_{0}g/2}+\Gamma^{(\pm\omega_{0}r)}(h,\gamma,r) in Eqs. (20) and (21). An example of the instability regions in the $h$ vs. $\gamma-2\Omega_{r}$ plane is shown in Fig. 1. For concreteness, we show the instability regions for small $g$, which we choose $g=10^{-2}$, and $r=0.2$, and for (a) $\phi=0$, (b) $\phi=\pi/2$ and (c) $\phi=\pi$. For $\phi=0$, the instability region consists of one cone centred at $\gamma=2\Omega_{r}$. For nonzero $\phi$, two additional outer instability regions appear centred at $\gamma-2\Omega_{r}=\pm\omega_{0}r=\pm 0.2\,\omega_{0}$ around the central region. For $\phi=\pi$, the resonance at $\gamma=2\Omega_{r}$ is completely suppressed, and only the instability regions around $\gamma-2\Omega_{r}=\pm\omega_{0}r$ are found.

A better insight regarding the properties of the system is given by studying the behaviour of the real and imaginary parts of the relevant roots as we vary $\gamma$ and $h$. This is shown in Fig. 2, focusing in particular on the instability phase diagram computed at $\phi=\pi/2$ [panel (b) of Fig. 1], which is also reported in Fig. 2 for completeness. The blue dots represent the four relevant roots [Eqs. (20) and (21)], plotted by showing their imaginary parts as a function of their real parts from which we subtract $\gamma/2$. We do not show the other four roots since they do not contribute to the instabilities of the system and therefore are not relevant for the present discussion.

We show the relevant roots in five prototype cases in Fig. 2: (a) inside the stable region [${\rm Im}(\mu)<0$, i.e., the white area (S) of the insets], all roots have negative imaginary part, which is equal to {\color[rgb]{0,0,0}-\omega_{0}g/2} in the small $g$ and $r$ limit. This case corresponds to an exponential damping in time for both oscillators. Inside the instability regions, one or two pairs of roots acquire in addition an equal and opposite imaginary part, depending on which one of the instability regions is entered. In the case of the outer instability regions [panel (c)], only two roots acquire an additional factor $\Gamma^{(\pm\omega_{0}r)}$ while having a real part that is locked to $\gamma/2$ [see also Eq. (21)], which is highlighted in the figure by the cyan dashed vertical line. Instead, inside the central instability region [panel (e)], four roots acquire an additional imaginary part $\Gamma^{(0)}$ while having a real part locked to $\gamma/2\pm\omega_{0}r/2$ [see also Eq. (20)], which is instead highlighted by the red vertical lines. These two types of instabilities correspond, respectively, to a Pitchfork bifurcation and to a Hopf instability Strogatz (2007). Parametric amplification occurs when the unstable region [${\rm Im}(\mu)>0$, i.e., the yellow area (U) in the insets] is entered, i.e., when the overall imaginary part is such that either {\color[rgb]{0,0,0}-\omega_{0}g/2}+\Gamma^{(\pm\omega_{0}r)}>0 for the outer regions, or {\color[rgb]{0,0,0}-\omega_{0}g/2}+\Gamma^{(0)}>0 for the central region, which identifies the threshold for parametric amplification.

The key result of this analysis is that, in our system, when the parametric linear instability is met, the two modes are amplified and oscillate with a frequency that is locked to half of the frequency of the pump. The parametric amplification can occur with or without beats, depending on which region of linear instability (the central or the outer ones, respectively) is entered. Such linearly unstable regions are the precursors of the stable regions of limit cycle and synchronized oscillations, in which nonlinear effects eventually stabilize the long-time dynamics. This will be the topic of the next sections (see also Appendix .3).

As a final remark, we stress that the advantage of using the perturbative method here presented is that it grants us a good analytical control. Quantitatively, the result presented in this section are valid strictly speaking in the limit of $h\rightarrow 0^{+}$. For a not too large finite value of $h$, there will be corrections to our findings, but the qualitative picture remains valid. For the sake of completeness, we mention that the full numerical solution can be obtained by resorting to the formalism of fundamental matrices Chicone (2006).

III Nonlinear case - Perturbative multiple-scale analysis

The method based on Floquet’s theorem presented in Sec. II allows us to systematically study systems of linear coupled parametric oscillators, but it cannot be applied in the presence of nonlinearities. When a nonlinear term is included in the equations of motion, one can resort to a multiple-scale perturbation method Kevorkian and Cole (1996) in order to determine the long-time dynamics of the oscillators. The goal of this section is to apply such method in order to study the dynamics of the system of coupled oscillators [Eq. (12)] in the specific case where a quadratic nonlinearity is included in the model.

III.1 Equations for the long-time dynamics

In the actual physical context, there are different sources of nonlinearities that can appear in the equations of motion [Eq. (12)], such as saturation, Kerr effects or pump depletion. In this section, we focus on one type of nonlinearity, namely, the pump depletion, which is in many experimental contexts the most relevant type of nonlinearity. We postpone the discussion of other types of nonlinearities to Sec. IV.

The pump depletion accounts for the fact that the pump intensity is depleted within the nonlinear medium, by means of down-conversion processes to the signal (and idler) field. In the limit of small depletion, we can therefore write the equations of motion as

[TABLE]

Here, $\beta$ quantifies the pump depletion. We now focus on the resonant case $\gamma=2\Omega_{r}$ where the system is more affected by the parametric instability. In order to determine the non-trivial long-time dynamics of the system, we proceed with a multiple-scale perturbative expansion Kevorkian and Cole (1996). The details of the calculation are reported in Appendix .1 for the sake of completeness.

In Eq. (22), we identify $\omega_{0}$ as the largest frequency scale, which identifies the fastest time scale of the system $t=2\pi/\omega_{0}$. We assume that the coupling constants, $r$, $g$ and $h$ are much smaller than unity, and influence the dynamics of $x_{1}(t)$ and $x_{2}(t)$ only on time scales which are much longer than $2\pi/\omega_{0}$. In these conditions, the full dynamics can be separated in fast-varying and slow-varying degrees of freedom. If we work at fixed $g$, which we take as the small expansion parameter of the theory, we can identify the characteristic time scale of the slow-varying degrees of freedom by $\tau=gt$. Therefore, we can write $x_{1}(t,\tau)=A(\tau)\,e^{i\omega_{0}t}+A^{*}(\tau)\,e^{-i\omega_{0}t}$ and $x_{2}(t,\tau)=B(\tau)\,e^{i\omega_{0}t}+B^{*}(\tau)\,e^{-i\omega_{0}t}$, where $e^{\pm i\omega_{0}t}$ describes fast oscillations at frequency $\omega_{0}=\gamma/2$, and $A(\tau)$ and $B(\tau)$ represent the slow-varying complex amplitudes for $x_{1}$ and $x_{2}$, respectively. In the following, we express time in units of $\omega_{0}$, i.e., we define $\tilde{\tau}=\omega_{0}\tau$, and it is convenient to redefine $h$ and $r$ with respect to $g$, i.e., by introducing $\tilde{h}=h/g$ and $\tilde{r}=r/g$.

By separating the real and imaginary parts of each complex amplitude, i.e., $A=A_{R}+i\,A_{I}$ and $B=B_{R}+i\,B_{I}$, the main result is that the dynamics of the slow-varying amplitudes is described by a set of four coupled equations:

[TABLE]

The system in Eq. (23) encodes the dynamics of $A$ and $B$. According to the standard analysis of nonlinear systems Strogatz (2007), all the informations that we need in order to describe the long-time dynamics of $x_{1,2}(t)$ can be found by studying the configuration of the fixed points of Eq. (23), which are found by imposing $\partial A/\partial\tilde{\tau}=\partial B/\partial\tilde{\tau}=0$. Their stability is determined by the eigenvalues of the Jacobian matrix at a specific point (see Appendix .5). Because we were not able to find an analytic expression for the fixed points, we resorted only to the numerical solution of Eq. (23). For the configuration of the fixed points in the decoupled case ($\tilde{r}=0$) the reader is referred to Appendix .2.1

III.2 Phase diagram of two coupled oscillators

The key result for $\phi=0$ is reported in Fig. 3, in which we show the phase diagram in the $h/(2g)$ vs. $r/g$ plane. According to our analysis, three main regions are found: when the system is below the threshold for parametric amplification [region (I)], the origin $A=B=0$ is the only stable attractor, and each trajectory of $A(\tau)$ and $B(\tau)$ is attracted into the origin. In this region, oscillations are suppressed in the long-time limit.

For any finite $\tilde{r}$, as the pump intensity is increased up to a threshold value identified by the red dash-dotted line in Fig. 3, the origin becomes a saddle point giving birth to a stable limit cycle in its surrounding via a supercritical Hopf bifurcation [region (II)]. In this region, the two oscillators display everlasting beats, whose frequency, close to the threshold, is determined by $\omega_{0}r$ and whose shape changes as $\tilde{h}$ is increased. For the analytical derivation of the boundary between region (I) and region (II) the reader is referred to Appendix .3.

As the pump intensity is further increased, stable attractor and saddle nodes are born in pairs via saddle-node bifurcations. Specifically, depending on the value of $\tilde{r}$, by increasing $\tilde{h}$, a region with either four [region (III) ] or eight [region (IV)] stable fixed points is entered. Inside these regions, the limit cycle disappears and the amplitude of the oscillations becomes constant in time (synchronized). The transition line between the two regimes has been numerically determined and is identified by the blue dashed line in the figure. We find that such a transition can occur in two different ways. For small couplings, we find that the period of the limit cycle, which represents the period of the beats on top of the fast oscillations at half the pump frequency, diverges as region (III) is approached (see also Sec. III.3). At the boundary between region (II) and region (III), eight fixed points (four attractors and four saddle points) are born on the limit cycle via a saddle-node bifurcation, causing the extinction of the limit cycle as region (III) is entered. This phenomenology is customary referred to as an infinite-period bifurcation. For larger values of the coupling, we find that there is a first area inside region (IV), in which the limit cycle can coexist with the stable attractors and therefore fast oscillations occur either displaying beats or with constant amplitude, depending on the initial conditions. After this region, the limit cycle collapses into one of the attractors, and therefore only synchronized oscillations are found. For even larger values of $\tilde{h}$, a region with sixteen stable fixed points [region (V)] is found.

As $\phi$ is increased from [math] to $\pi$, the region (II) in which the limit cycle is found tends to become smaller and smaller, and eventually completely disappears for $\phi=\pi$. Figure 4 shows the two boundaries for the supercritical Hopf bifurcation from region (I) to region (II), and for the infinite-period bifurcation from region (II) to region (III) discussed in Fig. 3, for different values of $\phi$. For $\phi=\pi$, the system directly passes from the below-threshold region to the synchronization one. Importantly, for the range of parameters considered here, the synchronization region after the limit cycle region, for small values of the coupling, is always found with four stable fixed points for all values of $\phi$.

Before concluding this section, we comment on the physicality of the model. As shown in Fig. 3, our model predicts a large number of fixed points. Because we are describing a system of two parametric oscillators, one would expect that, in the synchronized regime, the spin picture holds when the system has only four stable fixed points on the real axis (twice as many as a single parametric oscillator). Indeed, for experimental purposes, the model that we consider is relevant only for values of the pump that are not too far away from the oscillations threshold (i.e., when the single oscillator has two stable fixed points, see Appendix .2.1). This is precisely the situation discussed in Ref. Bello et al. (2019). In order for the model to be physical, it is important to verify that the condition $1-\beta({\bar{x}_{1}}^{2}+{\bar{x}_{2}}^{2})/2>0$ holds inside the full phase diagram, where $\bar{x}^{2}_{1,2}$ is the long-time average of $x^{2}_{1,2}(t)$. This condition ensures that, on average, the energy of the pump is always down converted to the optical fields. We have verified that such condition holds inside the numerically explored phase diagram, even for very large values of the pump intensity.

Notice that the experimentally accessed region is usually up to $h/(2g)\simeq 2$, which corresponds to the maximum of the conversion efficiency (i.e., unity) Schiller et al. (1999); Martinelli et al. (2001); Sturman and Breunig (2011); Breunig et al. (2011); Breunig (2016). To the best of our knowledge, the regions that are found for larger pump intensities, in which the model displays a number of stable points larger than four, remain experimentally unexplored. Our theory indicates the existence of interesting dynamics also in this high-pump intensity range. A deeper analysis of such regions is beyond the aim of the present manuscript and remains a subject of future studies.

III.3 Critical scaling by a three-scale analysis

In this section, we determine the critical exponent for the radius of the limit cycle close to the supercritical Hopf bifurcation boundary between regions (I) and (II), and of the period of the limit cycle close to the infinite-period bifurcation between region (II) and (III) of the phase diagram in Fig. 3.

In order to determine the scaling of the radius of the limit cycle, we focus on the points in region (II) which are close to the threshold $\tilde{h}_{\rm th}=2$ and sufficiently far from the infinite-period bifurcation. In this case, the condition $h\ll r\ll 1$ holds and the dynamics of the oscillators is determined by three different characteristic frequencies: $\omega_{0}$ (fast oscillations), $\omega_{0}r$ (medium-scale beats) and $\omega_{0}h$ (long-time overall amplitude). In this condition, it is natural to perform a multiple-scale analysis by introducing three different time scales. If we redefine $r=r_{0}\,\tilde{r}$ and $h=g\,\tilde{h}$, where $r_{0}$ is a characteristic scale for the dynamics of the beats, we can distinguish three different time scales in the expansion: $t$ (for fast oscillations), $\sigma=r_{0}t$ (for the medium-scale dynamics, i.e., beats) and $\tau=gt$ (for the slow dynamics).

Focusing on the resonant case $\gamma=2\omega_{0}$, it is convenient to rewrite the equations of motion in Eq. (22) using the $x_{\pm}=x_{1}\pm i\,x_{2}$ basis, as in Sec. II:

[TABLE]

We now expand $x_{\pm}=x^{(0)}_{\pm}+r_{0}\,x^{(B)}_{\pm}+g\,x^{(1)}_{\pm}$. By proceeding as in Sec. III.1 (see also Appendix .1), one obtains the equations for the fast-varying and medium-scale modes:

[TABLE]

and similarly the equations for the slow-varying modes

[TABLE]

From Eqs. (25) and (26), since $x_{\pm}=x_{\mp}^{*}$ and since the equations for the slow-varying modes of $x_{+}$ and $x_{-}$ are mutually complex conjugated, we can use solutions of the form

[TABLE]

where $C_{\pm}=C_{\mp}^{*}$ describe the medium-scale modes, and the slow-varying modes are described by the same complex amplitude $A_{S}$. By plugging these expressions into Eq. (25b), one has the solvability condition for the medium-scale dynamics:

[TABLE]

from which one obtains the beating factor $C_{\pm}(\sigma)=C(0)e^{\mp i\omega_{0}\tilde{r}\sigma/2}=C(0)e^{\mp i\omega_{0}rt/2}$, where $C(0)=|C(0)|e^{iu}$ is a complex number. Here, $|C(0)|$ is a normalization factor and $u$ determines the initial phase of the beats. By using these expression for $x_{\pm}^{(0)}$ in Eq. (26), and by neglecting oscillating factors as $e^{\pm i2\omega_{0}rt}$ that are strongly oscillating on the slow time scale, one obtains the solvability condition for the slow-varying amplitude $A_{S}$:

[TABLE]

which is nothing but Eq. (23a) without the $B$ term, and with the replacement $\beta\rightarrow 3\beta/4$. For its solution, the reader is referred to Appendix .2.1. For completeness, we recall below the main results.

Basing on the notation used in Appendix .2.1, close to the threshold value $\tilde{h}\simeq 2$, the origin is a saddle point and other four fixed points (two saddle points and two stable nodes) are found in its surroundings. Such points appear on the imaginary and real axes, respectively. We call the stable fixed point $W^{\prime}_{+}$, whose polar coordinates in the ${\rm Im}(A_{S})$ vs. ${\rm Re}(A_{S})$ plane are [Eq. (A7)] $\varphi_{W^{\prime}_{+}}=0$ and $R_{W^{\prime}_{+}}=\sqrt{(2-4/\tilde{h})/(3\beta)}$. Such a point is the only one that can stabilize the long-time dynamics of $x_{\pm}$. In particular, its coordinates determine the amplitude of the beats. Therefore, by recalling that, for $\phi=0$, $\tilde{h}_{\rm th}=2$, the radius of the limit cycle from the long-time dynamics of $A_{S}$, which we call $R_{\rm LC}\coloneqq\lim_{t\rightarrow\infty}|A_{S}(gt)|$, is readily determined:

[TABLE]

At the onset of the supercritical Hopf bifurcation, the limit cycle grows from zero amplitude with the critical exponent $1/2$. In terms of the $A$ and $B$ amplitudes determined by Eq. (23), the limit cycle is therefore identified by the dynamics

[TABLE]

In the long-time limit, the limit cycle is therefore a perfect circular arc whose frequency is determined by $r$ and whose radius solely depends on $(\tilde{h}-\tilde{h}_{\rm th})/\tilde{h}$ [Eq. (30)].

In order to explicitly show the critical exponent $1/2$ of the supercritical Hopf bifurcation, we compute the radius of the limit cycle as a function of $\tilde{h}$ for a fixed $\tilde{r}\neq 0$ by numerically solving Eq. (23) for $\phi=0$ and $\beta=10^{-2}$ close to the boundary of the supercritical Hopf bifurcation $\tilde{h}_{\rm th}=2$ (red dash-dotted line in Fig. 3). The result is shown in Fig. 5. We then superimpose the numerically determined data with the expected behaviour [Eq. (30)] found by the three-scale analysis. The agreement between the two behaviours confirms the prediction found in Eq. (30). Such square-root scaling can be further highlighted by rescaling the data and the analytical behaviour by introducing $x=\ln(1-\tilde{h}_{\rm th}/\tilde{h})$ and, from Eq. (30), $f(x)=\ln\sqrt{2/(3\beta)}+x/2$, which is shown in the inset.

The same critical exponent is found by studying the behaviour of the period of the limit cycle $T_{\rm LC}$ close to the infinite-period bifurcation Strogatz (2007) (blue dash line in the phase diagram in Fig. 3). The result of the simulation is shown in Fig. 6. We show the numerically determined value of $T_{\rm LC}$ as a function of the distance from the critical line $\tilde{h}-\tilde{h}_{c}$ at a fixed $\tilde{r}=5\times 10^{-3}$. The continuous line represents the best fit of the form $T_{\rm LC,fit}=c_{1}+c_{2}/\sqrt{\tilde{h}_{c}-\tilde{h}}$, where $c_{1}$ and $c_{2}$ are fit parameters. As evident from the figure, the agreement between the numerical data and the fit confirms the fact that the period of the limit cycle diverges as $T_{\rm LC}\sim{(\tilde{h}_{c}-\tilde{h})}^{-1/2}$ as the infinite-period bifurcation is approached.

Before concluding this section, we mention that the whole analysis remains valid if a different form of coupling is considered, i.e., $\omega_{0}^{2}r\,x_{2,1}$ in the equation of motion (22) of $x_{1}$ and $x_{2}$, respectively Wang et al. (2013). The proof of this statement is discussed in Appendix .4.

IV Different types of nonlinearity

In this section, we comment on the effects of different nonlinearity of the two oscillators. We show that the phenomenology discussed in Sec. III is not a consequence of the specific choice of the model, but it is common to other models which can be relevant in different experimental contexts. The properties that we discuss in this section are found by using exactly the same tools discussed in the previous sections. We therefore report the main results without explicitly showing all the technical details.

In order to ease the notation, we rewrite the equations of motion in a more compact and generic form as

[TABLE]

in which $F_{\rm NL;1,2}$ identify the nonlinear terms. In Eq. (22), we considered the pump-depletion nonlinearity $F_{\rm NL,1}=-h\beta\omega^{2}_{0}\sin(2\omega_{0}t)\,x^{3}_{1}$, $F_{\rm NL,2}=-h\beta\omega^{2}_{0}\sin(2\omega_{0}t+\phi)\,x^{3}_{2}$. As mentioned before, another possible nonlinearity arises from a Kerr or saturation effect. In this case, the cubic term in Eq. (22) will not be coupled to the pump, and the equations of motion in Eq. (32) are now written with $F_{\rm NL;1,2}=-\beta\omega^{2}_{0}\,x^{3}_{1,2}$ (see Table 1), whose multiple-scale equations are obtained as done for Eq. (23).

The phase diagram that we obtain from the configurations of the fixed points, which is shown in Fig. 7, panels (b), displays the same phases found for the phase diagram in Fig. 3, which is reported for completeness in panels (a): (I) a region in which the system is below threshold, (II) an extended region in which only a stable limit cycle can stabilize the dynamics, therefore yielding also in this case everlasting beats in the time evolution of $x_{1}$ and $x_{2}$, and (III a)-(III) a region in which stable attractors stabilize the dynamics. However, differently from the model analyzed in Fig. 3, in addition to region (III) in which four attractors are found, for $\phi>0$ [Fig. 7, panel (b2)] there is an additional intermediate region, which we call region (III a), in which only two attractors are found.

As discussed in Fig. 4, for $\phi=\pi$ [panels (a3) and (b3)], the region with the limit cycle [region (II)] disappears, and one passes directly from the below-threshold region to the region with stable attractors, which are four in the case of the pump-depletion nonlinearity and two in the case of the Kerr nonlinearity. In the latter case, for larger values of the pump, the region with four stable attractors is found above the one with two attractors only. In a more physical situation in which both nonlinearities are found, one always finds, for $\phi>0$, a small region with two stable fixed points before the one with four stable points, as the pump intensity is increased. Interestingly, when the Kerr nonlinearity is considered, the behaviour of the the radius of the limit cycle at the onset of the supercritical Hopf bifurcation (for $\phi=0$) is found to grow from zero with a critical exponent equal to $1/4$ in contrast to the critical exponent $1/2$ found in the case of the pump-depletion nonlinearity [see also Eq. (A13)]. Such difference can be exploited in experiments to distinguish between the two types nonlinearities. This point is left for future work.

From this analysis, apart from the specific quantitative details that depend on the specific model that we consider, it is therefore seen that the presence of a wide region in which the system displays everlasting beats comes solely from the interplay between parametric gain, losses, nonlinearity, coupling and, apart from extremely fine-tuned phase differences between the two oscillators, it always emerges for any nonzero coupling as the oscillation threshold is crossed (see Sec. V for the extension to the case of dissipative coupling). When the system is within such phase, the system of two coupled parametric oscillator therefore does not match the picture of two Ising spins discussed in previous work Wang et al. (2013); Inagaki et al. (2016); Yamamoto et al. (2017); Hamerly et al. (2019).

V Dissipative coupling and CIM

So far, we considered only the effect of an energy-preserving coupling, which is found when the energy exchange rate between the two oscillator is balanced. In this section, we consider the effects of a dissipative coupling, which is instead found whenever the two oscillators exchange energy with different rates, in order to connect our model to the one discussed in Ref. Wang et al. (2013) in the context of CIMs. We eventually a generic experimental implementation of such a coupling.

In order to show this connection, we first consider the linear case. We rewrite Eq. (12) as

[TABLE]

where all quantities are as in Eq. (33), and $\alpha\geq 0$ represents the strength of the dissipative part of the coupling that quantifies the unbalancing between the energy exchange rates between the two oscillators. As we will show, depending on the relation between $r$ and $\alpha$, the system undergoes a transition between the CIM behaviour Wang et al. (2013) and the beating phenomenology discussed in the previous sections.

In order to diagonalize Eq. (33), we introduce the basis ${(k_{-}\,\,k_{+})}^{T}=K(x_{1}\,\,x_{2})^{T}$, where $T$ denotes the transposition and the non-unitary matrix $K$ is

[TABLE]

In this basis, the equation of motion (33) becomes

[TABLE]

From Eq. (35) and using the same tools as in Sec. II, one can see that two different regimes arise:

For $r>\alpha$, i.e., when the nature of the coupling is mostly non dissipative, the term $\sqrt{r^{2}-\alpha^{2}}$ is real. This is the situation studied in Sec. II. When $\gamma=2\omega_{0}$, the solution displays beats at a frequency $\pm\omega_{0}\sqrt{r^{2}-\alpha^{2}}/2$. Two additional parametric resonances at $\gamma=2\omega_{0}\pm\omega_{0}\sqrt{r^{2}-\alpha^{2}}$ are found, at which parametric amplification occurs without beats; 2. 2.

For $r<\alpha$, i.e., when the dissipative part of the coupling dominates, the term $\sqrt{r^{2}-\alpha^{2}}=i\,\sqrt{\alpha^{2}-r^{2}}$ is imaginary. Now, only the parametric resonance at $\gamma=2\omega_{0}$ is found, for all values of $\phi$, and the solution never displays beats. Instead, the modes $k_{\pm}$ have now different loss terms $g\mp\sqrt{\alpha^{2}-r^{2}}$ for $k_{\pm}$, leading to different oscillation thresholds.

The systems therefore undergoes a transition between the CIM to the beating behaviour at $r=\alpha$. The analysis can be extended to the non-linear case by including the pump-depletion nonlinearity, and proceeding with the two-scale expansion as in Eq. (23), see Appendix .3 for more details. By using this method, we compute the phase diagram in the $h/(2g)$ vs. $r/g$ plane, as done in Figs. 3 and 7, for different values of $\alpha$ (which we also rescale as $\tilde{\alpha}=\alpha/g$) and $\phi$.

The result is shown in Fig. 8. The phase diagram for $\tilde{\alpha}=0$ and $\phi=0$ is the same as in Figs. 3 and 7, panel (a1), and is reported here for completeness. As explained there, when crossing the oscillation threshold $h=2g$, one always enters the beating region. For stronger pumps, the system undergoes a transition to a region with four (or eight) fixed points. The picture changes when $\tilde{\alpha}>0$. In particular, we choose $\tilde{\alpha}=0.15$ and first show the result for $\phi=0$. We see that, in this case, an additional phase with two stable fixed points emerges, and for $\tilde{r}<\tilde{\alpha}$, this phase is found directly above the threshold (see Appendix .3 for the analytical computation). For larger values of $\tilde{h}$, from the region with two stable points, the phase with four stable fixed points is found. This situation matches the one discussed in Ref. Wang et al. (2013), in which the two-oscillator system can be used as a CIM directly above the oscillation threshold. For $\tilde{r}>\tilde{\alpha}$, the limit cycle region discussed throughout this manuscript emerges between the below-threshold and the CIM regions. This analysis suggests that there are two different routes to reach the CIM regime, whose further analysis is left for future work.

For $\phi>0$, the picture remains qualitatively similar. The width of the region with two stable fixed points, as well as the limit cycle region, is reduced as $\phi$ is increased, as discussed in the previous sections for $\tilde{\alpha}=0$. For $\phi=\pi$, only the region with four stable fixed points is found above threshold.

Before concluding, we discuss a generic experimental setup to realize our system as in Eq. (33). A minimal setup is reported in Fig. 9, which can be implemented both by means of radio-frequency Bello et al. (2019) or optical components, thus ensuring the scalability of the setup. The two fields $x_{1}$ and $x_{2}$ are generated inside two cavities $A$ and $B$, respectively, and they are coupled by a power-splitter coupling. In the most general case, this component accounts for the following quantities: (i) the transmittance coefficients $c_{AA}$ and $c_{BB}$, which can be without loss of generality taken equal for both cavities, i.e., $c_{AA}=c_{BB}=c\geq 0$, whose effect is to renormalize the intrinsic loss of the cavities $g$, and (ii) the coupling coefficients $c_{AB}$ and $c_{BA}$ that when they have the same sign, determine the amount of energy that is transferred from $A$ to $B$, and from $B$ to $A$, respectively. Without loss of generality, we consider $c_{BA}>0$.

In this notation, the coupling between the two oscillators can be written as

[TABLE]

When $c_{AB}=c_{BA}$, one defines $r=c_{AB}=c_{BA}$, and this balanced coupling leads only to the presence of beats, without any CIM region. Instead, when the coupling is unbalanced, i.e. $c_{AB}\neq c_{BA}$, it is possible to achieve the CIM regime. In this case, one can write $r=(c_{BA}+c_{AB})/2$ and $\alpha=(c_{BA}-c_{AB})/2$, so that $c_{AB}=r-\alpha$ and $c_{BA}=r+\alpha$, and Eq. (36) becomes

[TABLE]

By taking the time derivative on both sides of Eq. (37), and by including this coupling in the equations of motion, Eq. (33) is obtained.

VI Conclusions

In this work, we reported a detailed analytical and numerical analysis of two parametric oscillator coupled by a power-splitting coupling, first focusing on the case in which the coupling was purely energy-preserving, and later discussing the relation of our model with CIMs in the case of a dissipative coupling.

We first studied in detail the linear case by resorting to the Floquet theorem. We analytically showed that the system displays three resonances, whose relative splitting in frequency depends on the coupling strength, and then we numerically determined the full stability phase diagram. We showed that, depending on what resonance is met, parametric amplification for both oscillators can occur with or without the beats. In the former case, the frequency of the beats is solely determined by the coupling strength

We then discussed the nonlinear case, first by studying in detail the model with one specific type of nonlinearity, namely, the pump depletion. Next, we corroborated the generality of our finings by discussing the validity of our results in different models, considering different types of nonlinearity and coupling. A single parametric oscillator, above the oscillation threshold, has two possible solutions that are identified by a relative time shift of $\pi$ (one period of the pump). For this reason, a single parametric oscillation is suitable for the simulation of a classical spin-1/2 degrees of freedom, the two states of the spin being identified by the two solutions. In contrast, we showed that two nonlinear coupled parametric oscillators display a wide region in parameter space in which, sufficiently not too far away from the oscillation threshold, only a stable limit cycle is found and oscillations occurs with everlasting beats whose shape and frequency depends on the system parameters. This phenomenology was found as long as the nature of the coupling was mostly non-dissipative, irrespective of the details of the nonlinearity and away from extremely fine-tuned values of the phase difference between the pumps.

Our findings, from a generic perspective, show a way to use parametric oscillators in order to preserve coherence indefinitely. On the other hand, they are immediately relevant to the context of CIMs. Indeed, given the richer physics that we found in the minimal building block of the two-oscillator system with respect to what has been previously addressed, it is crucial to understand how the interplay between two couplings of different nature affects a more structured network. For this reason, the extension of the study presented in this manuscript to more than two coupled parametric oscillators, specifically, studying the fate of the limit cycle when several oscillations are coupled is an important step in the analysis of large-scale CIMs. We leave this point as an outstanding perspective for future work.

**ACKNOWLEDGEMENTS **

We thank Joseph Avron, Ivan Bonamassa, Claudio Conti, Nir Davidson, Igor Gershenzon, Ron Lifshitz, Chene Tradonsky, and Yoshihisa Yamamoto for fruitful discussions. We are grateful to David A. Kessler for careful reading and invaluable comments on this manuscript. A. P. acknowledges support from ISF grant No. 46/14. M. C. S. acknowledges support from the ISF grants No. 231/14 and 1452/14.

Appendix

.1 Details on the derivation of Eq. (23)

In this appendix, we report the details of the derivation of the system in Eq. (23). As we discussed in the main text, one first separates the fast-varying time scale $t$ from the slow varying one $\tau=gt$. We now proceed with the perturbative expansion, treating $g$ as the small expansion parameter, and consider only terms in the expansion that are at most of the order of $g$. First, in the dynamics of $x_{1}$ and $x_{2}$, we can explicitly separate the fast-varying time scale from the slow-varying one, i.e., we write $x_{1,2}=x_{1,2}(t,\tau)$. We can therefore express the time derivative as $d/dt=\partial/\partial t+g(\partial/\partial\tau)$, and therefore $d^{2}/dt^{2}\simeq\partial^{2}/\partial\tau^{2}+2g(\partial/\partial\tau)(\partial/\partial t)$, where we neglect terms of the order of $g^{2}$. Similarly, we expand $x_{1,2}=x^{(0)}_{1,2}+g\,x^{(1)}_{1,2}$, where $x^{(0)}_{1,2}$ and $x^{(1)}_{1,2}$ represent the zero-order and first-order correction to $x_{1,2}$, respectively.

Using these definitions into Eq. (22) and the fact that $2\Omega_{r}\simeq 2\omega_{0}$, we can separate the terms that do not appear multiplied by $g$, which are

[TABLE]

from the terms that are proportional to $g$, which are

[TABLE]

From Eq. (A1), we can write $x_{1}(t,\tau)=A(\tau)\,e^{i\omega_{0}t}+A^{*}(\tau)\,e^{-i\omega_{0}t}$ and $x_{2}(t,\tau)=B(\tau)\,e^{i\omega_{0}t}+B^{*}(\tau)\,e^{-i\omega_{0}t}$, where $A(\tau)$ and $B(\tau)$ represent the slow-varying complex amplitudes for $x_{1}$ and $x_{2}$, respectively. If these expressions of $x^{(0)}_{1,2}$ are used into Eq. (A2a) and Eq. (A2b), one has [we show explicitly the calculation for Eq. (A2a) only, the one for Eq. (A2b) being essentially the same]

[TABLE]

where ${\rm c.c.}$ denotes the complex conjugation. The terms proportional to $e^{\pm i\omega_{0}t}$ in Eq. (A3), which are commonly referred to as secular terms, represent a resonant driving force applied to the $x^{(1)}_{1}$ oscillator. Such a force, will always cause the solution for $x^{(1)}_{1}$ to be unbounded. In order to ensure the solvability of Eq. (A3), we need to impose that such secular terms are zero. This gives the solvability condition for Eq. (A3):

[TABLE]

By separating real and imaginary part of $A$ and $B$, i.e., $A=A_{R}+i\,A_{I}$ and $B=B_{R}+i\,B_{I}$, we can write the two coupled equations for $A_{R}$ and $A_{I}$ shown in Eq. (23a). By repeating the same steps for Eq. (A2b), we therefore arrive to the set of four coupled equations for the real and imaginary parts of the complex amplitudes of the fields in Eq. (23).

.2 Stability analysis of the single parametric oscillator

In this appendix, we report for completeness the stability analysis of the nonlinear Mathieu’s equation for the single parametric oscillator in the presence of the pump-depletion or Kerr nonlinearity.

.2.1 Pump-depletion nonlinearity

We first focus on the case of the pump-depletion nonlinearity, i.e., the case discussed in Sec. III.1, for $\tilde{r}=0$. In this case, the two oscillators are decoupled and one can study only the dynamics of one of the two oscillators $A$ or $B$ in Eq. (23), since for $r=0$ the effect of $\phi$ is trivial and the two equations of motion describe exactly the same physics. In order to simplify the analytical calculation, we therefore study the equations for $A$ in Eq. (23), for $\phi=0$, which for completeness we recall below ($\tilde{\tau}=\omega_{0}\tau$):

[TABLE]

It is convenient to find the coordinates of the fixed points in the $A_{I}$ vs. $A_{R}$ plane in polar coordinates. We can therefore define $A_{R}=R\,\cos(\varphi)$ and $A_{I}=R\,\sin(\varphi)$. From Eq. (A5), the condition $\partial A/\partial\tilde{\tau}=0$ therefore yields the set of equations for $R\neq 0$:

[TABLE]

From Eqs. (A6), one sees that $\varphi=0$ and $\varphi=\pi/2$ are two possible solutions for the angular variable, which define two sets of fixed points that we call $W_{\pm}$, i.e., $\varphi_{W_{-}}=\pi/2$ and $\varphi_{W_{+}}=0$. The corresponding radial variables are found to be

[TABLE]

Instead, for $\varphi\neq 0,\pi/2$, one has from Eq. (A6a)

[TABLE]

and if this is used in Eq. (A6b), one has the solution for the angular variable

[TABLE]

for $\tilde{h}>4$. Equation (A9) identifies two additional fixed points that we call $P$ and $Q$, whose angular and radial variables are therefore

[TABLE]

We therefore have the following picture: for $\beta>0$, $W_{-}$ is found for any $\tilde{h}>0$, $W_{+}$ is found for $\tilde{h}>2$, whereas $P$ and $Q$ are found for $\tilde{h}>4$. From the expression of $\mathbf{J}_{1}$ in Appendix .5, the eigenvalues of the Jacobian are found to be

[TABLE]

for a given fixed point, i.e., for the specific values of $R$ and $\varphi$. The eigenvalues of the Jacobian matrix are independent of $\beta$, since $R\sim 1/\sqrt{\beta}$ for all fixed points with $R\neq 0$.

There are four main different situations depending on the values of the system parameters that one can consider: first, for $\beta=0$ and $\tilde{h}<2$, the origin is the only fixed point of the system and it is a stable node [Fig. 10, panel (a)], which becomes a saddle point when $\tilde{h}>2$ (not shown).

Second, for $\beta>0$ and $\tilde{h}<2$ [Fig. 10, panel (b)], two additional fixed points ($W_{-}$) appear in addition to the origin. We see that, at the origin ($R=0$) the eigenvalues of the Jacobian are $\lambda_{\pm;O}=-\omega_{0}/2\pm\omega_{0}\tilde{h}/4$, and therefore they are both real and negative if $\tilde{h}<2$, whereas the eigenvalues of the Jacobian for the $W_{-}$ points are $\lambda_{\pm;W_{-}}=-\omega_{0}/2\pm(\omega_{0}+3\tilde{h})/2$, which are always real and with opposite sign. The points $W_{-}$ are always saddle points, and therefore in the case $\tilde{h}<2$ only the origin is a stable point also for $\beta>0$.

Third, for $2<\tilde{h}<4$, two stable nodes ($W_{+}$) are born in pairs from the origin via a saddle-node bifurcation, after which the origin becomes a saddle point, independent of $\beta$ [Fig. 10, panel (c)]. This can be seen by looking at the eigenvalues of the Jacobian matrix: for the fixed points $W_{+}$ and $W_{-}$, the eigenvalues of the Jacobian are $\lambda_{\pm;W_{+}}=-\omega_{0}/2\pm\omega_{0}|3-\tilde{h}|/2$ and $\lambda_{\pm;W_{-}}=-\omega_{0}/2\pm\omega_{0}(3+\tilde{h})/2$. In this range of $\tilde{h}$, the point $W_{-}$ is always a saddle point, whereas the point $W_{+}$ is a stable node, with both eigenvalues of the Jacobian real and negative, whereas the origin (whose eigenvalues of the Jacobian matrix are $\lambda_{\pm;O}=-\omega_{0}/2\pm\omega_{0}\tilde{h}/4$, see above), is a saddle points when $\tilde{h}>2$. In this situation, for the trajectories flowing to the fixed point $W_{+}$, the imaginary part of the complex amplitude $A_{I}$ is suppressed ($A_{I}\rightarrow 0$), and the real part $A_{R}$ is stabilized to some nonzero value. This situations corresponds to squeezing.

Fourth, for $\tilde{h}>4$, two new stable attractors ($P$ and $Q$) are born from $W_{+}$ via a saddle-node bifurcation, after which the fixed point $W_{+}$ becomes a saddle point [Fig. 10, panel (d)]. For the fixed points $Q$ and $P$, the eigenvalues of the Jacobian matrix are $\lambda_{\pm;Q,P}=-\omega_{0}/2\pm(\omega_{0}/2)\sqrt{9-h^{2}/2}$. The fixed points $P$ and $Q$ are therefore stable nodes with both eigenvalues real and negative for $4<\tilde{h}<3\sqrt{2}$, and they are stable focuses (with eigenvalues with negative real part and nonzero imaginary part) for $\tilde{h}>3\sqrt{2}$.

One can see that the effect of having $\phi>0$ is to rigidly rotate the flow of the nonlinear equation by an angle of $\phi/2$. Therefore, in computing the position of the fixed points, one simply has to redefine the angles as $\varphi_{W_{\pm}}\rightarrow\varphi_{W_{\pm}}+\phi/2$ and $\varphi_{P,Q}\rightarrow\varphi_{P,Q}+\phi/2$, while the radial coordinates and the eigenvalues of the Jacobian matrix remain unaffected by $\phi$.

.2.2 Kerr nonlinearity

We here recall the stability diagram of the nonlinear Mathieu’s equation with Kerr nonlinearity in the resonant case (see for instance also Ref. Kidachi and Onogi (1997)). Focusing on the case of $\phi=0$, the nonlinear equations for the slow-varying amplitude are ($\tilde{\tau}=\omega_{0}\tau$)

[TABLE]

From Eq. (A12), one can see that there are two possible situations: when $\beta=0$, the origin is a stable node for $\tilde{h}<2$ and it is a saddle point for $\tilde{h}>2$, as in the case of Appendix .2.1. For $\beta>0$, the origin is the only fixed point for $\tilde{h}<2$ and it is a stable point. For $\tilde{h}>2$, the origin becomes a saddle node and two additional stable fixed points (which we denote by $P$) are born from the origin via a saddle-node bifurcation. As done in Appendix .2.1, we express the coordinates of the fixed point $P$ in polar coordinates, whose radial coordinate is

[TABLE]

and the angular coordinate is

[TABLE]

The eigenvalues of the Jacobian are found to be

[TABLE]

for a given fixed point. As in Eq. (A11), the eigenvalues of the Jacobian are independent of $\beta$. For the origin ($R=0$), one has as in Appendix .2.1 the eigenvalues $\lambda_{\pm;O}=-\omega_{0}/2\pm\omega_{0}\tilde{h}/4$. Instead, for the point $P$, one has $\lambda_{\pm;P}=-\omega_{0}/2\pm(\omega_{0}/2)\sqrt{5-\tilde{h}^{2}}$. The fixed points $P$ are the only fixed points in addition to the origin (in the resonant case $\gamma=2\omega_{0}$) that are found for all $\tilde{h}>2$, and they are stable nodes for $2<\tilde{h}<\sqrt{5}$ and stable focuses for $\tilde{h}>\sqrt{5}$. An example of the flow is shown in Fig. 11, in the case of $\beta>0$ and for the below-threshold case $\tilde{h}<2$ [panel (a)], and for the above-threshold case $\tilde{h}>2$ [panel (b)].

.3 Oscillation threshold and beats

In this appendix, we discuss the stability properties of the origin as a function of the system parameters for the system discussed in Sec. III, in the presence of both energy-preserving and dissipative coupling (see Sec. V). The results that we discuss in this section are valid also for $\beta=0$, since the stability properties of the origin are unaffected by the nonlinearity.

The informations regarding the position of the critical line for the bifurcation separating the region of the phase diagram in which the origin is a stable attractor (below threshold) from the one in which the origin is unstable (limit cycle, CIM or synchronization region) can be determined by studying the Jacobian matrix at the origin, i.e., $A=B=0$ (see Appendix .5). The eigenvalues of the Jacobian matrix at the origin are found to be

[TABLE]

The origin is a stable point if all the eigenvalues in Eq. (A16) have negative real part. In this case, the largest negative real part gives the decay rate to the trivial solution $A=B=0$ below the oscillations threshold: $\tau^{-1}_{\rm decay}=-g\,{\rm max}\{{\rm Re}[\lambda_{\pm,\pm}]\}=g\,{\rm min}|{\rm Re}[\lambda_{\pm,\pm}]|$. For sufficiently long times, below threshold, the decay to the trivial solution is then $x_{1}(t),x_{2}(t)\sim e^{-t/\tau_{\rm decay}}$. We can therefore have two regimes:

(i) Case $\tilde{r}>\tilde{\alpha}$ - Beating regime. By defining for convenience the function

[TABLE]

one can see that, in the $\tilde{h}$ vs. $\tilde{r}$ plane, the origin is stable point if $\tilde{h}<\tilde{h}_{\rm th}(\tilde{r},\tilde{\alpha},\phi)$, where $\tilde{h}_{\rm th}(\tilde{r},\tilde{\alpha},\phi)$ is identified by the contour line $Y(\tilde{h},\tilde{r},\tilde{\alpha},\phi)=8$, which yields

[TABLE]

Notice that, for generic $\phi$, one has from Eq. (A18) $\tilde{h}_{\rm th}(\tilde{r},\tilde{\alpha},\phi)\geq 2$, for all $\tilde{r}$, where the lower bound $\tilde{h}_{\rm th}=2$ is found for $\phi=0$, and it is independent of $\tilde{r}$, which is correctly the threshold condition $h/(2g)=1$ discussed in Eq. (11) at the parametric resonance ($\epsilon=0$). For a generic $\phi$, the threshold for parametric oscillations depends on the strength of the coupling.

The imaginary part of the eigenvalues in Eq. (A16), when present, is what determines the presence or absence of beats. From Eq. (A16), one can see that

[TABLE]

and therefore the angular frequency of the limit cycle at threshold is $\omega_{\rm LC,th}(\tilde{r},\tilde{\alpha},\phi)=\omega_{0}g\,\left|{\rm Im}[\lambda_{\pm,\pm}]\right|$.

It is now interesting to compare the results in Eq. (A18) and Eq. (A19) with the one discussed in the linear case in Sec. II.3. One can see that, above threshold for $\phi<\pi$, the origin is always an unstable point with $\left|{\rm Im}[\lambda_{\pm,\pm}]\right|>0$, and therefore the oscillators always display beats. In such situation, for $\beta=0$, no limit cycle can stabilize the amplitude of the beats and therefore the long-time dynamics of the oscillators is given by an exponential amplification with the beats superimposed. This situation corresponds to the one discussed in the linear case in Sec. II.3 by means of Floquet theorem, in the case in which the central instability region was present [see for example Fig. 1, panels (a) and (b)].

Instead, at $\phi=\pi$ and above the threshold identified by $\tilde{h}_{\rm th}(\tilde{r})=2\sqrt{\tilde{r}^{2}+1}$, one has $\left|{\rm Im}[\lambda_{\pm,\pm}]\right|=0$. In this case, the origin is unstable and amplification occurs without beats. This situation, for $\beta=0$, corresponds to the situation shown in Fig. 1, panel (c), in particular when the system is in the region in which the two outer instability regions overlap (around $\gamma=2\Omega_{r}\simeq 2\omega_{0}$). On top of these behaviours, found also in the linear case, the interplay between $\beta,g>0$ is what eventually stabilizes the beats with the presence of the limit cycle (for $\phi<\pi$) or the synchronization with the presence of additional stable attractors (for $\phi=\pi$) as the oscillation threshold [Eq. (A18)] is crossed.

(ii) Case $\tilde{r}<\tilde{\alpha}$ - CIM regime. In this case, in the regime of interest, the eigenvalues in Eq. (A16) are real. This first implies that, even when the origin becomes unstable, no limit cycle is found. The origin is an unstable point when $\lambda_{+,+}>0$, and therefore the threshold is identified by the condition

[TABLE]

from which for the requirement $\tilde{h}\in\mathbb{R}$ one obtains

[TABLE]

As discussed in Sec. V, above this threshold, the region with two stable fixed points is found.

.4 Alternative form of the coupling - Multiple-scale analysis

We report in this appendix the result of the calculation of the multiple-scale equations, using an alternative and commonly-used form of the coupling, and $\alpha=0$. We write the equation of motion [Eq. (22)] as

[TABLE]

[TABLE]

The multiple-scale equations governing the dynamics of the slow-varying amplitudes of $x_{1}$ and $x_{2}$ in Eqs. (A22) and (A23) are obtained as done for Eqs. (23). One obtains ($\tilde{\tau}=\omega_{0}\tau$)

[TABLE]

By comparing Eqs. (A24) with Eqs. (23), one can verify that the two set of equations describe the same dynamics if we perform in Eq. (23) the rotation $B_{R}\rightarrow-B_{I}$, $B_{I}\rightarrow B_{R}$ and redefine $\phi\rightarrow\pi+\phi$. This proves that, in the limit of small coupling, the long-time dynamics of the model with the coupling as in Eq. (22) and pump dephasing $\phi$ is equivalent to the long-time dynamics of the model with the coupling as in Eqs. (A22) and (A23) with pump dephasing $\pi+\phi$.

.5 Expression of the Jacobian matrix of the system in Eq. (23)

We here explicitly report for completeness the expression of the Jacobian matrix computed around a given point $A_{R},A_{I},B_{R},B_{I}$ from the system in Eq. (23). The Jacobian matrix can be compactly written as

[TABLE]

where we define the $2\times 2$ blocks as follows: one block for the oscillator A

[TABLE]

one phase-dependent block for the oscillator $B$

[TABLE]

and eventually the block describing the coupling between the oscillator $A$ and the oscillator $B$

[TABLE]