Coherent dynamics in frustrated coupled parametric oscillators

We explore the coherent dynamics in a small network of three coupled parametric oscillators and demonstrate the effect of frustration on the persistent beating between them. Since a single-mode parametric oscillator represents an analog of a classical Ising spin, networks of coupled parametric oscillators are considered as simulators of Ising spin models, aiming to efficiently calculate the ground state of an Ising network - a computationally hard problem. However, the coherent dynamics of coupled parametric oscillators can be considerably richer than that of Ising spins, depending on the nature of the coupling between them (energy preserving or dissipative), as was recently shown for two coupled parametric oscillators. In particular, when the energy-preserving coupling is dominant, the system displays everlasting coherent beats, transcending the Ising description. Here, we extend these findings to three coupled parametric oscillators, focusing in particular on the effect of frustration of the dissipative coupling. We theoretically analyze the dynamics using coupled nonlinear Mathieu's equations, and corroborate our theoretical findings by a numerical simulation that closely mimics the dynamics of the system in an actual experiment. Our main finding is that frustration drastically modifies the dynamics. While in the absence of frustration the system is analogous to the two-oscillator case, frustration reverses the role of the coupling completely, and beats are found for small energy-preserving couplings. This result is specifically relevant to Ising simulators, where residual energy-preserving couplings may lead to coherent beats and prevent the system from converging to the Ising ground state.


