Higgs time crystal in a high-$T_c$ superconductor

We propose to induce a time crystalline state in a high-$T_c$ superconductor, by optically driving a sum resonance of the Higgs mode and a Josephson plasma mode. The generic cubic process that couples these fundamental excitations converts driving of the sum resonance into simultaneous resonant driving of both modes, resulting in an incommensurate subharmonic motion. We use a numerical implementation of a semi-classical driven-dissipative lattice gauge theory on a three-dimensional layered lattice, which models the geometry of cuprate superconductors, to demonstrate the robustness of this motion against thermal fluctuations. We demonstrate this light-induced time crystalline phase for mono- and bilayer systems and show that this order can be detected for pulsed driving under realistic technological conditions.


I. INTRODUCTION
Optical driving of solids constitutes a new method of designing many-body states. Striking examples of this approach include light-induced superconductivity [1][2][3] as well as optical control of charge density wave phases [4]. For these states, the carefully tuned light field either renormalises the phase boundary of the equilibrium phase, as is the case for light-induced superconductivity, or renormalises a near-by metastable state into a stable state of the driven system, as is the case for lightcontrolled charge density waves.
These observations are part of a larger effort to determine the steady states of periodically driven manybody systems. In a parallel development in cold atom systems, serving as well-defined many-body toy models, the generic regimes that were proposed, see Refs. [5,6], firstly include renormalised equilibrium states, for which the above mentioned states are examples. Secondly, regimes beyond the equilibrium states emerge, in particular genuine non-equilibrium orders, which have no equilibrium counterpart, and only exist in the driven state. A striking example of a non-equilibrium order is time crystals [7][8][9][10][11][12][13], reported in systems such as ion traps or nitrogen-vacancy centers [14,15]. Thirdly, for strong driving, chaotic states emerge. These different regimes are achieved for different driving amplitudes and driving frequencies, which constitutes the dynamical phase diagram of the system.
In this paper, we propose to create a light-induced time crystalline state in a high-T c superconductor. This advances light control of superconductors towards genuine non-equilibrium orders, and furthers time crystals in the solid state domain [16]. We characterise the observed non-equilibrium state as a time crystal based on the following criteria [12]: (i) A time crystal spontaneously breaks time-translation symmetry, that is, it exhibits a * ghomann@physnet.uni-hamburg.de subharmonic response to the drive. (ii) The subharmonic response is robust against perturbations which respect the time-translation symmetry of the Hamiltonian. (iii) The subharmonic response emerges in a many-body system with a large number of locally coupled degrees of freedom, and it persists for an infinite time.
We call the novel dynamical phase a Higgs time crystal because we induce it via optical driving of a sum resonance of the Higgs mode and a Josephson plasma mode. The Higgs mode and the Josephson plasma mode correspond to the two fundamental collective excitations of a system with broken U (1) symmetry and with an underlying approximate particle-hole symmetry. The Higgs mode is an amplitude oscillation of the order parameter, as depicted in Fig. 1(a) for the |ψ| 4 theory used in the following. The Higgs mode is a gapped excitation due to the increase of the potential energy in the radial direction. The Josephson plasma mode is a phase oscillation, as indicated. This mode also has a gapped excitation spectrum owing to the electromagnetic interaction of the system. Because of the approximate particle-hole symmetry, these two oscillations are orthogonal to each other [17,18].
To identify the Higgs time crystalline phase, we map out the dynamical phase diagram of optically driven high-T c superconductors as a function of the driving frequency ω dr and the driving amplitude E 0 , which is shown in Fig. 2(a), for instance. The time crystalline state is induced by driving the superconductor via the non-linear coupling ∼ a 2 h of the electromagnetic field a and the Higgs field h. We demonstrate that driving at the frequency ω dr = ω H + ω J induces a time crystalline phase, where ω H is the Higgs frequency and ω J is the plasma frequency, as depicted in Fig. 1(d). We note that this nonlinear coupling has been confirmed in conventional superconductors [19][20][21][22], while a direct probe of the Higgs field is challenging due to its scalar nature. Further studies on the Higgs mode in high-T c cuprates and organic superconductors are reported in Refs. [23][24][25][26][27][28][29]. Persistent multi-frequency dynamics of the superconducting order parameter has been investigated in Ref. [30]. To describe the dynamics of optically driven superconductors, we develop a lattice gauge simulation that describes the motion of the order parameter of the superconducting state ψ(r, t) and the electromagnetic field A(r, t). We first utilise our method to show how to induce the time crystalline state, and to determine its regime in the dynamical phase diagram. Furthermore, we demonstrate the robustness of the time crystalline phase against thermal fluctuations, and show that it can be realised and identified under pulsed operation.

