Coherent dynamics in frustrated coupled parametric oscillators

We consider a system of two identical single-mode parametric oscillators, with equal proper frequency ω₀, driven by an external pump field at frequency 2ω₀ and with amplitude h, injected into a parametric amplifier (PA) [41] as depicted in figure 1. The field inside each oscillator 1 (or 2) is identified by a classical variable x₁ (x₂). The two oscillators are coupled by a power-splitter coupling [37], which accounts for: (i) transmission coefficients c₁₁ and c₂₂ for oscillator 1 and 2, respectively, which renormalize the intrinsic loss of each oscillator, providing an overall loss rate that we denote by g, and (ii) coupling coefficients c₁₂ and c₂₁, which give the rate of energy exchange between the two oscillators. In this framework, the fields x₁ and x₂ are coupled according to the equation (see appendix A)

$$\begin{equation}\left(\begin{matrix}\hfill {\dot {x}}_{1}\hfill \\ \hfill {\dot {x}}_{2}\hfill \end{matrix}\right)={\omega }_{0}\left(\begin{matrix}\hfill 0\hfill & \hfill {c}_{12}\hfill \\ \hfill -{c}_{21}\hfill & \hfill 0\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \end{matrix}\right)\enspace ,\end{equation} \tag{ 1 }$$

where the dot denotes the time derivative. The dynamics of the two-oscillator system is described by a pair of coupled Mathieu's equations [38]

$$\begin{equation}\begin{cases}{\ddot {x}}_{1}+{\omega }_{0}^{2}\left[1+h\left(1-\beta \enspace {x}_{1}^{2}\right)\mathrm{sin}\left(2{\omega }_{0}t\right)\right]\enspace {x}_{1}+{\omega }_{0}\enspace g\enspace {\dot {x}}_{1}-{\omega }_{0}\enspace {c}_{12}\enspace {\dot {x}}_{2}=0\quad \hfill \\ {\ddot {x}}_{2}+{\omega }_{0}^{2}\left[1+h\left(1-\beta \enspace {x}_{2}^{2}\right)\mathrm{sin}\left(2{\omega }_{0}t\right)\right]\enspace {x}_{2}+{\omega }_{0}\enspace g\enspace {\dot {x}}_{2}+{\omega }_{0}\enspace {c}_{21}\enspace {\dot {x}}_{1}=0\quad \hfill \end{cases}\enspace .\end{equation} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Equation (2) also includes a second-order nonlinearity in the amplitude of the pump field (hereafter referred to as 'pump-depletion nonlinearity'), whose strength is quantified by β. Such a nonlinearity describes the fact that the intensity of the pump field inside each oscillator is depleted by x₁ and x₂, and in many experimental contexts captures the most relevant nonlinear process [38]. In general, the rate of energy flow between the two oscillators can be unbalanced, i.e., c₁₂ ≠ c₂₁, indicating dissipation in the coupling itself. Without loss of generality, one can parametrize the coupling coefficients in terms of an antisymmetric and symmetric part with respect to the exchange x₁ ↔ x₂ in equation (2): c₁₂ = r − α and c₂₁ = r + α, where the antisymmetric part r ⩾ 0 represents the energy-preserving component of the coupling, whereas the symmetric part α ⩾ 0 is the dissipative one.

The energy-preserving coupling r induces a coherent exchange of energy between the two oscillators. Its energy-preserving nature follows from the fact that the equations of motion (2), with β = 0, g = 0, and c₁₂ = c₂₁ = r, can be derived from the Hamilton's equations [19] starting from the Hamiltonian

$$\begin{equation}H=\frac{{p}_{1}^{2}+{p}_{2}^{2}}{2\enspace m}+\frac{1}{2}\enspace m\enspace {\omega }_{0}^{2}\left[1+\frac{{r}^{2}}{4}+h\enspace \mathrm{sin}\left(2{\omega }_{0}t\right)\right]\left({x}_{1}^{2}+{x}_{2}^{2}\right)+\frac{{\omega }_{0}r}{2}\left({p}_{1}{x}_{2}-{p}_{2}{x}_{1}\right)\enspace ,\end{equation} \tag{ 3 }$$

