Classical many-body time crystals

TL;DR

This paper introduces a simple classical framework for realizing many-body time crystals using coupled resonators, clarifies the distinction between single-mode and many-body symmetry breaking, and demonstrates the concept experimentally with mechanical oscillators.

Contribution

It provides a pedagogical classical approach to many-body time crystals and experimentally demonstrates their realization with coupled mechanical oscillators.

Findings

Classical period-doubling bifurcation theory explains time crystal formation.

Experimental demonstration with coupled mechanical oscillators.

Clear distinction between single-mode and many-body time-translation symmetry breaking.

Abstract

Discrete time crystals are a many-body state of matter where the extensive system's dynamics are slower than the forces acting on it. Nowadays, there is a growing debate regarding the specific properties required to demonstrate such a many-body state, alongside several experimental realizations. In this work, we provide a simple and pedagogical framework by which to obtain many-body time crystals using parametrically coupled resonators. In our analysis, we use classical period-doubling bifurcation theory and present a clear distinction between single-mode time-translation symmetry breaking and a situation where an extensive number of degrees of freedom undergo the transition. We experimentally demonstrate this paradigm using coupled mechanical oscillators, thus providing a clear route for time crystals realizations in real materials.

Equations

$${\ddot{x}}_i+{\omega }_i^2[1$$

$$+{\alpha }_i{x}_i^3+{\eta }_i{x}_i^2{\dot{x}}_i-\sum _{i\ne j}{\beta }_{ij}{x}_j=0,$$

$$\dot{\mathbf{X}}=A(\mathbf{X})\mathbf{X},$$

$$A=\left(\begin{array}{cccc}{a}_1(\mathbf{X})& b_{12}& \cdots & b_{1N}\\ b_{1,2}& a_2(\mathbf{X})& \ddots & ⋮\\ ⋮& \ddots & \ddots & b_{(N-1)N}\\ b_{1N}& \cdots & b_{(N-1)N}& a_N(\mathbf{X})\end{array}\right),$$

$$a_i(\mathbf{X})=-\frac{1}{8\omega }\left(\begin{array}{cc}{a}_{i,1}& a_{i,2}\\ a_{i,3}& a_{i,4}\end{array}\right),b_{ij}=\left(\begin{array}{cc}0& \frac{{\beta }_{ij}}{2\omega }\\ -\frac{{\beta }_{ij}}{2\omega }& 0\end{array}\right),$$

$$\lambda_{\rm th}^{s/a}=\frac{4\omega}{\omega_{0}^{2}}\sqrt{\frac{\gamma^{2}}{4}+\left(\frac{\omega^{2}-\omega_{0}^{2}}{2\omega}+\left\{\begin{array}[]{lr}\frac{(N-1)}{N},&\text{s}\\ \frac{-1}{N},&\text{a}\end{array}\right\}\frac{\beta^{2}}{2\omega}\right)^{2}}.$$

$$\dot{\mathbf{X}}=A\mathbf{X},$$

$$A$$

$$a$$

$$b$$

$$A=Id_{N}\otimes a+\left(\begin{array}[]{cccc}0&1&\cdots&1\\ 1&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&1\\ 1&\cdots&1&0\end{array}\right)\otimes b=Id_{N}\otimes a+M_{N}\otimes b\,.$$

$$A{w}_m$$

$$=r_m\otimes (a+{\rho }_{m}b)s_m$$

$$={\sigma }_m{w}_m.$$

$$\displaystyle a+\rho_{m}b=\left(\begin{array}[]{cc}a_{1}&a_{2}-a_{3}\\ a_{2}+a_{3}&a_{1}\end{array}\right)\,,$$

$${\sigma }_{m,\pm }$$

$$s_{m,\pm }$$

$$w_{m,\pm }$$

$${\sigma }_{m,+}=a_1+\sqrt{a_2^2-a_3^2}$$

$${\lambda }_{th,m}=\frac{4\omega }{{\omega }_0^2}\sqrt{a{1}^2+a_3^2}$$

Full text

Classical many-body time crystals

Time crystals are readily obtained in the steady state of many-body classical systems that undergo period-doubling bifurcations.

Toni L. Heugel1,∗, Matthias Oscity1,2,∗, Alexander Eichler3, Oded Zilberberg2,

and R. Chitra1

1 Institute for Theoretical Physics, ETH Zürich, Wolfgang-Pauli-Straße 27, 8093 Zürich, Switzerland.

2 Fachhochschule Nordwestschweiz FHNW, Klosterzelgstrasse 2, CH-5210 Windisch, Switzerland.

3 Institute for Solid State Physics, ETH Zürich, Wolfgang-Pauli-Straße 27, 8093 Zürich, Switzerland.

∗ These authors contributed equally.

Discrete time crystals are a many-body state of matter where the extensive system’s dynamics are slower than the forces acting on it. Nowadays, there is a growing debate regarding the specific properties required to demonstrate such a many-body state, alongside several experimental realizations. In this work, we provide a simple and pedagogical framework by which to obtain many-body time crystals using parametrically coupled resonators. In our analysis, we use classical period-doubling bifurcation theory and present a clear distinction between single-mode time-translation symmetry breaking and a situation where an extensive number of degrees of freedom undergo the transition. We experimentally demonstrate this paradigm using coupled mechanical oscillators, thus providing a clear route for time crystals realizations in real materials.

In periodically modulated nonlinear systems, discrete time-translation symmetry can be spontaneously broken, leading to inherently slower dynamics than that of the drive  (?, ?, ?, ?, ?, ?). A rapidly expanding community is principally focused on such a phenomenon in periodically-driven closed quantum systems, where disorder and interactions are considered to be essential for so-called discrete time crystals  (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?). A time-crystalline phase of matter stabilized by many-body localization was first observed in a one-dimensional trapped-ion system  (?). Surprisingly, time crystals were also seen in three-dimensional ensembles of NV-centers  (?) and in spin-$1/2$ nuclei in phosphate materials  (?) where disorder-induced localization effects are absent. The latter results indicate a wider class of time-crystalline behavior, including classical counterparts  (?).

A natural arena for realizing time crystals is provided by parametric resonators. A parametrically-pumped resonator mode plays an important role in many areas of science and technology. In its best-known form, parametric pumping describes the modulation of a resonator’s potential at twice its natural frequency  (?, ?, ?, ?). When the modulation depth exceeds an instability threshold, the resonator undergoes a period-doubling bifurcation to a new regime stabilized by nonlinearities  (?). This time-translation symmetry breaking (TTSB) leads to two stable parametric phase states that have equal amplitude, opposite phase, and half the oscillation frequency of the parametric drive  (?, ?, ?, ?). Interestingly, these states can be associated with two states of a classical bit  (?, ?, ?, ?, ?) or with an Ising spin  (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?). Network of such coupled resonators have been proposed as simulation platforms for complex Ising-like models that are very hard to solve with conventional computers  (?, ?, ?, ?, ?, ?, ?, ?, ?).

In this work, we show that a many-body TTSB can be easily realized in a classical network of dissipative parametric resonators. We present a general theoretical analysis and derive conditions for the manifestation of many-body TTSB in this system. This is complemented by a simple tabletop experimental demonstration using two coupled resonators. Our setup allows us to tune the coupling strength, and we find a regime where the modes of the system jointly undergo TTSB into well-defined parametric phase state configurations. Our experiment thus realizes the simplest building block that highlights the plethora of accessible TTSB solutions. At the same time, we test our understanding of the general many-body model against a well-controlled and accessible experimental implementation. Our work lifts the ambiguity surrounding the concept of time crystals by establishing sufficient conditions for their generation.

We consider a classical network of $N$ coupled nonlinear parametric oscillators, whose dynamics is governed by $N$ equations of motion

[TABLE]

where dots mark differentiations with respect to time $t$, $x_{i}$ is the displacement, $\omega_{i}$ is the eigenfrequency, $\gamma_{i}$ the dissipation, $\alpha_{i}$ the quartic nonlinearity, and $\eta_{i}$ the nonlinear damping of the $i^{\rm th}$ mode. The system is excited by a single parametric pump of modulation depth $\lambda$ and frequency $\Omega_{p}$. Each mode $i$ couples to other modes $j\neq i$ in the form of a driving force in proportion to $x_{i}$ and with a coupling coefficient $\beta_{ij}$.

We can perturbatively solve the system using the slow-flow method  (?): we rewrite Eq. (1) as $2N$ first-order differential equations and perform a van der Pol transformation with frequency $\omega=\Omega_{p}/2$, followed by time-averaging, to obtain the slow-flow equation

[TABLE]

where $\bm{X}=(u_{1},v_{1},u_{2},v_{2},\ldots,u_{N},v_{N})^{T}$, with $u_{i}$ and $v_{i}$ the slowly varying phase-space quadratures of the individual resonators. This equation is valid if the dimensionless quantities $1-(\frac{\omega}{\omega_{i}})^{2}$, $\lambda$, $\gamma_{i}/\omega_{i}$, $\frac{\eta_{i}}{\omega_{i}}x_{i}^{2}$, $\frac{\beta_{ij}}{\omega_{i}^{2}}$, and $\frac{\alpha_{i}}{\omega_{i}^{2}}x_{i}^{2}$ are of order $\epsilon$, where $0<\epsilon\ll 1$  (?). These conditions are easily satisfied for a network of nearly identical oscillators. The matrix $A$ can be written as

[TABLE]

where the $a_{i}$ and $b_{ij}$ are given by,

[TABLE]

with $i\neq j$ and $i,j=1,2,\ldots,N$, and using the definitions $a_{i,1}=a_{i,4}=4\gamma_{i}\omega+\eta_{i}\omega X_{i}^{2}$, $a_{i,2}=2\left(\lambda\omega_{i}^{2}+2(\omega_{i}^{2}-\omega^{2})\right)+3\alpha_{i}X_{i}^{2}$, $a_{i,3}=2\left(\lambda\omega_{i}^{2}-2(\omega_{i}^{2}-\omega^{2})\right)-3\alpha_{i}X_{i}^{2}$, and $X_{i}=u_{i}^{2}+v_{i}^{2}$. In general, the number of steady-state solutions, both stable and unstable, to this $N$-body problem varies from $1$ to $5^{N}$ depending on the parameter regimes  (?).

In the absence of nonlinearities, $\alpha_{i}=\eta_{i}=0$, the natural description of the resonator network is given by $N$ normal modes with eigenfrequencies $\nu_{k}$, $k=1,\ldots,N$. The dynamics of the normal modes is determined by the eigenvalues and the eigenvectors of $A$. The $2N$ eigenvectors define the positions and momenta of the $N$ normal modes. The time evolution of the $k$-th normal mode is given by $e^{\Lambda_{k}t}$, with $\Lambda_{k}$ the respective eigenvalue. The motion will be bounded for negative ${\rm Re}\{\Lambda_{k}\}$ and manifest parametric instability, i.e., unbounded dynamics when ${\rm Re}\{\Lambda_{k}\}>0$. Each normal mode exhibits a corresponding parametric stability phase diagram known as ‘Arnold tongues’, delineating regions where dissipation stabilizes the motion and regions where the linear system shows unbounded dynamics, see Fig. 1(a). In the following, we will focus on the dominant instability lobe occurring around twice the natural frequency of the normal mode $k$, $\Omega_{p}\sim 2\nu_{k}$, when the parametric drive exceeds a threshold $\lambda\geq\lambda_{\rm th}^{k}$  (?).

In general, it is not dissipation but the underlying nonlinearities $(\alpha_{i},\eta_{i})$ that stabilize the normal-mode oscillations against unbounded growth  (?). At the boundary of its main instability lobe, each normal mode undergoes a period-doubling bifurcation alongside a spontaneous $Z_{2}$ symmetry breaking between the two parametric phase states, see Fig. 1(b). This is a simple manifestation of TTSB in the steady state of an effective single parametric mode. The parametric phase states define an effective Ising-like phase bit  (?, ?, ?, ?, ?). It is important to note that although a single normal mode can involve an extensive number of resonators of the network, it does not give rise to a many-body TTSB because it does not involve an extensive number of independent degrees of freedom.

A many-body TTSB phase is realized in the resonator network in a region where an extensive number of normal modes undergo the aforementioned period-doubling transition. A simple recipe to realize a many-body TTSB consists of finding the parametric pumping amplitude $\lambda_{\rm th}^{\rm MB}(\omega_{P})=\min_{\lambda}\{\lambda>\lambda_{\rm th}^{k}(\omega_{P})\,,\,\forall k\}$ at which all normal modes are driven above their respective instability thresholds, see Fig. 1(c). There, each normal mode finds itself in a parametric phase state, see Fig. 1(d). Note that the many-body threshold holds in the limit of weak nonlinearities and does not include corrections stemming from nonlinear inter-normal mode coupling. In the mean-field limit of $N$ identical resonators, i.e., $\omega_{i}\equiv\omega_{0}$ and $\gamma_{i}=\gamma$, with all-to-all coupling $\beta_{ij}=\beta/\sqrt{N}$, apart from the symmetric mode, all other instability lobes coincide with that of the antisymmetric ($a$) mode. The respective instability thresholds ($\lambda>\lambda_{\rm th})$ are given by  (?):

[TABLE]

The overlap region of $\lambda_{\rm th}^{s/a}$ defines $\lambda\geq\lambda_{\rm th}^{\rm MB}(\Omega_{p})$.

In the following we discuss two limits, ‘strong’ and ‘weak’ coupling, that are defined relative to the parametric modulation strength, $\lambda$. For weak $\beta_{ij}$, the normal modes closely resemble the underlying constituent resonators. However, as $\beta_{ij}$ increase, the normal modes become collective in nature. In the many-body TTSB phase the system can choose one of $2^{N}$ to $3^{N}$ configurations: in the weak coupling regime, these correspond to the possible configurations of the $N$ individual resonators  (?, ?), while in the strong coupling regime, they correspond to the configurations of collective normal modes. In both cases, all these configurations manifest TTSB and the chosen configuration will depend on initial conditions, noise and the strength of the nonlinearities. To summarize, we predict that an array of coupled dissipative parametric resonators realizes a stable TTSB phase in its steady state. This phase endures in a wide region of parameter space and is robust to fluctuations.

We now report on an experimental demonstration of many-body TTSB in a system of two coupled mechanical modes. Our setup is based on the lowest transverse vibrational modes of two macroscopic strings. The strings are clamped onto a fixed frame at one end, while the other end is attached to a stiff plate that has two purposes; firstly, the plate can be driven into vibrations parallel to the string axes by an electric motor. These vibrations modulate the tension inside the strings and generate parametric pumping of both string modes. Secondly, the plate transmits vibrations between the strings, which leads to weak intrinsic coupling between the modes. In some experiments, we introduce strong mode coupling by way of a mechanical connection close to the mode antinodes, see Fig. 2(a).

The motion of each string is independently measured with a dedicated piezo detector embedded into one clamping point. We use a lock-in amplifier (Zurich Instruments HF2 LI) to actuate plate vibrations and to read out the electrical signals from the two piezo detectors, which are proportional to the strings’ displacements. All measurements in this work were carried out in the form of frequency sweeps, where the actuation frequency $\Omega_{p}=2\omega$ and the detection frequency $\omega$ were swept slowly to capture the steady-state response of the modes.

We use weak external driving for calibration of the modes, similarly to the procedure outlined in Ref.  (?, ?). In these experiments, the vibration amplitude is kept low and the influence of the intrinsic coupling is negligible. From the Lorentzian response of each mode, we extract typical values for $\omega_{1,2}/2\pi=155\pm 10$ Hz (depending on ambient temperature) and $Q_{1,2}\sim 1200$, while we calculate the effective mass $m=4.14\times 10^{-4}$ kg from the geometry of the strings. By fitting to the large-amplitude response under strong parametric pumping, we obtain the coefficients of the nonlinear potential term, $\alpha_{1}=11.93$ mV${}_{\rm{RMS}}^{-2}$s*-2* and $\alpha_{2}=6.24$ mV${}_{\rm{RMS}}^{-2}$s*-2*, as well as those of the nonlinear damping, $\eta_{1}=7.1\,\mu$V${}_{\rm{RMS}}^{-2}$s*-1* and $\eta_{2}=3.9\,\mu$V${}_{\rm{RMS}}^{-2}$s*-1* (in the strong coupling case, we find $\eta_{1}=3.55\,\mu$V${}_{\rm{RMS}}^{-2}$s*-1* and $\eta_{2}=1.95\,\mu$V${}_{\rm{RMS}}^{-2}$s*-1*)  (?). Finally, in the presence of strong coupling, we use the normal mode frequency splitting to estimate $\beta=36.2\pm 0.1$ Hz.

Strong coupling [Fig. 2(a)]: we first explore the regime where the two instability lobes corresponding to the symmetric and antisymmetric normal modes are well separated, see Fig. 2(b). In Figs. 2(c) and (d), we show the measured amplitudes and phases of both strings under a common parametric modulation as a function of frequency $\omega/2\pi$, respectively. As the frequency is slowly swept upwards, both resonators oscillate with the same phase from $146$ Hz up to $149$ Hz. As the frequency is ramped further, the resonators are in opposing phase states from $155$ Hz up to $157.5$ Hz. The modes exhibit identical symmetries (s/a) when the frequency is swept downwards. These qualitative observations were consistent over many sweeps. The small peaks around $\omega/2\pi=153$ Hz correspond to an unidentified eigenmode in the experimental setup that does not appear to affect the modes of interest.

We model the system with Eq. (2) for $N=2$ using the parameters extracted from the experiment. The results of our calculations provide a simple understanding of the experimental observations: as the frequency is swept, either the symmetric or antisymmetric normal modes undergo TTSB at their respective instability thresholds, recreating the effective single-mode TTSB discussed earlier. The coupling between the normal modes induced by nonlinearities is irrelevant in this regime as one mode is strongly off-resonant with the other. The experimental results are well described by the phase-space bifurcation diagrams for each resonator plotted in Figs. 2(e) and (f). Despite the fact that both resonators participate in the TTSB of the symmetric or antisymmetric modes, many-body TTSB is not observed in this strong-coupling limit as the two instability lobes do not overlap for experimentally accessible parametric excitation strengths.

Weak coupling: next, we remove the connection between the strings and rely on the driving plate to provide weak coupling between the string modes [Fig. 3(a)-(b)]. The experimental data look very different in this regime [Fig. 3(c)-(d)]. Both strings have nearly identical natural frequencies (within $50$ mHz from each other) and exhibit hysteresis when sweeping the frequency upwards and downwards. The frequencies where the oscillation drops to zero (during upsweeps) or jumps to a finite amplitude (during downsweeps) are precisely the same for both resonators. The strings oscillate in phase during the upsweep and out of phase during the downsweep. All of these features were reproduced over many sweeps.

The theoretical model corresponds to normal modes that are split by a very small coupling $\beta$, such that their instability lobes overlap strongly [Fig. 3(b)]. Since both normal modes exhibit TTSB and are weakly coupled by nonlinearities, we witness the realization of two-body TTSB. As before, the experimental results for the amplitude and phase are consistently explained by the weak coupling bifurcation diagram for both strings shown in Figs. 3(e)-(f). In comparison with the strong coupling scenario of Figs. 2(e)-(f), the weakly coupled system exhibits richer behavior. The selection of symmetric and antisymmetric solutions as a function of the sweeping direction may be explained in terms of the phase response of a linear resonator to a periodic external force. Below its natural frequency, a harmonic resonator oscillates with almost no phase lag in response to an external force. As the two string modes drive each other, they prefer to move in phase. In contrast, since the harmonic resonator response has a phase lag of $\sim\pi$ above the natural frequency, the string modes preferably oscillate out of phase during the downsweep. This many-body TTSB state is stable against small detunings $\omega_{1}\neq\omega_{2}$ and robust to noise (as seen in the experiment). Increasing noise levels are expected to preserve the underlying TTSB, but to induce transitions between the different stable solutions.

Coupled parametric resonators provide the simplest platform to realize macroscopic states with robust discrete time-translation symmetry breaking. Period-doubling bifurcations in stable steady-states provide a rich space of solutions that manifest such many-body phenomena. This can be readily generalized to the quantum realm, where the bifurcations physics is replaced by dissipative first- and second-order phase transitions  (?, ?, ?, ?). Furthermore, higher-period TTSB can also be realized in these systems through a judicious choice of modulated nonlinearities  (?). In the weak coupling limit, the classical network can be viewed as an Ising machine that simulates complex problems, where the system parameters can be tuned to engineer desired ‘spin configurations’ of the Ising-like phase states. The analogous quantum network comprising dissipative Kerr parametric resonators is expected to manifest an equivalent TTSB phase  (?, ?, ?, ?, ?, ?). Such networks have been proposed as quantum annealers  (?), and following this work can now be used as quantum simulators of many-body time crystals.

I Derivation of $\lambda_{th}$

We present here the derivation of the many-body parametric driving threshold amplitude for $N$ resonators that are equally coupled to one another. For the prupose of this calculation, it suffices to consider $N$ coupled linear resonators. The slow-flow equation describing this system (cf. Eq.(2) in the paper with coupling $\beta_{ij}=\beta$ for all $i\neq j$) is given by:

[TABLE]

where the matrix $A$ is given by:

[TABLE]

and the individual matrix entries $a$ and $b$ are given by:

[TABLE]

The dynamics of the linear system can be deduced by decomposing the initial state $\bm{X}$ into the eigenvectors of $A$. The time evolution of each eigenvector is then determined by $e^{\Lambda t}$, where $\Lambda$ is the corresponding eigenvalue. ${\rm Im}\Lambda\neq 0$ imposes an oscillatory behavior whose envelope decreases exponentially for ${\rm Re}\Lambda<0$ and increases exponentially for ${\rm Re}\Lambda>0$. To evaluate these eigenvalues and eigenvectors, it is useful to rewrite the matrix $A$ as:

[TABLE]

Based on the structure of $A$, the eigenvectors obey the ansatz $w_{m}=r_{m}\otimes s_{m}$, where $r_{m}$ are the $N$ dimensional eigenvectors of $M_{N}$ with eigenvalue $\rho_{m}$ and $s_{m}$ are the 2-dimensional eigenvectors of $a+\rho_{m}b$ with eigenvalues $\sigma_{m}$. Since $M_{N}$ has $N$ eigenvectors, $r_{m}$, with $2$ corresponding $s_{m}$, this ansatz describes all the $2N$ eigenvectors and eigenvalues of the matrix $A$. The 2-vector $s_{m}$ describes the amplitude and momentum $(u_{m},v_{m})$ in a particular mode’s phase-space and $r_{m}$ generically describes the relative amplitudes and the phase configuration, e.g., $r_{m}=(1,-1)$ means that the two oscillators have opposite phases. We can readily show that this ansatz is indeed an eigenvector of $A$:

[TABLE]

Since $M_{N}$ has a simple structure, we see that the eigenvectors $\{r_{m}\}$ take the form $r_{0}=(1,1,\cdots)$ with eigenvalue $\rho_{0}=N-1$ and $r_{m}=(0,\cdots 0,1,-1,0\cdots)^{T}$, where the $+1$ is the $m^{th}$ entry, are eigenvectors of $M_{N}$ with eigenvalues $\rho_{m}=-1$ ($m\in\mathbb{N}$, $0<m\leq N-1$). Note that the eigenvectors $r_{m}$ effectively determine the normal mode transformations of the problem. Next, we evaluate the eigenvectors $s_{m}$ and eigenvalues $\sigma_{m}$ of

[TABLE]

where $a_{1}=-\frac{\gamma}{2}$, $a_{2}=-\frac{\lambda\omega_{0}^{2}}{4\omega}$ and $a_{3}=\frac{\left(\omega_{0}^{2}-\omega^{2}\right)}{2\omega}-\rho_{m}\frac{\beta}{2\omega}$. These are given by,

[TABLE]

To summarize, the $2N$ eigenvectors of the matrix $A$ are given by

[TABLE]

with corresponding eigenvalues $\sigma_{m,\pm}$ and $0\leq m\leq N-1$.

If ${\rm Re}\sigma_{m,\pm}>0$, the corresponding $w_{\pm,m}$ grows exponentially indicating a parametric instability. We obtain the parametric driving threshold $\lambda_{th,m}$ for this instability by imposing the condition:

[TABLE]

Note that we have $a_{1}<0$, whereas $\sqrt{a_{2}^{2}-a_{3}^{2}}$ can be either real-valued and positive or complex-valued. Solving Eq. SI.19, we obtain

[TABLE]

For identical oscillators, we see that there are primarily two instability thresholds corresponding to (i) the instability of the symmetric normal mode, $w_{0,+}$, and (ii) to the instability of all other normal modes: $w_{m,+}$ including the antisymmetric mode.

II Calibration measurements

In Fig. S1, we present test measurements that we have performed to ensure that the weakly coupled strings were degenerate in frequency. On timescales of hours, thermal drift sometimes caused detuning between the strings, which we balanced by adjusting the tension of the strings separately. In Fig. S2, we show the fits used to extract the nonlinear coefficients of the two weakly coupled strings. Please refer to Ref. [23] of the main text for details regarding the model of a nonlinear parametric oscillator.