II. THREE-DIMENSIONAL LATTICE GAUGE MODEL
We represent the layered structure of high-T c superconductors via the lattice geometry illustrated in Fig. 1(b). We note that this geometry of CuO 2 layers perpendicular to the c-axis has motivated a low-energy description of stacks of Josephson junctions [31][32][33], which captures the appearance of Josephson plasma excitations reported in Refs. [34][35][36]. Each layer is represented by a square lattice, leading to a discretisation of the fields of the form ψ(r, t) → ψ l,m,n (t) ≡ ψ r (t). The in-plane discretisation length d ab constitutes a short-range cut-off well below the in-plane coherence length. In doing so, we generalise the modelling of layered cuprates to a three-dimensional (3D) lattice of Josephson junctions. Each component of the vector potential A i,r (t) is located between a lattice site r and its nearest neighbour r (i) in the i-direction, where i ∈ {x, y, z}. According to the Peierls substitution, it describes the averaged electric field along the bond of a plaquette in Fig. 1 We focus on temperatures below T c , where the dominant low-energy degrees of freedom are Cooper pairs. We describe the Cooper pairs as a condensate of interacting bosons with charge −2e, represented by the complex field ψ r (t). To construct the Hamiltonian of the lattice gauge model, we discretise the Ginzburg-Landau free energy [37] on a layered lattice and add time-dependent terms. We explicitly simulate the coupled dynamics of the condensate and the electromagnetic field. We discretise space by mapping it on a lattice, as mentioned, but implement the compact U (1) lattice gauge theory in the time continuum limit [38]. The particle-hole symmetry inherent to our relativistic model creates stable Higgs oscillations, even in bilayer cuprates where the Higgs frequency is between the two longitudinal Josephson plasma frequencies.
We consider mono-and bilayer cuprate superconductors. For bilayer cuprates, we assign the strong (weak) junctions to the even (odd) layers. The corresponding tunnelling coefficients are t 2n = t s and t 2n+1 = t w . The interlayer spacings d 2n,2n+1 = d s,w are the distances between the CuO 2 planes in the crystal. Note that we suppose the z-direction to be aligned with the c-axis of the crystal throughout this paper. The Hamiltonian of the lattice gauge model is The first term is the |ψ| 4 model of the superconducting condensate in the absence of Cooper pair tunnelling: where π r = K 2 ∂ t ψ * r is the conjugate momentum of ψ r , µ is the chemical potential, and g is the interaction strength. This Hamiltonian is particle-hole symmetric due to its invariance under ψ r → ψ * r . The coefficient K describes the magnitude of the dynamical term.
The electromagnetic part H em is the discretised form of the free field Hamiltonian, modified by tunable interlayer permittivities s,w to capture the screening due to bound charges in the material: where E i,r denotes the i-component of the electric field. The vector potential is located on the bonds between the superconducting sites. Consequently, this applies to the electric field as well. Note that we choose the temporal gauge for our calculations, i.e., E i,r = −∂ t A i,r . Meanwhile, the magnetic field components B i,r = ijk δ j A k,r are centred about the plaquettes. This arrangement is consistent with the finite-difference time-domain (FDTD) method for solving Maxwell's equations [39]. We calculate the spatial derivatives according to δ j A k,r = (A k,r (j) − A k,r )/d j,r , where d j,r is the length of the bond. The dielectric permittivities are x,r = y,r = 1 and z,r = n . The other prefactors in Eq. (3) account for the anisotropic lattice geometry. They are defined as κ x,r = κ y,r = 1 and κ z,r = d n /d c , while β x,r = β y,r = 2ed ab d n / and β z, The non-linear coupling between the Higgs field and the electromagnetic field derives from the tunnelling term The unitless vector potential a i,r = −2ed i,r A i,r / couples to the phase of the superconducting field, ensuring the gauge-invariance of H kin . The in-plane tunnelling coefficient is t ab , and the c-axis tunnelling coefficients are t s,w . We solve the equations of motion for ψ r (t) and A r (t) obtained from the Hamiltonian numerically, employing Heun's method with an integration step size ∆t = 1.6 as. Thermal fluctuations are included by adding dissipation and Langevin noise to the equations of motion for both fields. For example, the time evolution of the superconducting field is given by where γ is a damping constant and ξ r represents white Gaussian noise with zero mean, see Table I for noise correlations. We note that the inclusion of in-plane dynamics and arbitrarily strong amplitude fluctuations constitutes a qualitative advance of previous descriptions, such as 1D sine-Gordon models [40,41]. We determine the response of the superconductor to periodic driving of the electric field along the c-axis. The external drive E dr (t) has the frequency ω dr and the effective field strength E 0 . We consider the long-wavelength limit such that the external drive is assumed to be homogeneous in the bulk of the sample. Thus, the time evolution of E z,r (t) reads where η z,r is white Gaussian noise with zero mean. The equations of motion for E x,r (t) and E y,r (t) are analogous to Eq. (6), except for the driving term. We characterise the response by evaluating the sample averages of the condensate amplitude |ψ(t)| and the supercurrent density J(t), see also Appendix B.
By applying the optical driving as described, we obtain the full dynamical phase diagram due to the direct coupling of the electromagnetic field to the superconducting order parameter. We note that resonant optical driving of phonon modes has been utilised and discussed in Refs. [1][2][3][40][41][42]. Here, we ignore the phononic resonances, so that our predictions are valid away from these resonances. A combined description will be given elsewhere.