where p₁ and p₂ are the canonical momentum variables for x₁ and x₂, respectively. Such an Hamiltonian is analogous to that of a charged particle (charge q = mω₀r/2) moving on a two-dimensional plane identified by the spatial coordinates (x₁, x₂, z = 0), subject to a vector potential $\mathbf{A}={\left(-{x}_{2},{x}_{1},0\right)}^{\text{T}}$, where T denotes the transposition (for the details of the derivation, see appendix B).

Zoom In Zoom Out Reset image size

The dissipative coupling α, in contrast, introduces additional loss or gain terms [38], which give rise to the Ising-type coupling between the oscillators that is usually considered in the standard analysis of PO-CIMs [21–24, 26–30]. This coupling guides the convergence of the two-oscillator system to the desired Ising ground state [38]. In the long-time limit, the two oscillators will prefer to lock according to the sign of α: in-phase for 'ferromagnetic' coupling (α > 0), yielding the two 'ferromagnetic' configurations (00) or (ππ), or in anti-phase for 'anti-ferromagnetic' coupling (α < 0), yielding the two 'anti-ferromagnetic' configurations (0π) or (π0), where the first (second) label denotes the corresponding phase solution the first (second) oscillator.

We now review the effect of the interplay between r and α on the long-time dynamics of the system. We focus on the dynamics of the slow-varying amplitudes that modulate the fast-varying oscillations at half the pump frequency, which are instead integrated out, by employing the multiple-scale perturbative expansion in [38]. We take the intra-cavity loss g as a small expansion parameter, and identify the fast-varying and slow-varying time scales as t = 2π/ω₀ and τ = gt, respectively. By writing ${x}_{1}\left(t,\tau \right)=A\left(\tau \right)\enspace {\mathrm{e}}^{i{\omega }_{0}t}+{A}^{{\ast}}\left(\tau \right)\enspace {\mathrm{e}}^{-i{\omega }_{0}t}$ and ${x}_{2}\left(t,\tau \right)=B\left(\tau \right)\enspace {\mathrm{e}}^{i{\omega }_{0}t}+{B}^{{\ast}}\left(\tau \right)\enspace {\mathrm{e}}^{-i{\omega }_{0}t}$, where A(τ) and B(τ) are the complex amplitudes that encode the slow-varying dynamics, and by rescaling $\tilde {h}=h/g$, $\tilde {r}=r/g$, $\tilde {\alpha }=\alpha /g$, and $\tilde {\tau }={\omega }_{0}\tau$, one finds that A and B obey the following set of coupled first-order differential equations [38]:

$$\begin{align}\hfill \begin{aligned}\hfill \frac{\partial A}{\partial \tilde {\tau }}& =\frac{\tilde {h}}{4}\enspace {A}^{{\ast}}-\frac{\tilde {h}\enspace \beta }{4}\left(3{\vert A\vert }^{2}{A}^{{\ast}}-{A}^{3}\right)-\frac{A}{2}+\frac{\tilde {r}+\tilde {\alpha }}{2}\enspace B\hfill \\ \hfill \frac{\partial B}{\partial \tilde {\tau }}& =\frac{\tilde {h}}{4}\enspace {B}^{{\ast}}-\frac{\tilde {h}\enspace \beta }{4}\left(3{\vert B\vert }^{2}{B}^{{\ast}}-{B}^{3}\right)-\frac{B}{2}-\frac{\tilde {r}-\tilde {\alpha }}{2}\enspace A\enspace .\hfill \end{aligned}\end{align} \tag{ 4 }$$