Introduction
Parametric oscillators are a viable experimental platform to study the physics of time crystals, i.e., systems that can spontaneously break time translational symmetry [1,2]. The possibility of the existence of such a phase of matter at equilibrium was first proposed in 2012 by Frank Wilczek and collaborators [3,4], both for quantum and classical systems. The original proposal evokes the possibility for a system to break continuous time translational symmetry, in analogy with the formation of space crystals in condensed matter where space translational symmetry is broken. Shortly after arXiv:2005.03562v1 [nlin.CD] 7 May 2020 its proposal, it became clear that this kind of time-crystalline phase cannot exist at equilibrium [5][6][7]. However, following Wilczek's original idea, it was understood that time crystals can be realized out of equilibrium, in periodically-driven system, also referred to as Floquet systems. This new type of time crystals, dubbed Floquet time crystals, accounts for the fact that, under certain conditions, a periodically-driven system can break the discrete time translational symmetry enforced by the external drive [8][9][10][11][12][13][14][15][16][17][18]: Instead of merely following the external drive, the system undergoes a periodic motion at a frequency that is different from that of the drive (see ref. [1] for a review).
The periodically driven single-mode classical parametric oscillator is the canonical example of period-doubling instability (see refs. [19,20] for an introduction), and represents the simplest case of a classical Floquet time crystal. Indeed, when excited above the amplification threshold, the parametric oscillator oscillates at half the frequency of the drive and admits only two distinct phase solutions, dubbed "0" and "π", with a relative shift in time by one period of the drive. One of the two solutions is chosen by the system depending on the initial conditions, a phenomenology analogous to a spontaneous breaking of a Z 2 (Ising) symmetry. Because of this, a single degenerate parametric oscillator may be regarded as a classical bit, or an Ising spin, where the two states "up" or "down" of the spin are given by the two distinct "0" and "π" solutions. Exploiting this property, networks of many coupled parametric oscillators have been proposed as a platform, called coherent Ising machine (CIM) [21], to simulate the behaviour of a network of many coupled Ising spins. Such a machine, whose experimental realization has been reported in refs. [22][23][24], is envisioned to solve the NP-hard problem of finding the ground state of the classical Ising model [25].
In the last years, the analysis of various issues related to the computational performance of CIMs has been the focus of a remarkable amount of work [26][27][28][29][30][31][32]. While the underlying assumption in CIMs is that a system of coupled parametric oscillators behaves as a set of coupled Ising spins, we pointed out recently that already a pair of coupled parametric oscillators may display a much richer dynamics, beyond the Ising description, depending on the nature of the coupling (energy-preserving or dissipative) between the oscillators [33,34]. Specifically, we studied in detail both theoretically and experimentally in a radio-frequency experiment a pair of coupled parametric oscillators, which is the minimal system to explore nontrivial coupling effects. Our main finding was that, when driven above the amplification threshold at the parametric resonance condition, the two oscillators can either display persistent coherent beats when the coupling is mostly energy preserving, or behave as a CIM [21] when the coupling is mostly dissipative.
The existence of such a nontrivial dynamics in just a pair of coupled parametric oscillators opens the question on how the nature of the coupling affects the dynamics of a larger network composed by more than two parametric oscillators, with potential implications in the specific context of CIMs, and also in the broad view of exploiting large-scale networks of coupled parametric oscillators to realize classical many-body time crystals [35]. Motivated by these perspectives, we present in this paper a detailed theoretical and numerical analysis of three coupled parametric oscillators, which is the minimal system where nontrivial connectivity effects, such as frustration, can be studied. The coupling between any two oscillators is parametrized by two coupling components -energy-preserving and dissipative. The main focus of this paper is to analyze the effects of frustration in the dissipative components of the coupling, which turns out to be dramatic.
To reach this goal, we model each parametric oscillator as a classical variable and describe the system by three coupled nonlinear Mathieu's equations [34], in the presence of an external pump, intrinsic dissipation, and pump depletion nonlinearity, to analyze the phase diagram of the system for different values of the system parameters, as detailed hereon. Our theoretical predictions are confirmed by low-level numerical simulation of the field propagation within the parametric oscillators both in time and space, as close as possible to an actual experimental setup. Our numerical scheme emulates directly the dynamics of the field inside a cavity with parametric gain, with no explicit mention of the equations of motion studied in our analytical model.
Our main finding is that frustration totally inverts the dynamical picture of the coupled system. While in the absence of frustration the system behaves similar to the two-oscillator case, where beats are observed only when the energy-preserving coupling is larger than the dissipative one, in the presence of frustration we find two main differences: First, the system shows coherent everlasting beats for small energy-preserving couplings. This finding can be reasoned by the fact that a frustrated system cannot distinguish between two (or more) degenerate Ising states that are found when the coupling is purely dissipative. Thus, any non-vanishing value of the energypreserving coupling induces beating between those degenerate states. Second, for large energy-preserving couplings and large frustration, the network converges to a phaselocked oscillation, which however is not the Ising ground state. This paper is organized as follows. In section 2, we shortly review the case of two coupled parametric oscillators, introducing our model and notations. We then present our theoretical analysis for the case of three coupled parametric oscillators in section 3. We discuss in section 4 a possible experimental implementation of our system, and present the results of the low-level numerical simulation of such an experiment. We then draw our conclusions in section 5, and report some relevant details on the calculations in the appendixes.

Two parametric oscillators
Before moving to the analysis of three-coupled oscillators, let us introduce the relevant notation and analytical tools, via a review of the simple case of two degenerate coupled parametric oscillators, summarizing the main findings of refs. [33,34].

Ising Beats
Sub-threshold Figure 1. Schematic representation of the two-parametric-oscillator system. Each parametric oscillator is described by a classical field x k , with k = 1, 2, driven by an external pump h(t) = h sin(2ω 0 t) injected into a parametric amplifier (PA). The coupling between the oscillators is described in the general case by a coupling matrix that accounts for (i) transmittance coefficients c 11 and c 22 , which renormalize the intrinsic losses of the oscillators, and (ii) coupling coefficients c 12 and c 21 , which determine the rate of energy flow from oscillator 1 to 2, and vice versa.