III. TWO-MODE MODEL
Before we present the full numerical simulation, we identify the main resonant phenomena of the system. We consider the zero-temperature limit, where the in-plane dynamics can be neglected and the model simplifies to a 1D chain along the c-axis. Furthermore, we restrict ourselves to weak driving and a monolayer structure with t s = t w ≡ t J and d s = d w ≡ d. For periodic boundary conditions, the time evolution then reduces to two coupled equations of motion. Keeping only linear terms except for the lowest order coupling between the Higgs field and the unitless vector potential, we find where h = (ψ − ψ 0 )/ψ 0 is the Higgs field with ψ 0 being the equilibrium condensate amplitude, γ is the damping constant, and α is the capacitive coupling constant of the junction. Note that the unitless vector potential a equals the phase difference between adjacent planes in this setting. The external drive appears through the current j dr . The Higgs and plasma frequencies are ω H = 2µ/K 2 and ω J = t J /αK 2 , respectively. The main finding of this work is the emergence of a time crystalline phase by driving at the sum of the Higgs and plasma frequencies, ω dr = ω J + ω H . A cubic interaction process, visualised in Fig. 1(d), allows for simultaneous resonant driving of both the Higgs and the plasma modes [43].
In addition to the sum resonance, we identify various other resonances from the simplified equations of motion. For a response of the vector potential at the driving frequency, i.e., a = a 1 cos(ω dr t), Eq. (8)  oscillator with a resonance at ω dr = ω H /2. This recovers the sub-gap Higgs resonance [22]. The sub-gap resonance and the sum resonance originate from the same cubic coupling term ∼ a 2 h, as illustrated in Figs. 1(c) and 1(d). Next, we consider the range of driving frequencies where the Higgs field exhibits a second-harmonic response, that is, the external drive induces Higgs oscillations of the form h = h 0 + h 1 cos(2ω dr t) through the a 2 term in Eq. (8). For small driving amplitudes, the ah term in Eq. (7) can be neglected so that the equation reduces to a forced oscillator with a resonance at ω dr = ω J . However, the response is modified once the coupling to the Higgs field becomes significant. Then, Eq. (7) approaches a parametrically driven oscillator. The parametric resonances emerge at ω dr = ω J /k, where k ∈ N.