Equation (4) can be further recast in terms of the real and imaginary parts of the complex amplitudes, A = A_R + iA_I and B = B_R + iB_I, where A_R (B_R) and A_I (B_I) are, respectively, the real and imaginary parts of A (B). The long-time dynamics is determined by the configuration and stability of the fixed points $\left({\overline{A}}_{\text{R}},{\overline{A}}_{\text{I}},{\overline{B}}_{\text{R}},{\overline{B}}_{\text{I}}\right)$ of equation (4), where the overline is used to denote the steady-state values of the quadratures. Note however that, sufficiently close to the oscillation threshold, A_I and B_I decay very quickly $\left({\overline{A}}_{\text{I}}={\overline{B}}_{\text{I}}=0\right)$ [38] due to the phase dependent amplification and squeezing in parametric oscillators, allowing to focus the discussion only on the dynamics of the real parts A_R and B_R.

While in general the configuration of the fixed points depends on the form of the nonlinearity, especially far from the amplification threshold, most of the interesting physics throughout this paper occurs close to the threshold, where nonlinear effects are negligible and the system is almost linear. The properties of the system at threshold can be therefore found by focussing on the spectrum of the Jacobian matrix around the origin A = B = 0, analysing specifically the eigenvalue with largest real part (λ_max, which we dub 'most efficient eigenvalue' from now on). Importantly, when the system exceeds the amplification threshold ${\tilde {h}}_{\text{th}}$, defined by the condition Re[λ_max] = 0⁺, the imaginary part of λ_max determines the frequency of the beats at threshold ω_B = ω₀g|Im[λ_max]| between the two oscillators. The beat frequency ω_B is the key observable to describe the behaviour of the system as the oscillators are driven above the oscillation threshold. Specifically, when ω_B = 0, A and B reach eventually constant values ${\overline{A}}_{\text{R}}$ and ${\overline{B}}_{\text{R}}$, and the system behaves as a time crystal, and can simulate Ising spins. Indeed, for positive (negative) ${\overline{A}}_{\text{R}}$, x₁ converges to the '0' ('π') solution, and analogously for x₂. Instead, for ω_B > 0, A and B display persistent coherent beats. The presence of the beats implies that each oscillator x_1,2 periodically flips between the '0' and 'π' solutions, and therefore the system neither obeys the Ising description, nor it behaves as a time crystal.

A concrete calculation of the phase diagram of the system in equation (4) is shown in figure 2, where we identify the following main phases:

• (a)  
  the sub-threshold phase, for a pump amplitude $\tilde {h}{< }{\tilde {h}}_{\text{th}}$, where the origin A = B = 0 is the only stable fixed point of equation (4);
• (b)  
  the Ising or PO-CIM phase, for $\tilde {h}{ >}{\tilde {h}}_{\text{th}}$, where two stable fixed points are found, the origin being unstable. Phase locking occurs at (00) or (ππ), for 'ferromagnetic' $\tilde {\alpha }{ >}0$, or at (0π) or (π0) for 'antiferromagnetic' $\tilde {\alpha }{< }0$ (not shown);
• (c)  
  the beating phase, for $\tilde {h}{ >}{\tilde {h}}_{\text{th}}$, where the long-time dynamics of the oscillators' amplitudes is attracted into a stable limit cycle;
• (d)  
  a phase with four stable fixed points, corresponding to all the four Ising configurations (00), (ππ), (0π), and (π0), in which the system behaves as two decoupled spins. This behaviour matches the one previously discussed in reference [21].

Slightly above the amplification threshold ${\tilde {h}}_{\text{th}}$, the PO-CIM phase is found for $\tilde {r}{< }\tilde {\alpha }$ (ω_B = 0), whereas the beating phase is found for $\tilde {r}{ >}\tilde {\alpha }$ (ω_B > 0). At threshold, the system therefore undergoes a transition between the PO-CIM and the beating behaviour when $\tilde {r}=\tilde {\alpha }$.

A schematic representation of a three-oscillator system is shown in figure 3. Now, the coupling matrix c between the oscillators, according to equation (1), is