Model and notation
We consider a system of two identical single-mode parametric oscillators, with equal proper frequency ω 0 , driven by an external pump field at frequency 2ω 0 and with amplitude h, injected into a parametric amplifier (PA) [36] as depicted in figure 1. The field inside each oscillator 1 (or 2) is identified by a classical variable x 1 (x 2 ). The two oscillators are coupled by a power-splitter coupling [33], which accounts for: (i) transmission coefficients c 11 and c 22 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 12 and c 21 , which give the rate of energy exchange between the two oscillators. In this framework, the fields x 1 and x 2 are coupled according to the equation (see Appendix A) where the dot denotes the time derivative. The dynamics of the two-oscillator system is therefore described by a pair of coupled Mathieu's equations [34]     ẍ 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 1 and x 2 , and in many experimental contexts captures the most relevant nonlinear process (see ref. [34] for a critical discussion). In the general case, the rate of energy flow between the two oscillators can be unbalanced, i.e., c 12 = c 21 , indicating dissipation in the coupling itself. Without loss of generality, one can parametrize the coupling coefficients as c 12 = r − α and c 21 = r + α, where r ≥ 0 represents the energy-preserving component of the coupling, whereas α ≥ 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 12 = c 21 = r, can be derived from the Hamilton's equations [19] starting from the Hamiltonian where p 1 and p 2 are the canonical momentum variables for x 1 and x 2 , respectively. Such an Hamiltonian is analogous to that of a charged particle (charge q = mω 0 r/2) moving on a two-dimensional plane identified by the spatial coordinates (x 1 , x 2 , z = 0), subject to a vector potential A = (−x 2 , x 1 , 0) T , where T denotes the transposition (for the details of the derivation, see Appendix B). The dissipative coupling α, in contrast, introduces additional loss or gain terms [34], which give rise to the Ising-type coupling between the oscillators, and guides the convergence of the two-oscillator system to the desired Ising ground state [34], as required by a CIM [21]. 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 "antiferromagnetic" 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.

Multiple-scale expansion -Phase-locking and beats
We now review the effect of the interplay between r and α on the long-time dynamics of the system. To reach this goal, we employ a perturbative analysis that gives the solution in terms of the slowly-varying dynamics of the unperturbed oscillators. That technique is known as multiple-scale expansion in the context of nonlinear dynamics [37], and as the slow-varying amplitude approximation in quantum and nonlinear optics. We here review the main steps, referring the interested reader to ref. [34] for complete description.
We take the intra-cavity loss g as a small expansion parameter, and identify the fast-varying time scale as t = 2π/ω 0 -the period of the fast oscillations at half the pump frequency, whereas the slow-varying time scale is identified by τ = gt. One then studies the dynamics only of the slow-varying degrees of freedom, integrating out the fast-varying ones.
In the unperturbed case, the two fields x 1 and x 2 are then explicitly expressed by separating the two time scales as x 1 (t, τ ) = A(τ ) e iω 0 t + A * (τ ) e −iω 0 t and x 2 (t, τ ) = B(τ ) e iω 0 t +B * (τ ) e −iω 0 t , for some complex amplitudes A(τ ) and B(τ ), in which the longtime dynamics is encoded, and A * (B * ) is the complex conjugate of A (B). One then studies the equations of motion (2) in the limit in which all the coupling constants affect the slow-varying motion alone. In order to proceed with the perturbative expansion, h, r, and α are taken proportional to g, i.e., one definesh = h/g,r = r/g, and α = α/g. By defining the dimensionless time asτ = ω 0 τ , one finds that the slow-varying complex amplitudes A and B obey the following set of coupled first-order differential equations [34]: Equation (4) (4), and by studying whether they are stable or unstable by looking at the eigenvalues of the Jacobian matrix J(A, B) of the system in equation (4) around a specific fixed point. Note however that, sufficiently close to the oscillation threshold, the imaginary part of both oscillation amplitudes decays very quickly (A I = B I = 0) [34] 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 found by studying the spectrum of the Jacobian matrix around the origin A = B = 0. Such an analysis is sufficient to understand if the parametric oscillators converge to a steady-state oscillation, corresponding to reaching the "0" or "π" solutions for each oscillator, or not. It therefore provides informations on whether the overall system mimics an Ising network, and behaves as a time crystal, or not [34].
The general procedure is as follows: For a given set of parametersh,r andα, we look at the eigenvalue of the Jacobian matrix J(0, 0) with largest real part (λ max , which we dub "most efficient eigenvalue" from now on). For a pump amplitude below threshold h <h th , we have Re[λ max ] < 0, implying that oscillations decay: lim t→∞ x 1,2 (t) = 0. The value ofh th is defined as the value ofh such that Re[λ max ] = 0. Forh >h th , we have Re[λ max ] > 0, indicating that parametric amplification sets in, and the amplitudes A and B of the fields x 1,2 (t) exponentially grow in time, until eventually stabilized by the pump depletion nonlinearity.
Importantly, forh =h th , the imaginary part of λ max determines the frequency of the beats at threshold ω B = ω 0 g|Im[λ max ]| between the two oscillators, which modulate the aforementioned exponential growth of the amplitudes. The beat frequency ω B is the key observable to describe the behaviour of the system as the oscillators are driven above the oscillation threshold: • When ω B = 0, A and B reach eventually constant values A R and B R , due to the presence of the nonlinearity. Hence, in the long-time limit, x 1 (t) ∼ 2A R cos(ω 0 t) and similar for x 2 (t). This is the situation where the system behaves as a time crystal, and can simulate Ising spins. Indeed, if A R > 0, x 1 converges to the "0" solution, whereas if A R < 0, x 1 converges to the "π" solution, and analogously for x 2 . The two-oscillator system then converges to some of the four possible "Ising" configurations (00), (ππ), (0π), and (π0), depending on the initial conditions and on the coupling parameters; • Instead, for ω B > 0, A and B display persistent coherent beats. In such a situation, one has, in the long-time limit, 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, since the system periodically jumps between all the four Ising configurations (00), (ππ), (0π), and (π0). In the phase diagram of figure 2, we identify the following main phases:

Phase diagram
(i) The sub-threshold phase, for a pump amplitudeh <h th , where the origin A = B = 0 is the only stable fixed point of equation (4); (ii) The Ising or CIM phase, forh >h th , defined by the presence of two stable fixed points, related by the inversion symmetry A ↔ −A and B ↔ −B of equation (4), the origin being unstable. The solutions x 1 (t) and x 2 (t) display a fast oscillation at half the pump frequency and are phase-locked either at (00) or (ππ), for "ferromagnetic"α > 0, or either at (0π) or (π0) for "antiferromagnetic"α < 0 (not shown); (iii) The beating phase, forh >h th , in which no stable fixed point is found, and the long-time dynamics is attracted into a stable limit cycle, which is identified by an attractive isolated closed orbit [20] in the space of the amplitudes A and B. In this phase, x 1 and x 2 display fast oscillations at half the pump frequency whose amplitude is modulated (beats) with a frequency that depends onr,α, andh [34].
(iv) 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

Ising Beats
Sub-threshold  (4) is marked by red curves. The phase diagram consists of four main phases, characterized by different configurations of the fixed points, as shown in right panels: (i) below the oscillation threshold, where only the origin A = B = 0 is a stable attractor; (ii) the Ising or CIM region, where the system has two stable fixed points, corresponding to the two ground-state Ising solutions (00) and (ππ) (for "ferromagnetic"α > 0); (iii) a phase in which a stable limit cycle stabilizes the long-time dynamics, and the system displays coherent beats; (iv) a phase with four stable fixed points, corresponding to both ground-state (00) and (ππ), and excited-state (0π) and (π0) configurations. The red dashed line in the phase diagram is the boundary for the oscillation thresholdh th . Other phases with more than four fixed points [34] are not labelled and not relevant for the present discussion, and additional unstable fixed points different from the origin are not shown for the sake of clarity.
decoupled spins. This behaviour matches the one previously discussed in ref. [21]. As in phase (ii), the amplitude of the fast oscillations converge to a constant value, and the system behaves as a time crystal.
Slightly above the amplification thresholdh th , the CIM phase is found forr <α (ω B = 0), whereas the beating phase is found forr >α (ω B > 0). At threshold, the system therefore undergoes a transition between the CIM and the beating behaviour whenr =α.