IV. DYNAMICAL PHASE DIAGRAM
We now present our numerical results in two steps. Firstly, we verify our analytical predictions for the resonances and, in particular, the Higgs time crystal by mapping out the dynamical phase diagrams of mono-and bilayer cuprate superconductors at zero temperature. We will show how the sum resonance is modified in a bilayer system, which has two plasma modes. Secondly, we test the robustness of this phase against thermal fluctuations using finite-temperature simulations.

A. Monolayer cuprate superconductor
Here, we consider a monolayer cuprate with ω H /2π ≈ 6.3 THz, ω J /2π ≈ 16.0 THz, γ/2π = 0.5 THz, and α = 0.33, see Table I for full parameter set. The sys-tem is continuously driven at various amplitudes and frequencies in the terahertz regime. In each realisation, the drive is applied for 20 ps and the relevant frequency spectra are computed using the final 10 ps, which amounts to 5 < M tot < 300 driving cycles in the frequency range of interest. The dynamical phase diagram in Fig. 2(a) is mapped out by analysing the normalised power spectra of |ψ(t)| and J(t) defined as , and T s = 10 ps is the sampling interval. Specifically, we obtain the spectral entropy for the dynamics of the condensate amplitude, S |ψ| = − dωP |ψ| (ω)lnP |ψ| (ω).
The heating regime, which is characterised by a strong depletion of the condensate, is identified based on the threshold S |ψ| > 2.2 × 10 −2 . It indicates the appearance of resonant phases associated with the Higgs and plasma excitations. We note that the two dominant heating tongues are weakly red-detuned from the expected resonance frequencies ω H /2 and ω J , respectively. Such a renormalisation of the fundamental frequencies is inherent to strongly driven non-linear systems [44]. This effect is further amplified by the damping terms present in our model. We identify the small tongue at ω dr /2π ≈ 4.8 THz as the third order parametric resonance of the Josephson plasma mode around ω J /3.
For intermediate driving intensity, we observe several dynamical regimes due to resonances. The resonance with the lowest frequency is the Higgs resonance at ω dr = ω H /2. In general, resonant excitation of the Higgs mode is marked by strong modulation of the condensate amplitude as exemplified in Fig. 2(b). Moreover, the Higgs resonance exhibits a commensurate and superharmonic response of |ψ(t)| with respect to the driving E dr (t) as seen from the closed trajectory in Fig. 1(c) and the sharp peak at 2ω dr of the condensate amplitude spec-trum in Fig. 2(c). We emphasise that driving away from any noticeable resonance, indicated as the blue regime in Fig. 2(a), induces only a single sharp peak in the supercurrent spectrum, namely at the driving frequency. The condensate amplitude oscillates at twice the driving frequency in the blue regime. This also applies to the regime near the Josephson plasma resonance at ω dr = ω J , where the system responds with strong oscillations of the supercurrent.
The red regime in Fig. 2(a), identified via the condition 10 −4 < S |ψ| < 2.2 × 10 −2 , is the Higgs time crystal introduced earlier. We emphasise that its resonance condition ω dr = ω J + ω H differs from the sub-gap frequencies ω dr ω H /2 used in standard Higgs spectroscopy. The sum resonance simultaneously couples to the Higgs and plasma resonances as evident from the exemplary meanfield trajectory in Fig. 1(d), where the amplitude oscillation is accompanied by a strong oscillation of the phase difference between the junctions. Despite a smaller driving amplitude E 0 , the plasma mode is excited with larger amplitude than for the Higgs resonance. The strong activation of the plasma mode results in a partial depletion of the condensate as visible in Fig. 2(b), where the time average of the oscillatory motion of the condensate amplitude is below 1. The key feature of the novel phase is the subharmonic response of the condensate amplitude as |ψ(t)| oscillates at ω H when the superconductor is driven at ω dr = ω J +ω H . This phenomenon is highlighted in Fig. 2(b) and in the strong subharmonic peak in the power spectrum of |ψ(t)| shown in Fig. 2(c). The other dynamical phases respect the time-translation symmetry imposed by the external drive as evidenced by Figs. 2(b) and 2(c).
The subharmonic collective motion is one of the defining features of a time crystal. In addition to being subharmonic, the response of the time crystalline state is also incommensurate to the external driving. That is, the phase-space trajectory traces an open loop for any number of driving cycles, see also Fig. 1(d). Therefore, and more specifically, the state that we propose to create is an incommensurate time crystal in high-T c superconductors. We will confirm its robustness against perturbations of the drive and thermal fluctuations for the bilayer case. We note that the subharmonic response can be expected to be rigid as it emerges for a broad regime of driving parameters rather than a fine-tuned point in the dynamical phase diagram. In addition, our finitetemperature calculations with a large number of lattice sites will highlight the many-body nature of the Higgs time crystal.