$$\begin{equation}\left(\begin{matrix}\hfill {\dot {x}}_{1}\hfill \\ \hfill {\dot {x}}_{2}\hfill \\ \hfill {\dot {x}}_{3}\hfill \end{matrix}\right)={\omega }_{0}\left(\begin{matrix}\hfill 0\hfill & \hfill {c}_{12}\hfill & \hfill {c}_{13}\hfill \\ \hfill -{c}_{21}\hfill & \hfill 0\hfill & \hfill {c}_{23}\hfill \\ \hfill -{c}_{31}\hfill & \hfill -{c}_{32}\hfill & \hfill 0\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill {x}_{3}\hfill \end{matrix}\right)\enspace .\end{equation} \tag{ 5 }$$

The system is now described by a set of three coupled nonlinear Mathieu's equations, which we write in a compact form for the sake of simplicity (j, k = 1, 2, 3):

$$\begin{equation}{\ddot {x}}_{j}+{\omega }_{0}^{2}\left[1 + h\left(1 - \beta \enspace {x}_{j}^{2}\right) \mathrm{sin}\left(2{\omega }_{0}t\right)\right]{x}_{j}+{\omega }_{0}g{\dot {x}}_{j}-{\omega }_{0} \sum _{k\ne j}\mathrm{sgn}\left(k - j\right){c}_{jk}{\dot {x}}_{k}=0\enspace ,\end{equation} \tag{ 6 }$$

where sgn(⋅) denotes the sign function. From equation (6), one obtains the corresponding multiple-scale equations for the slow-varying amplitudes of the fields {x_j} as in equation (4). Here, we renormalize each element of the coupling matrix as ${\tilde {c}}_{jk}={c}_{jk}/g$, and to ease the notation, we use the symbols {X_j(τ)} of the amplitudes, such that ${x}_{j}\left(t,\tau \right)={X}_{j}\left(\tau \right){\mathrm{e}}^{i{\omega }_{0}t}+{X}_{j}^{{\ast}}\left(\tau \right){\mathrm{e}}^{-i{\omega }_{0}t}$. The equations for the complex amplitudes are then determined (j, k = 1, 2, 3):

$$\begin{equation}\frac{\partial {X}_{j}}{\partial \tilde {\tau }}=\frac{\tilde {h}}{4}{X}_{j}^{{\ast}}-\frac{1}{2}{X}_{j}-\frac{\tilde {h}\beta }{4}\left(3{\vert {X}_{j}\vert }^{2}{X}_{j}^{{\ast}}-{X}_{j}^{3}\right)+\frac{1}{2}\sum _{k\ne j}\mathrm{sgn}\left(k-j\right){\tilde {c}}_{jk}{X}_{k}\enspace .\end{equation} \tag{ 7 }$$

As in section 2, we decompose the coupling c_jk in equations (5)–(7) between any two oscillators in terms of an energy-preserving (antisymmetric) and dissipative (symmetric) part, respectively r_jk and α_jk. Due to this increase of parameter space with respect to the case in section 2 (the coupling matrix now has in general six independent components), we focus on a specific choice of the coupling matrix, with the ambition to highlight the role of frustration in the dissipative components of the coupling matrix. We choose r_jk = r for all j and k, and introduce two different dissipative couplings: α₁₂ = η and α₁₃ = α₂₃ = α, so that the coupling matrix in equation (5) reads

$$\begin{equation}\left(\begin{matrix}\hfill {\dot {x}}_{1}\hfill \\ \hfill {\dot {x}}_{2}\hfill \\ \hfill {\dot {x}}_{3}\hfill \end{matrix}\right)={\omega }_{0}g\left(\begin{matrix}\hfill 0\hfill & \hfill \tilde {r}+\tilde {\eta }\hfill & \hfill \tilde {r}+\tilde {\alpha }\hfill \\ \hfill -\tilde {r}+\tilde {\eta }\hfill & \hfill 0\hfill & \hfill \tilde {r}+\tilde {\alpha }\hfill \\ \hfill -\tilde {r}+\tilde {\alpha }\hfill & \hfill -\tilde {r}+\tilde {\alpha }\hfill & \hfill 0\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {x}_{1}\hfill \\ \hfill {x}_{2}\hfill \\ \hfill {x}_{3}\hfill \end{matrix}\right)\enspace .\end{equation} \tag{ 8 }$$