Three coupled parametric oscillators
The results discussed in the previous section pointed out the existence of a persistent coherent beating dynamics in coupled parametric oscillators, not considered in the standard analysis of CIMs. A natural question that arises is how the presence of such Figure 3. Schematic representation of the three-parametric-oscillator system. Similar to figure 1, the parametric oscillators are described by a classical field x k (k = 1, 2, 3), driven by the same external pump h(t) = h sin(2ω 0 t). The coupling between the oscillators, describing the connectivity of the system, includes all the mutual couplings c jk for j, k = 1, 2, 3. The transmittance coefficients c jj for j = 1, 2, 3 have been omitted from the figure for the sake of simplicity.
a dynamics affects a more structured network with more than two coupled parametric oscillators. Here, we begin to address this question by extending the previous discussion to the case of three degenerate coupled parametric oscillators, which is the simplest configuration where one can systematically study the role of connectivity. Specifically, our main focus is to study how frustration of the dissipative coupling (which reflects the underlying Ising model) affects the coherent dynamics of the system.

Model
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 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): 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 section 2.3 [equation (4)]. Here, we renormalize each element of the coupling matrix as c jk = c jk /g, and to ease the notation, we use the symbols {X j (τ )} of the amplitudes, such that x j (t, τ ) = X j (τ )e iω 0 t + X * j (τ )e −iω 0 t . The equations for the complex amplitudes are then determined (j, k = 1, 2, 3): As in section 2, we decompose the coupling matrix c jk in equations (5)- (7) in terms of an antisymmetric and symmetric part, respectively r jk and α jk . Physically, the antisymmetric part corresponds to the energy-preserving coupling, and the symmetric part corresponds to the dissipative one. 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: α 12 = η and α 13 = α 23 = α, so that the coupling matrix in equation (5) reads   ẋ In the rest of the paper, our goal is to study the physics of the system near threshold as a function of the coupling parametersr,α, andη, where as before all quadratures are real (X 1,I = X 2,I = X 3,I = 0). Before discussing this general case, we first focus on two main configurations of interest for the coupling matrix in equation (8). Namely, assumingr,α ≥ 0: (i) the non-frustrated case, forη =α, and (ii) the fully-frustrated case, forη = −α. 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.

Non-frustrated network
First, we analytically discuss the threshold properties of the non-frustrated network, forη =α in equation (8). The most efficient eigenvalue λ max can be found by explicit inspection of the Jacobian matrix J(X 1 , X 2 , X 3 ) around the origin X 1 = X 2 = X 3 = 0. We analyze separately the cases ofr >α andr <α.
Since λ max is complex, beats are found. From the imaginary parts of λ max , one can find the expression of the beat frequency: In particular, forr approachingα, the frequency of the beats reduces towards zero with the critical behaviour of ω B ∼ (r −α) 1/3 , differently from the critical exponent of 1/2 in the two-oscillator case [34]. Whenr <α, we find that λ max = −1/2 +h/4 + F (α, −r) + G(α, −r). 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 forr = 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).
Forr <α, instead, λ max reads as λ max = −1/2+h/4+e i π/3 F (α,r)+e −i π/3 G(α,r). Therefore, the parametric oscillation occurs with beats, where the beat frequency is One can see by inspection that, whenr <α, the presence of the limit cycle at threshold makes the system periodically flip between six possible phase-locked configurations, The experimental simulation agrees exceptionally well with the predicted theoretical behaviour, apart from small deviations (some noisy regions and small discrepancies near the phase boundaries) that we ascribe mostly to the difficulty of estimating the oscillation threshold in the experimental simulation for some values of the system parameters.
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 whenr = 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 forr <α, 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 forr <α, and converges to phase-locked oscillations otherwise. For fixedα, the frequency of the beats increases linearly from zero for smallr, and goes to zero asr approachesα with the same critical exponent 1/3 as before.