B. Bilayer superconductor
We now focus on bilayer cuprates. Due to the staggered tunneling coefficients t s and t w along the c-axis, the system has two fundamental longitudinal plasma excitations with frequencies ω J1 and ω J2 . The dynamical  Table I for full parameter set.
phase diagram at zero temperature in Fig. 3(a) displays a regime, in which a Higgs time crystal is induced by optical driving at a sum resonance. Here, the resonance condition is ω dr = ω H + ω J2 , so it is the sum of the Higgs frequency and the upper plasma frequency.
First, we examine how perturbing the optical drive itself affects the subharmonic response. To excite the sum resonance, we initially drive the bilayer superconductor with E 0 = 0.1 MV cm −1 and ω dr /2π = 21 THz. At some instant of time t 0 , the driving is altered so that the oscillation amplitude of the field strength depends on its sign for t > t 0 : for cos(ω dr t) ≥ 0, (E 0 + δE) cos(ω dr t) for cos(ω dr t) < 0.
After allowing the system to relax to a steady state, we take the power spectrum of the condensate amplitude and determine the dominant frequency ω peak . The robustness of the subharmonic response is demonstrated by Fig. 3(b), where perturbations of the driving amplitude between δE/E 0 = −0.4 and δE/E 0 = 1 do not destroy the sum resonance. We have also verified the persistence of the subharmonic response for 10 5 cycles of continuous driving at T = 0 [43]. Because of experimental and numerical limitations in accessible timescales (∼ 10 2 driving cycles for our finite-temperature calculations), we will not distinguish here between a 'true' time crystal and a slowly decaying time crystal [13,14]. We note that the time crystalline response is stabilised by the non-linear coupling between the Higgs and plasma modes, which further highlights the collective nature of the Higgs time crystal. Furthermore, the amplitudes of the oscillations are saturated by non-linear processes in (d) Temperature dependence of the optimal time crystalline fraction for a bilayer cuprate superconductor, rescaled by its value at T = 0. The optimal crystalline fraction at a given temperature corresponds to the maximum value of PJ (ω dr + ωH) in the relevant section of the dynamical phase diagram, as exemplified in (c). The error bars in (d) arise from the standard errors of Lorentzian fits to the blue-detuned side peaks. The parameters for the bilayer system are the same as in Fig. 3. The resonance frequencies are shifted at finite temperature.
the system, see Ref. [45] for example, while the dissipative coupling to the environment limits heating.
Next, we demonstrate the robustness of the Higgs time crystal against thermal fluctuations modelled as Langevin noise in the dynamics of the fields. These fluctuations are a natural test for the rigidity of the subharmonic response against temporal perturbations [13]. When considering thermal fluctuations, we include the in-plane dynamics of the fields in a full 3D simulation. The complete parameter set is summarised in Table I, implying the Higgs frequency ω H /2π ≈ 6.3 THz and the two longitudinal Josephson plasma frequencies ω J1 /2π ≈ 2.0 THz and ω J2 /2π ≈ 14.3 THz at T = 0. For simplicity, we keep the chemical potential fixed in the following finite-temperature calculations, µ(T ) ≡ µ. We choose the parameters within the CuO 2 planes to yield a critical temperature of T c ∼ 30 K. We find that a discretisation of 48 × 48 × 4 lattice sites with periodic boundaries is sufficient to obtain fully converged results with respect to the system size. Note that both the Higgs and Josephson plasma frequencies are renormalised at finite temperature [43].
Examples of the power spectra of the condensate amplitude and the supercurrent density at T = 3 K are shown in Figs. 4(a) and 4(b), respectively. When the sum resonance is driven, the condensate amplitude exhibits strong subharmonic modulation as evidenced by a sharp peak in the amplitude spectrum in Fig. 4(a). Moreover, we observe in Fig. 4(a) how the modulation of the condensate amplitude is suppressed as the driving frequency is tuned away from the resonance frequency. As shown in Fig. 4(b), we identify an experimentally relevant signature of the superconducting time crystalline phase, which is the appearance of two side peaks at ω dr ± ω H in the power spectrum of the supercurrent density. The side peaks vanish as the driving frequency is tuned away from the resonance frequency. Coherent dynamics of supercurrents can be experimentally probed using secondharmonic measurements [46,47].
To quantify the time crystalline fraction, we use the height of the blue-detuned side peak in the power spectrum of the supercurrent density, P J (ω dr + ω H ). Figure 4(d) displays the temperature dependence of the optimal crystalline fraction for a bilayer cuprate, normalised to the optimal time crystalline fraction at T = 0. The optimal driving parameters at each temperature were inferred from coarse scans such as that in Fig. 4(c). As we expect for time crystals under increasingly strong perturbation, the crystalline fraction decreases with temperature. Nevertheless, the subharmonic response is robust against thermal noise for temperatures up to T = 6 K ∼ 0.2 T c .