In the rest of the paper, our goal is to study the physics of the system near threshold as a function of the coupling parameters $\tilde {r}$, $\tilde {\alpha }$, and $\tilde {\eta }$, where as before all quadratures are real $\left({\overline{X}}_{1,\text{I}}={\overline{X}}_{2,\text{I}}={\overline{X}}_{3,\text{I}}=0\right)$. Before discussing this general case, we first focus on two main configurations of interest for the coupling matrix in equation (8). Namely, assuming $\tilde {r},\tilde {\alpha }{\geqslant}0$: (i) the non-frustrated case, for $\tilde {\eta }=\tilde {\alpha }$, and (ii) the fully-frustrated case, for $\tilde {\eta }=-\tilde {\alpha }$. The reason why we focus on these two fine-tuned cases is because they are, on one hand, easily analytically tractable, and on the other hand, they capture the dramatic effect of frustration in the dissipative coupling. We discuss the general case later in section 3.4.

Zoom In Zoom Out Reset image size

First, we analytically discuss the threshold properties of the non-frustrated network, for $\tilde {\eta }=\tilde {\alpha }$ in equation (8). We analyse separately the cases of $\tilde {r}{ >}\tilde {\alpha }$ and $\tilde {r}{< }\tilde {\alpha }$.

For $\tilde {r}{ >}\tilde {\alpha }$, we have ${\lambda }_{\mathrm{max}}=-1/2+\tilde {h}/4-{\mathrm{e}}^{i\pi /3}\enspace F\left(\tilde {r},\tilde {\alpha }\right)+{\mathrm{e}}^{-i\enspace \pi /3}\enspace G\left(\tilde {r},\tilde {\alpha }\right)$, where

$$\begin{equation}F\left(\tilde {r},\tilde {\alpha }\right)=\frac{{\tilde {r}}^{2}-{\tilde {\alpha }}^{2}}{4\enspace G\left(\tilde {r},\tilde {\alpha }\right)}\quad G\left(\tilde {r},\tilde {\alpha }\right)=\frac{1}{2}\enspace {\left[\left({\tilde {r}}^{2}-{\tilde {\alpha }}^{2}\right)\left(\tilde {r}+\tilde {\alpha }\right)\right]}^{1/3}\enspace .\end{equation} \tag{ 9 }$$

Since λ_max is complex, beats are found. From the imaginary parts of λ_max, one can find the expression of the beat frequency:

$$\begin{equation}{\omega }_{\text{B}}=\frac{g{\omega }_{0}\sqrt{3}}{4}{\left({\tilde {r}}^{2}-{\tilde {\alpha }}^{2}\right)}^{1/3}\left[{\left(\tilde {r}-\tilde {\alpha }\right)}^{1/3}+{\left(\tilde {r}+\tilde {\alpha }\right)}^{1/3}\right]\enspace .\end{equation} \tag{ 10 }$$

In particular, for $\tilde {r}$ approaching $\tilde {\alpha }$, the frequency of the beats reduces towards zero with the critical behaviour of ${\omega }_{\text{B}}\sim {\left(\tilde {r}-\tilde {\alpha }\right)}^{1/3}$, differently from the critical exponent of 1/2 in the two-oscillator case [38].

When $\tilde {r}{< }\tilde {\alpha }$, we find that ${\lambda }_{\mathrm{max}}=-1/2+\tilde {h}/4+F\left(\tilde {\alpha },-\tilde {r}\right)+G\left(\tilde {\alpha },-\tilde {r}\right)$. Now, λ_max is real, and above the oscillation threshold, parametric amplification occurs without beats, ω_B = 0. In this case, the phase-locked steady-state oscillations correspond to the two 'ferromagnetic' configurations (000) or (πππ), as shown in the left panel figure 4 (in the figure specifically for $\tilde {r}=0$). As one may expect, the behaviour of the non-frustrated network is qualitatively the same as the behaviour of two coupled oscillators (section 2).