Interpolating case
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 Label size for LaTeXiT: 22 (000), (πππ) (π00), (π0π) (0ππ), (0π0) (00π), (ππ0) Figure 6. Ising energy levels relative to the ground-state energy as a function of η/α, from the classical Ising energy that the dissipative components of the coupling reflect: is the Ising spin variable corresponding to the phase solution (respectively "0" or "π") of x j . The three energy levels are shown in different colors for different states: blue for (000) and (πππ), red for (π00), (π0π), (0ππ), and (0π0), and green for (00π) and (ππ0). The Ising gap goes to zero at the fully-frustrated point η/α = −1, where the six states (π00), (π0π), (0ππ), (0π0), (000) and (πππ) become degenerate. 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α, as a function ofr andη, and discern regions in parameter space in which 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). 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 colors). 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 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 suggest that 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. We therefore conclude that our findings provide a first explicit example where the presence of an energy-preserving coupling (even if small) represents a source of error for CIMs.
In addition, we find that, while the three-oscillator systems correctly behaves as a 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 3 (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. A deeper analysis on these points requires further analysis of larger spin models, which is beyond the scope of the present manuscript, and it is left for future work.

Numerical simulation of the experimental implementation
To corroborate the previous results and confirm our analytical predictions, we conducted a direct numerical simulation of the dynamics of the field inside a parametric-oscillator cavity with three (or more) modes, aiming to emulate as closely as possible the dynamics of a future experimental setup. In such an experiment, we intend to couple between parametrically driven modes of a multi-mode radio-frequency cavity in a fully-controlled and tunable manner. The dynamical coupling will be controlled by a field-programmable gate array (FPGA). Full details of the planned actual implementation will be reported in future work. Our numerical approach is completely distinct from the analysis presented above, and makes no explicit mention (or use) of the coupled Matheiu's equations of motion. In what follows, we first discuss our numerical procedure, and then compare the results of the simulated experiment with the analytical results presented in the previous sections.

Numerical procedure
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 S = (S 1 , S 2 , S 3 ) 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 parametric amplifier, the signal inside the cavity consists of noise alone. A round trip inside the cavity is identified by the following steps: First, the pump field and the signal are injected into the parametric amplifier. 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 parametric amplifier, 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 c = 1 − T 2 c . The transmitted signal T c S 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 c c S and the reflected signal R c S are combined on another coupling device, again with transmission and reflection coefficients T c and R c . At this step, the reflected signal T c R c S + R c F is blocked (contributing to the overall cavity losses), and the transmitted signal R 2 c S + T c F is fed back into the cavity. Last, the output coupler with transmission and reflection coefficients T out = 1/2 and R out = 1 − T 2 out allows to couple out the field T out (R 2 c S + T c F) from the cavity and analyze it both it in time and frequency, and the remaining field R out (R 2 c S + T c F), fed with noise, is input together with the pump field into the parametric amplifier to start the next round trip.  Figure 7. Schematic of the simulated experiment. At each round trip, the field inside the cavity is fed with noise and injected together with the pump field into a parametric amplifier. After the parametric amplification, part of the signal is sent into the coupling mechanism that couples between the time slots, and injected back into the cavity. At the end of the round trip, part of the signal is extracted from the cavity and measured. The green backslashed boxes denote a coupler device, identified by reflection and transmission coefficients (R, T ), whose values are chosen differently depending on the simulation step [(R c , T c ) for mode coupling, or (R out , T out ) for output coupling and noise feeding].
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.

Numerical results
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 fullyfrustrated 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  (10) and (11)].
frequency (divergence of the period of the beats), eventually inducing phase-locking as h is increased above the oscillation threshold [33]. 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.