V. PULSED EXCITATION OF THE HIGGS TIME CRYSTAL
While significant progress has been made in generating continuous-wave terahertz sources [48], typical experiments in optically driven superconductors utilise pulsed excitation, as in most pump-probe experiments. We now point out that the time crystalline phase can be detected when the system is driven with a short pulse, rather than the steady state discussed so far. We consider a pulsed driving scheme by introducing a Gaussian envelope of the periodic driving, that is, E dr (t) = E 0 cos(ω dr t) exp(−t 2 /2σ 2 ) with the pulse width σ. In Fig. 5, we present an example of the dynamical response of the bilayer system under pulsed excitation. The response shown in Fig. 5(b) is approximately the Fourier broadened form of Fig. 4(b). The similarity between the two results suggests that the defining features 1 of the Higgs time crystal of continuously driven superconductors are detectable for pulsed driving protocols with realistic pulse lengths. The response can be clearly distinguished from normal dynamical phases by probing the coherent dynamics of the supercurrent. Thus, the Higgs time crystal predicted here can be observed in current state-of-the-art experiments with optically driven high-T c superconductors.

VI. DISCUSSION
In conclusion, we have demonstrated the emergence of a time-crystalline phase in a high-T c superconductor, which is induced by optical driving of a sum resonance of the Higgs mode and a Josephson plasma mode. Using a newly developed lattice gauge simulator, we demonstrate this time crystal for mono-and bilayer cuprates, and show its robustness against thermal fluctuations, for up to ∼ 20% of the critical temperature. As an experimentally accessible signature we observe the emergence of two side peaks at ω dr ±ω H in the supercurrent spectra. This signature is also visible in pulsed operation, which mimics realistic experimental conditions.
The emergent time crystalline order that we propose to induce, constitutes a qualitative departure from previous light-induced states in solids, because it is a genuine nonequilibrium state with no equilibrium counterpart. The realisation of such a state expands the scope of the scientific effort to design many-body states by optical driving beyond the paradigm of renormalising equilibrium orders. While even this existing paradigm has been and continues to be thought-provoking and stimulating, the work presented here urges the design and exploration of light-induced non-equilibrium states beyond that framework, and thereby expands the scope of the effort to design quantum matter on demand.