Zoom In Zoom Out Reset image size

We now move to the case of the fully-frustrated network, i.e., $\tilde {\eta }=-\tilde {\alpha }$ in equation (8). By proceeding as before, we find that, for $\tilde {r}{ >}\tilde {\alpha }$, the most efficient eigenvalue is ${\lambda }_{\mathrm{max}}=-1/2+\tilde {h}/4+F\left(\tilde {r},-\tilde {\alpha }\right)-G\left(\tilde {r},-\tilde {\alpha }\right)$. Now, in stark contrast to the non-frustrated case, λ_max is real, and above the oscillation threshold parametric amplification occurs without beats, ω_B = 0.

For $\tilde {r}{< }\tilde {\alpha }$, instead, λ_max reads as ${\lambda }_{\mathrm{max}}=-1/2+\tilde {h}/4+{\mathrm{e}}^{i\enspace \pi /3}\enspace F\left(\tilde {\alpha },\tilde {r}\right)+{\mathrm{e}}^{-i\enspace \pi /3}\enspace G\left(\tilde {\alpha },\tilde {r}\right)$. Therefore, the parametric oscillation occurs with beats, where the beat frequency is

$$\begin{equation}{\omega }_{\text{B}}=\frac{g{\omega }_{0}\sqrt{3}}{4}{\left({\tilde {\alpha }}^{2}-{\tilde {r}}^{2}\right)}^{1/3}\left[-{\left(\tilde {\alpha }-\tilde {r}\right)}^{1/3}+{\left(\tilde {\alpha }+\tilde {r}\right)}^{1/3}\right]\enspace .\end{equation} \tag{ 11 }$$

One can see by inspection that, when $\tilde {r}{< }\tilde {\alpha }$, the presence of the limit cycle at threshold makes the system periodically flip between six possible phase-locked configurations, namely, (0ππ), (0π0), (000), (π00), (π0π), and (πππ), corresponding to the six degenerate ground-state configurations of the frustrated Ising model. One of these configurations stabilizes the long-time dynamics only when $\tilde {r}=0$ (see figure 4, right panel). The frustration of one of the dissipative components of the coupling has therefore a dramatic effect on the coherent dynamics of the network. In the non-frustrated case, the system at threshold converges to a phase-locked configuration for $\tilde {r}{< }\tilde {\alpha }$, and displays persistent beats otherwise. Instead, the behaviour of the fully-frustrated network is reversed with respect to the non-frustrated case: it presents beats for $\tilde {r}{< }\tilde {\alpha }$, and converges to phase-locked oscillations otherwise. For fixed $\tilde {\alpha }$, the frequency of the beats increases linearly from zero for small $\tilde {r}$, and goes to zero as $\tilde {r}$ approaches $\tilde {\alpha }$ with the same critical exponent 1/3 as before.

We now expand the discussion to consider the general case of equation (8), which interpolates between the non-frustrated and the frustrated cases, and generalizes the analysis in sections 3.2 and 3.3. Here, one can find the frequency of the beats at threshold ω_B from the imaginary part of the most efficient eigenvalue, for a fixed $\tilde {\alpha }$, as a function of $\tilde {r}$ and $\tilde {\eta }$, and discern regions in parameter space where phase-locked oscillations or beats are observed. We present our findings in figure 5 for both theory [panel (a)] and low-level simulation of the experiment [panel (b)], which will be discussed in detail in the next section. To ease the comparison between the analytical prediction and the low-level simulation results, we express the frequency of the beats in units of α. This makes the frequency of the beats ω_B/α be a function of r/α multiplied by pure numbers and independent of g (see sections 3.2 and 3.3, and appendix C).

Zoom In Zoom Out Reset image size