Conclusions
We analyzed the behaviour of three coupled degenerate parametric oscillators -the minimal case to study nontrivial coupling and connectivity effects. By extending our previous work on two coupled parametric oscillators, we modelled the system as three coupled Mathieu's equations, where the coupling between any two oscillators is comprised of both energy-preserving and dissipative components. We focused in particular on two main cases of connectivity, namely, for frustrated and non-frustrated dissipative coupling. Our theoretical predictions, obtained by linearizing the effective equations of motion, were confirmed by a direct numerical simulation in time of the dynamics inside a parametric oscillator cavity, as it would be implemented in an actual experiment. The good agreement between the results obtained by these two different approaches strengthens the fact that the coupled nonlinear Mathieu's equations capture the relevant dynamics of coupled parametric oscillators.
Our main finding was that frustration of the dissipative component of the coupling has a dramatic effect on the coherent dynamics of the system. While in the nonfrustrated case the system phase locks once the dissipative coupling exceeds the energy-preserving one, and behaves as a CIM, in qualitative agreement with the behaviour of two coupled parametric oscillators, the frustrated case shows a totally reversed behaviour. In particular, when the dissipative coupling dominates close to full frustration, the Ising gap is vanishingly small and any vanishingly small energy-preserving coupling induces coherent beats between the (quasi) degenerate Ising configurations. For large values of the frustration parameter and for large energypreserving coupling, the system phase locks. Interestingly, in this phase-locking region the system of three coupled oscillators does not obey the Ising description.
Our results provide an additional piece of evidence of the highly nontrivial dynamics in networks of coupled parametric oscillators, which is considerably richer than an Ising network of spins. In the view of using coupled parametric oscillators to simulate Ising models, our results hint that, in situations where the energy gap of the corresponding Ising model is very small or vanishes, the presence of even a small energy-preserving coupling between the oscillators will induce coherent beats, causing the CIM operation to fail. Alternatively, in the beating phase the network continuously and coherently samples a complete subspace of the Ising configuration space, indicating that one may be able to engineer the combination of energy-preserving and dissipative couplings to harness the beats to guide the convergence to the Ising ground state. Because of these intriguing implications, the theoretical and experimental investigation of largescale networks is now highly desirable in order to establish our results for larger sets of parametric oscillators. We are currently planning the experimental implementation of such a large-scale network, whose analysis will be reported in future work. Ultimately, in light of our findings and the correspondence between the phase-locked behavior and classical time-crystals, it will be important to understand how nontrivial connectivities affect the stability of classical many-body time crystals. times, t n = n τ rt , for n = 0, 1, 2, . . ., via the splitter matrix as We consider c 11 = c 22 ≡ c, where 0 ≤ c ≤ 1 represents the transmittance coefficient of the coupling. With this choice, equation (A.1) is equivalently recast as , and by rewriting c .

Appendix C. Jacobian matrix spectrum in the interpolating case
In this appendix, we report the expressions of the eigenvalues of the Jacobian matrix around the origin in the interpolating case discussed in section 3.4 (see also figure 5). By generalizing the functions in equation (9), one has G(r,α,η) = 1 2  η r 2 −α 2 + r 2 − 2α 2 +η 2 3 By finding the most efficient eigenvalue λ max (r,α,η) according to the usual condition (the eigenvalue with largest real part), the frequency of the beats at threshold reads ω B (r,α,η) = ω 0 g |Im[λ max (r,α,η)]|.

Appendix D. Additional details on the choice of the round-trip time
In this appendix, we provide some details on the choice of τ rt discussed in figure 8. In figure D1, we show the comparison between the theoretical expressions of the frequency of the beats at threshold [equations (10) and (11)] ω B /α, as a function of r/α, and the data from the simulated experiment for different values of the round-trip time τ rt , as in the legends. The effect of changing τ rt is to renormalize the frequency units for the simulated experiment. Because of the excellent agreement between theory and data from the simulated experiment in the non-frustrated case for τ rt = 0.0195, this value of τ rt was chosen to quantitatively match the theoretical phase diagram in figure 5. Indeed, as evident, for this τ rt , theory and data are essentially overlapped.