ACKNOWLEDGMENTS
We thank Andrea Cavalleri, Junichi Okamoto, and Kazuma Nagao for fruitful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of SFB 925 and the Cluster of Excellence "Advanced Imaging of Matter" (EXC 2056), Project No. 390715994.

Appendix A: Noise correlations
The fluctuation-dissipation theorem requires for the noise term of the superconducting field, where V 0 = d 2 ab d c is the discretisation volume of a single superconducting site. The noise correlations for the electric field are Appendix B: Characterisation of the response We characterise the response of the system to the periodic driving by studying the dynamics of the sample averages of the condensate amplitude and the supercurrent along the c-axis. The supercurrent along a single junction in the c-direction is given by the Josephson relation J z l,m,n = 4et n d c Im ψ * l,m,n+1 ψ l,m,n e ia z l,m,n .
The sample average of the supercurrent density along the c-axis can be obtained from where J s,w (t) denotes the spatial average of Josephson currents along either strong or weak junctions. In the case of non-zero temperatures, we average the power spectra P |ψ| (ω) and P J (ω) over an ensemble of trajectories. We find that 100 trajectories are enough to obtain convergent results for sampling thermal fluctuations at non-zero temperatures.   Table I summarises the parameters of our numerical calculations for mono-and bilayer systems, respectively. In both cases, our parameter choice of µ and g implies an equilibrium condensate density n 0 = µ/g = 2 × 10 21 cm −3 at T = 0. The bilayer system has two longitudinal c-axis plasma modes. Their eigenfrequencies are as follows from a sine-Gordon analysis at T = 0 [32,33].
Here we introduced the bare plasma frequencies of the strong and weak junctions where d c = (d s + d w )/2. The capacitive coupling constants are given by Besides, there is a transverse c-axis plasma mode with the eigenfrequency We have α s = 0.5, α w = 1, ω J1 /2π ≈ 2.0 THz, ω J2 /2π ≈ 14.3 THz, and ω T /2π ≈ 11.8 THz for the parameters specified in Table I. The in-plane plasma frequency amounts to 154 THz.
To study the behavior near the sum resonance, we write with the detuning δ. Inserting this into the second order equations induces secular terms, which we demand to vanish: The conditions (18) and (19) imply solutions of the form Using this ansatz, we find If the real part of r is positive, the amplitudes C J and C H grow exponentially. Such a behaviour signals the excitation of the sum resonance. It requires a sufficient driving amplitude given by the condition Let us consider the case ω dr = ω J + ω H , i.e., δ = 0. In this case, the required driving amplitude to induce the sum resonance is for the parameters specified in the main text. Higher order terms play an important role in saturating the amplitude of oscillations [1], which can be understood from the perspective of non-linear oscillators having amplitude dependent eigenfrequencies.
In the case of driving close to the difference frequency, we find Here the real part of r is always negative. Hence, there is no difference resonance in the system.

II. RIGIDITY OF THE HIGGS TIME CRYSTAL
The following zero-temperature simulations refer to the bilayer cuprate superconductor specified in the main text. We take the driving as where E 0 is the strength of the external field effectively penetrating the sample. Additionally, the external drive is characterised by the frequency ω dr and the rise time τ .
To realise the sum resonance of the Higgs mode and the upper Josephson plasmon, we drive the electric field with E 0 = 0.2 MV cm −1 and ω dr /2π = 21 THz. The long-time persistence of the time-translation symmetry breaking is exemplified in Fig. 1, where the subharmonic oscillations in the condensate amplitude are found to survive even after 10 5 driving cycles.
As discussed in Ref. [2], a signature of a phase transition to a time crystalline order in classical systems is the hysteretic behaviour across a critical point. Here, we demonstrate an indicator of such hysteresis in the response of the condensate amplitude. This can be seen in Fig. 2 as we tune the driving amplitude across the time crystalnormal response transition from E 0 = 0.08 MV cm −1 to E 0 = 0.1 MV cm −1 and vice versa, while keeping the driving frequency fixed at ω dr /2π = 21 THz. In particular, there is a clear difference in the time that it takes the system to enter and leave the time crystalline phase.