Panel (a) of figure 5 shows ω_B/α in the r/α vs η/α plane, where the red lines mark the special cases of fully frustrated and totally non-frustrated Ising coupling (η/α = ±1) that were discussed in sections 3.2 and 3.3. Green boundaries separate the regions where phase-locking is found (ω_B = 0, dark blue) from those where persistent beating is manifested (ω_B > 0, other colours). We see that, for η/α > 0, a region of beats is found when the energy-preserving coupling dominates over the dissipative coupling (r > α), and phase-locking is found otherwise, similar to the two-oscillator analysis and to the non-frustrated case. However, when η/α < 0, a 'tooth-shaped' region of beats appears when the dissipative coupling dominates (r < α), and phase locking is found otherwise. As r is lowered towards r = 0, the width of the tooth region of beats decreases until it collapses to a point when r = 0, at η = −α, i.e., the fully-frustrated point of section 3.3.

This peculiar behaviour of beating for small r/α may have important implications in the context of PO-CIMs. For example, fixing r and scanning η from positive to negative allows to interpolate between the ferromagnetic non-frustrated (η = α), and fully frustrated (η = −α) Ising models, where the transition to the fully-frustrated case occurs at η/α = −1. At the transition point, the Ising gap (i.e., the energy difference between the ground- and first-excited configurations) closes, causing the multiplicity of the ground state to increase (in our case, from two to six, see blue and red curves in figure 6). Our findings show that, in the system of three parametric oscillators coupled by the matrix in equation (8), any vanishingly small energy-preserving coupling induces coherent beating between the oscillators as the Ising gap becomes vanishingly small, preventing the system from converging to the Ising ground-state configuration.

Zoom In Zoom Out Reset image size

In addition, we find that, while our three-oscillator system correctly behaves as a PO-CIM at the fully-frustrated point and in the phase-locking 'Ising' region in figure 5, in the other phase-locking ('Non Ising') region, the system does not yield the expected Ising behaviour. Indeed, the Ising model predicts a four-fold degenerate ground-state when η < −α (see red curve in figure 6). However, in the 'Non Ising' region in figure 5, the oscillator system slightly above the threshold converges only to two fixed points. In particular, we find the following behavour:

• for r = 0 and r = α, the two fixed points are found on the X_3,R = 0 plane, implying that lim_t→∞x₃(t) = 0, and the other two oscillators converge to (0π) or (π0). Clearly, the suppression of one oscillator in the long-time limit is not a valid Ising configuration;
• for 0 < r < α, the two fixed points correspond to the states (0ππ) and (π00), which are only two of the four ground states of the Ising model;
• for r > α, the two fixed points correspond to the states (0π0) and (π0π). These two configurations, for η < −α, are two of the four ground states, as before, but for −α < η < 0, they correspond to two of the four excited states of the Ising model.

This finding hints that frustration may cause phase-locked oscillation at threshold that however transcend the Ising description. Before concluding, we stress that the presence and details of the beating region for vanishingly small energy-preserving coupling strongly depend on the form of the coupling matrix, as well as on the number of oscillators. Indeed, while frustrated spin models with a small number of spins, simulated with PO-CIMs, have been experimentally studied in previous work [42, 43], coherent beats were not reported. This fact can be due to both the different coupling matrix considered in [42, 43], as well as to the fact that the energy-preserving coupling was possibly suppressed. A deeper analysis on how our results translates to general number of coupled oscillators and general coupling topology requires further analysis of larger spin models, which is beyond the scope of the present manuscript, and it is left for future work.

We consider a multimode cavity where each temporal slot acts as an independent parametric oscillator, dynamically coupled to the other modes. In our simulation, the field inside the cavity propagates as illustrated in the block diagram in figure 7. At each round trip inside the cavity, the time signal is partitioned into N = 3 time slots, and parametrized as a three-dimensional vector $\mathbf{S}={\left({S}_{1},{S}_{2},{S}_{3}\right)}^{\text{T}}$. In each such interval, the field is assumed to vary slowly and is amplified independently of the other time slots. This is a reasonable assumption since the parametric gain is an instantaneous process and the pump for each time slot is uncoupled from that of the other slots. Thus, each time slot defines a distinct parametrically driven mode that without coupling evolves independently from the others. Furthermore, additive noise is fed inside the cavity at each round trip through an output coupler device to simulate thermal noise or vacuum fluctuations. During the first round trip, before being injected into the PA, the signal inside the cavity consists of noise alone.

Zoom In Zoom Out Reset image size

A round trip inside the cavity is identified by the following steps: first, the pump field and the signal are injected into the PA. Importantly, since our goal is to probe the linear, near threshold, properties of the system, the pump intensity is set slightly above the oscillation threshold (section 3). At the output of the PA, the residual pump field is blocked (dark grey parallel lines in figure 7) and the signal is injected into a coupler, which splits the field according to transmission and reflection coefficients T_c = 1/4 and ${R}_{\mathrm{c}}=\sqrt{1-{T}_{\mathrm{c}}^{2}}$. The transmitted signal T_cS is sent into a coupling mechanism, which implements the coupling matrix c of equation (8). In a future experiment, such a coupling mechanism will be implemented by an FPGA. After the coupling, the coupled signal F = T_ccS and the reflected signal R_cS are combined on another coupling device, again with transmission and reflection coefficients T_c and R_c. At this step, the reflected signal T_cR_cS + R_cF is blocked (contributing to the overall cavity losses), and the transmitted signal ${R}_{\mathrm{c}}^{2}\mathbf{S}+{T}_{\mathrm{c}}\mathbf{F}$ is fed back into the cavity. Last, the output coupler with transmission and reflection coefficients T_out = 1/2 and ${R}_{\text{out}}=\sqrt{1-{T}_{\text{out}}^{2}}$ allows to couple out the field ${T}_{\text{out}}\left({R}_{\mathrm{c}}^{2}\mathbf{S}+{T}_{\mathrm{c}}\mathbf{F}\right)$ from the cavity and analyse it both it in time and frequency, and the remaining field ${R}_{\text{out}}\left({R}_{\mathrm{c}}^{2}\mathbf{S}+{T}_{\mathrm{c}}\mathbf{F}\right)$, fed with noise, is input together with the pump field into the PA to start the next round trip.

For a given set of parameters, we run the simulation over a sufficient number of round trips n_rt in order for the oscillation to reach a steady state. We then examine the steady-state dynamics to obtain the slow-varying amplitude of all oscillators and numerically extract the frequency of the beats.

The numerical frequency of the beats is measured by computing the fast Fourier transform (FFT) of the signal at the output to identify the frequency component with largest amplitude, for different values of r/α and/or η/α. In the numerics, the frequencies obtained from the FFT are given in units of 1/τ_sim, where the total simulation time is τ_sim = n_rtτ_rt, and round-trip time τ_rt in our numerical context is an arbitrary time scale. In order to quantitatively compare the numerical results with the analytical prediction, we express the frequency of the beats in units of α (section 3.4) and use τ_rt = 0.0195 (see appendix D for more details).

Figure 8 shows the frequency of the beats as a function of r/α for both the numerical simulation of the experiment and the theoretical analysis. We compare the results in the special cases of the non-frustrated network (η = α) and the fully-frustrated one (η = −α). As evident, the numerical results in both cases agree well with the theoretical curves, with some slight deviations especially in the fully-frustrated cases, which we ascribe to the difficulty of correctly estimating the oscillation threshold and therefore choosing the proper value of h, due to both noise and nonlinearities. Indeed, it has been shown that pump-depletion nonlinearity tends to lower the beating frequency (divergence of the period of the beats), eventually inducing phase-locking as h is increased above the oscillation threshold [37].

Zoom In Zoom Out Reset image size

Last, in figure 5, panel (b), we evaluate the beating frequency ω_B/α as a function of r/α and η/α, in order to verify the theoretical phase diagram in panel (a). The analytical phase boundary, marked by the green line, which is the same for both panels of figure 5, is superimposed to the numerical phase diagram to ease comparison between theory and simulated experiment. As evident, our numerical data agree exceptionally well with the analytical prediction also in the interpolating case.