III. OPTICAL CONDUCTIVITY OF THE HIGGS TIME CRYSTAL
The optical conductivity σ(ω) is a crucial quantity to characterise the electric transport properties of a superconductor in the linear response regime. It is a macroscopic observable that can be measured in pump-probe experiments [3,4]. Here, we investigate how the emergence of the Higgs time crystal alters the c-axis optical conductivity of a bilayer cuprate superconductor with parameters as specified in the main text. For this purpose, the system is driven into the time crystalline phase with E 0 = 0.2 MV cm −1 and ω dr /2π = 21 THz at T = 0. Then, we add a probing term to the external drive, where t pr = 10 ps. The probing amplitude E pr has to be one order of magnitude smaller than E 0 to enter the linear response regime. We evaluate σ(ω) = J tot (ω)/E(ω) from a Fourier analysis over 50 ps in the steady state. The average electric field along the c-axis is given by where E s,w (t) denotes the spatial average of electric fields along either strong or weak junctions. The total current is the sum of the average supercurrent and the average displacement current inside the sample, that is, where J(t) is the supercurrent given in the main text. As visible in Fig. 3(a), the real part of the optical conductivity acquires additional resonance peaks in the time crystalline phase, especially at ω L = ω dr − ω H and ω R = ω dr + ω H . These frequencies correspond to the side peaks previously observed in the supercurrent spectra. Remarkably, the current response is amplified at the left side peak while attenuated at the right side peak. The counterparts of the peaks in σ 1 are sharp edges in σ 2 as evidenced by Fig. 3(b). The depletion of the condensate tends to reduce the plasma frequencies in the time crystalline phase. This effect is most apparent for the transverse Josephson plasmon shifting from 11.8 THz to 11.4 THz. For the same reason, we find a smaller prefactor of the 1/ω divergence of σ 2 at low frequencies.

IV. THERMAL PHASE TRANSITION
Here, we elaborate on the thermal phase transition of the simulated bilayer cuprate superconductor, see main text for parameters. The thermal equilibrium at a given temperature is established as follows. We initialise the system in its ground state at T = 0 and let the dynamics evolve without external driving, influenced only by thermal fluctuations and dissipation. After 10 ps, the average condensate density n = 1 N r |ψ r | 2 and the phase coherence are converged, indicating thermal equilbrium. To characterise the phase transition, we introduce the order parameter O = 1 n 1 N/2 l,m,n∈odd ψ * l,m,n+1 ψ l,m,n e ia z l,m,n .
The order parameter measures the phase coherence of the condensate across different bilayers. For each trajectory, it is evaluated from the average of 200 measurements within a time interval of 2 ps. Finally, we take the ensemble average of 100 trajectories. As shown in Fig. 4(a), the temperature dependence of the order parameter is reminiscent of a second order phase transition. Due to the finite size of the simulated system, the order parameter converges to a plateau with non-zero value for high temperatures. Instead of a sharp discontinuity, one finds a distinct crossover at T c ∼ 30 K. We also note that the lower Josephson plasmon vanishes in this temperature regime, which agrees with experimental observations [5]. Figure 4(b) reveals that the condensate density does not drop below 0.4 n 0 through the phase transition. Strikingly, the condensate density decreases almost linearly with temperature below T c . By contrast, it undergoes a nearly linear increase above the transition temperature. We see in Figs. 4(c) and 4(d) that the phase transition is only weakly modified by increasing the system size. with the corrections The estimate in Eq. (43) is compared to the purely numerical results in Fig. 6. For both curves, we take the ensemble average of 100 trajectories. The discrepancy between the the numerical and semi-analytical values can be ascribed to the approximations made in equations (37) and (38). Moreover, we have ignored the c-axis dynamics and higher order terms as present in the Mexican hat potential, for example. Nonetheless, our estimate distils the effect of the fourth order coupling terms ∼ a 2 h 2 on the Higgs frequency. For fixed model parameters, the Higgs frequency is shifted to a higher value because of thermally activated phase fluctuations in the ab-plane. In realistic systems, however, the chemical potential µ(T ) decreases with temperature such that the Higgs frequency does not necessarily increase.