Space-time crystalline order of a high-critical-temperature superconductor with intrinsic Josephson junctions

We theoretically demonstrate that the high-critical-temperature (high-Tc) superconductor Bi2Sr2CaCu2O8+x (BSCCO) is a natural candidate for the recently envisioned classical space-time crystal. BSCCO intrinsically forms a stack of Josephson junctions. Under a periodic parametric modulation of the Josephson critical current density, the Josephson currents develop coupled space-time crystalline order, breaking the continuous translational symmetry in both space and time. The modulation frequency and amplitude span a (nonequilibrium) phase diagram for a so-defined spatiotemporal order parameter, which displays rigid pattern formation within a particular region of the phase diagram. Based on our calculations using representative material properties, we propose a laser-modulation experiment to realize the predicted space-time crystalline behavior. Our findings bring new insight into the nature of space-time crystals and, more generally, into nonequilibrium driven condensed matter systems.

I n recent years, the notion of a space−time crystal (STC) has attracted a great deal of attention [1][2][3][4][5][6][7][8][9][10][11][12][13] . While there has been considerable discussion about what is truly outstanding in such a system, at present, it is mostly agreed that the STC refers to a nonequilibrium phase of matter displaying long-range order in both space and time 11 . Specifically, the (nonlinear) many-body interaction makes the system exhibit long-lived oscillations at a period longer than the period of the driving source, and the oscillation patterns show rigidity against perturbation from the environment 7,10,12,14 . These extraordinary behaviors have been theoretically conceived and experimentally observed in atomicmolecular-optical systems 3,[15][16][17][18][19] . But few candidates from condensed matter systems have been explored 20,21 .
In this paper, we present theoretically that the high-criticaltemperature (high-T c ) cuprate superconductor Bi 2 Sr 2 CaCu 2 O 8+x (BSCCO) is a natural candidate of a classical discrete STC that was recently envisioned 22 . This material, as illustrated in Fig. 1, acts as a stack of intrinsic Josephson junctions (IJJs) along its crystallographic c axis with s = 1.5 nm layer period [23][24][25] . Each junction is formed by the insulating BiO and SrO planes sandwiched between the superconducting CuO 2 planes. In the crystallographic ab plane, the junctions can be macroscopically large (≳1 mm) and behave as a long (and wide) Josephson junction. Adjacent junctions are coupled via phase gradients of the superconducting wave function, produced by currents flowing along the CuO 2 planes. The number N of junctions in a stack can vary from a few to thousands.
Our calculations indicate that when the critical current density j c of the junctions is subject to a parametric modulation, which is made periodic in time and uniform in space, BSCCO spontaneously develops half-harmonic oscillations of the Josephson currents in time and broken (continuous) translational symmetry in space. The newly formed order goes beyond a direct product of separate spatial and temporal orders and embodies space−time coupled symmetry 11 . At a fixed temperature T and thermal noise, the modulation frequency and amplitude span a phase diagram. Within a specific region with clear boundaries in this phase diagram, a nonzero spatiotemporal order parameter emerges, indicating the necessary robustness in phase formation.
Although parametric driving and half-harmonic generation of nonlinear oscillators such as Josephson junctions or junction arrays have been investigated previously [26][27][28] , their conceptual connection to the space−time crystalline order was only recently proposed 22 . Our study verifies that half-harmonic generation is accompanied with the spatiotemporal order formation in any spatial dimensions of BSCCO. But most notably, the STC phase is stable only at nonzero spatial dimension; the phase rigidity appears to increase with the spatial dimensionality. Therefore, our study of the intrinsic Josephson junction stack in BSCCO suggests the possible realization of a classical discrete STC from a naturally existing condensed matter system.

Results
Josephson plasma and photon modes in BSCCO. The dynamics of a BSCCO stack follows the inductively coupled sine-Gordon equations [29][30][31][32] . These are nonlinear differential equations for the gauge-invariant Josephson phase differences, parametrically dependent on the Josephson critical current density j c (T), the inplane and out-of-plane resistivities ρ ab (T) and ρ c (T), and the Cooper pair density n s (T). All of these quantities are T-dependent. The BSCCO stack supports collective oscillations of the supercurrents-the Josephson-plasma oscillations-across the junctions 33,34 . These oscillations can couple with the electromagnetic waves in the stack in the frequency range from ≲0.1 THz to ≳1 THz 24,25,[35][36][37] . In practice, BSCCO stacks are routinely fabricated into rectangular bars with tens-of-micron length L along x, several-micron width W along y, and a height H = Ns associated with tens or hundreds of junctions along z (cf. Fig. 1). This stack, resembling a slab waveguide cavity, hosts a series of resonance modes at gigahertz-to-terahertz frequencies.
For not too large oscillating amplitudes, the electric field and tunneling current along z on the nth junction take the form 38,39 , where x ∈ [0, L], y ∈ [0, W], and n = 1, 2, …, N. The mode indices l and m take zero or any positive integers, and q takes a positive integer from 1 to N. The associated resonance frequency reads with f p / ffiffiffiffiffiffiffiffiffiffi j c ðTÞ p being the Josephson plasma frequency and c q ¼ c= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À 2 s cosðπq=ðN þ 1ÞÞ p being the effective speed of light, where c is the Swihart velocity and s % 0:5 is the interlayer coupling constant. C ≤ 1 is a correction factor which we introduce here to account for a frequency shift generated by parametric driving. The values of j c , f p and c depend strongly on the charge carrier concentration. For convenience, we use j c ≈ 250 A cm −2 , f p ≈ 47 GHz, and c % 3:3 10 5 m/s at 4 K, which are typical for slightly underdoped BSCCO 31,32 . These quantities drop to zero at the critical temperature T c ≈ 85 K.
If l = m = 0, E z and j z are uniform in the x and y directions, independent of the lateral size L and W, and irrespective of the values of n and q in the thickness direction. When L and W are Each discretized element is characterized by the Josephson critical current I c , c-axis resistance R c , c-axis capacitance C c , ab-plane resistance R ab , and ab-plane inductance L ab , which is associated with the in-plane supercurrents and dominantly kinetic in origin. Additional geometric inductances could be added in series to R ab but are neglected. The discreteness along z is naturally set by the layered structure, whereas the element sizes along x and y are in theory infinitesimal, dx, dy → 0. In the numerical simulation, dx and dy are finite but smaller than the Josephson length λ J . Here we chose dx = dy = 0.5 μm, which is roughly a factor of 2 below λ J . c A typical rectangular-shaped 3D junction stack of length L, width W, and height H.
both small (≲2 μm), the only modes in the accessible frequency range are the uniform modes. In this case, if N = 1, we have an effectively 0-dimensional (0D) single-layer short junction. On the other hand, when L is large (≳10 μm) but W is small, we can differentiate two cases: If N = 1 and q = 1, we have an effectively 1D single-layer long junction; and if N > 1 and q = 1, 2, …, N, we have an effectively 2D multi-layer long junction stack. In these situations, l can take nonzero integers but m remains zero. Lastly, when L and W are both large and N > 1, we return to a 3D multilayer long and wide junction stack. Below, we show the emergence of space−time crystalline order in BSCCO with increasing dimensionality from 0D to 3D, under a temporally periodic but spatially uniform parametric modulation to j c . While the case of a 0D single-layer short junction is merely a replica of the children's swing problem 40,41 , exhibiting half-harmonic generation in time, the involvement of spatial dimensions from 1D to 3D not only exhibits half-harmonic generation but also stabilizes the temporal pattern on a much shorter time scale, and moreover, results in spontaneous (continuous) translational symmetry breaking in space. The parametric driving bandwidth and the robustness to thermal noise enhance significantly with the increase of spatial dimension. It is worthwhile to point out that there are three essential ingredients leading to the nontrivial behavior of this system: the discrete electromagnetic modes as described by Eq. (1), the nonlinear Josephson coupling of the superconducting planes, and the parametric drive through a periodic modulation of the Josephson critical current density. Superconductivity plays a key role in all these ingredients. Referring to Fig. 1b, in the normal state when superconductivity disappears, the inductive elements L ab would be absent (L ab → ∞) and the Josephson critical currents would be zero. Then the equivalent circuit would turn into a highly damped resistor network with linear capacitive couplings between the metallic layers. The mode spectrum as in Eq. (2) would be highly damped and, formally, f res would go to zero. For purely linear oscillators, parametric excitation still works. However, the oscillation amplitude may grow unbounded. We have explicitly tested this for the parameters of our system in the 0D and 1D cases by linearizing the Josephson current−phase relation. We found that the Josephson plasma oscillations or the cavity resonances were either not excited at all or grew exponentially. The sinusoidal Josephson current−phase relation removes this divergence.
Effectively 0D single-layer short junction. Let us first look at an effectively 0D single-layer short junction, whose governing equation is identical to that of a nonlinear pendulum 42,43 . If the critical current density j c is modulated in time, then the system is analogous to a children's swing. Our calculation takes account of the thermal (Nyquist) noise, which is related to the temperature T by the fluctuation-dissipation theorem. In the presence of some initial fluctuations, the current density j z across the junction starts to oscillate, provided that the modulation frequency is about twice the intrinsic Josephson plasma frequency f p and the modulation amplitude is large enough 40 . No inhomogeneous external force is applied to the system. Figure 2 shows our calculation for a single-layer short junction of BSCCO with the junction area L × W = 4 μm 2 . The parametric modulation takes the form where j c (T) is the unmodulated critical current density, f mod is the modulation frequency and A mod is the modulation amplitude. Figure 2a, b shows the time traces of j z (t), normalized to j c at 4 K, after many steps of initial relaxation with the same modulation frequency f mod ¼ 85 GHz % 2f p but two different modulation amplitudes A mod ¼ 0:01 (small) and 0.1 (large). For small modulation, the current amplitude is irregular. But for large modulation, the system can reach a Floquet steady state 14,44 , where the oscillation amplitude in time turns into an ordered pattern. Figure 2c, d gives the Fourier transform of j z (t) in terms of the amplitude |j z (f)| for the two modulation cases. A sharp half-harmonic peak around f p (along with a few other peaks at higher harmonics) can be observed in Fig. 2d for large modulation.
To quantify the temporal order, we define a (dimensionless) current−current correlation function, where the time average is taken in the Floquet steady state. The overline denotes additional ensemble averaging. Figure 2e, f shows g(δt), which is irregular for small modulation amplitude and almost periodic (like a sinusoidal function, as shown in the inset of Fig. 2f) for large modulation amplitude. To see the dependence of temporal order on the modulation parameters, we define a temporal order parameter, with the integer ν denoting even (2ν) and odd (2ν + 1) number of modulation period T mod in the time difference δt. The averaging is taken over both t and ν. The so-defined Δ t emphasizes oscillations of j z occurring at twice the modulation time 2T mod ¼ 2=f mod , since the odd and even parts have the same moduli but different signs. By contrast, it suppresses oscillations that are either periodic in T mod or incommensurate. Figure 2g gives a "phase diagram" of Δ t versus f mod and A mod at the temperature T = 4 K. We see that, like for the children's swing 40 , only if the modulation frequency falls in a narrow band near 2f p and the amplitude is neither too small nor too large, the system can produce ordered half-harmonics. If the modulation amplitude is overly large, the system can become chaotic. Figure 2h is a line-scan of Fig. 2g at the modulation amplitude A mod ¼ 0:1, showing a parametric excitation band within 80-90 GHz, near and slightly below 2f p = 94 GHz 40 . An important observation unique to the 0D case is that the correlation function and order parameter appear to persistently decay even after an exceedingly longtime drive. This means that the Floquet steady state here is in fact a quasi-steady state. In contrast, as shown below, after spatial dimensions are involved, temporal patterns are more efficiently stabilized to a true steady state. We deem this as a key difference between the well-known 0D single parametric oscillator and the more nontrivial multidimensional STC. (See "Methods" for a quantitative comparison between the 0D and 1D cases.) Effectively 1D single-layer long junction. The investigation above merely recovers the Floquet dynamics of a nonlinear parametric oscillator. However, by extending the studies into 1D, we find that even if the modulation still follows Eq. (3), the system permits spontaneous translational symmetry breaking. The system tends to pick favorable modes within a broadened modulation band that can be higher than 2f p . Figure 3 shows our calculations for an effectively 1D single-layer long Josephson junction in BSCCO with L = 25 μm and W = 2 μm. The spontaneously developed space−time crystalline order, especially, with rhombohedral unit cells, is remarkable. This symmetry is classified as C 2 m x m t in ref. 11 . To our knowledge, BSCCO provides the first material realization of this new symmetry class. Figure 3a gives the space−time trace j z (x, t) at f mod ¼ 140 GHz and A mod ¼ 0:1. Figure 3c gives its spatiotemporal Fourier transform in terms of the amplitude |j z (β, f)|, where β is the wavenumber (reciprocal of wavelength) along x. It shows a strong peak at the half-harmonic frequency f mod =2 and a wavenumber corresponding to mode index l = 9. This is manifested by the line projection of the Fourier amplitude across the half-harmonic peak versus frequency in Fig. 3b and versus wavenumber in Fig. 3d, respectively. We then define a spatiotemporal correlation function, where the average is taken for both space and time. Figure 3e gives a color plot for g(δx, δt), displaying nearly perfect sinusoidal oscillations along both δx and δt. To see the dependence of space −time crystalline order with the modulation parameters, we define a coupled spatiotemporal order parameter, as a natural extension of purely the temporal order parameter Eq. (5). This order parameter cannot be separated into a direct product of a spatial order parameter and temporal order parameter. It correlates the out-of-plane supercurrents at the most far distant points in space-the boundaries at x = 0 and x = L in 1D (and diagonally opposite corners in higher dimensions), consistent with our rectangular spatial system geometry. We find, because of the particular complexity of our system, this definition of the order parameter (or quasi order parameter) is practically the most rational and universal choice, owing to the fact that whenever a spatiotemporal order emerges and collective modes have formed, the long-range correlation function across the system must be large. In Fig. 3f, we plot the phase diagram of Δ xt versus f mod and A mod . Figure 3g gives a linescan of Fig. 3f when A mod ¼ 0:1. We can see that the half-harmonic response now extends into a broad frequency band from 120 to 160 GHz. Comparing with the 0D case, exhibiting an only 10 GHz wide band near and slightly below 2f p , the 1D system is more robust to the choice of the driving source. The peaks with Δ xt > 0.5 seen in Fig. 3g correspond to mode indices l between 3 and 10. Note that the fundamental mode with l = 0 when f mod % 2f p , which one may anticipate to show up, is suppressed.
Effectively 2D multi-layer long junction stack. Next, we include the spatial dimension along z by increasing the junction number N from 1 to 20. This dimension is distinct from x and y, and is unique for the naturally grown BSCCO crystal compared with the traditionally fabricated planar Nb or Al junction arrays 24,36 . Figure 4 gives our calculated results in this scenario. Figure 4a, b shows the space−time traces in the xt and zt plane at the fixed points z = 0 and x = 0, respectively. The modulation frequency is f mod ¼ 150 GHz and the amplitude is A mod ¼ 0:1. 2D crystalline order emerges in both directions. In Fig. 4c, d, we show the space −time Fourier transform for 4a, b respectively and plot the Fourier amplitudes in the wavenumber−frequency planes. In addition to the half-harmonic generation in t, the system spontaneously chooses the mode l = 6 in x and q = 6 in z. We can define the further generalized spatiotemporal correlation function, gðδx; δz; δtÞ ¼ hj z ðx; z; tÞj z ðx þ δx; z þ δz; t þ δtÞi x;z;t ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hj 2 z ðx; z; tÞi x;z;t hj 2 z ðx þ δx; z þ δz; t þ δtÞi x;z;t by taking averages in all dimensions and noting z = sn and δz = sδn. In Fig. 4e, f, we plot g(δx, 0, δt) and g(0, δz, δt), which again verifies the nearly perfect crystalline order between each pair of space−time dimensions. We can also define a further generalized coupled spatiotemporal order parameter by hj z ð0; 0; tÞj z ðL; H; t þ 2νT mod Þi t;ν ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hj 2 z ð0; 0; tÞi t;ν hj 2 z ðL; H; t þ 2νT mod Þi t;ν q À hj z ð0; 0; tÞj z ðL; H; t þ ð2ν þ 1ÞT mod Þi t;ν ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi hj 2 z ð0; 0; tÞi t;ν hj 2 z ðL; H; t þ ð2ν þ 1ÞT mod Þi t;ν q ; ð9Þ and plot its dependence on f mod at the fixed A mod ¼ 0:1, as shown in Fig. 4g. With the addition of the z-dimension, the system shows an almost continuous modulation bandwidth between 100 and 160 GHz.
3D multi-layer long and wide junction stack. Finally, we study the full 3D situation by assigning the y dimension a finite width W = 5 μm, and allowing for spatial variations of the Josephson phase differences along y. Such a width can be easily obtained in experiments 24 . Since in BSCCO, x and y are equivalent directions, we do not expect drastic changes of physics from the 2D case except for further enhanced robustness of space−time crystalline order due to the increased dimensionality. Our calculation does verify this expectation. Instead of repeating same sets of plots as above, we consider a possible experimental setup (see "Methods") and perform the 3D calculation with experimental parameters. We assume a laser beam, modulated in intensity at a frequency f mod , is absorbed by the BSCCO stack where it deposits an ac power P ac ¼ P ac;0 sin 2 ð2πf mod t=2Þ. If the BSCCO crystal is considered at thermal equilibrium, the result of the laser irradiation would just be a rise in temperature, without appreciable ac oscillations in the stack temperature, or phononic temperature. However, if the crystal is cooled sufficiently, one may see oscillations in the electronic temperature 45 , which is what we need.
In our calculation, we assume that the stack is glued to a 30 μm thick substrate, whose lower surface is kept at the bath temperature T b , which in practice is often lower than the temperature T of the junctions. The thickness of the glue layer is 30 μm. We assume a local equilibrium between the electrons and phonons in the crystal. To mimick oscillations in the electronic temperature we use extremely low values of the heat capacities of 5 J/m 3 K. We find, when starting at a bath temperature of 4 K and using P ac,0 = 0.5 mW, a time-averaged stack temperature of 29 K and an oscillation amplitude of 5.3 K, leading to a 10% modulation of j c , a 2% modulation of n s , a 20% modulation of ρ c and a 10% modulation of ρ ab . The time-averaged stack temperature is actually not constant across the stack but decays from 29.2 K in the center of the stack to 26.8 K at its edges. Figure 5a displays projections of the time-averaged intensity gð0; 0; 0; 0Þ ¼ hj 2 z ðx; y; z; tÞi=j c 2 of the normalized out-of-plane current densities j z flowing in the stack. The projections are to the (x, y), (x, z) and (y, z) planes, respectively. The modulation frequency is f mod = 154 GHz. The mode indices of the standing wave pattern are (l, m, q) = (8, 2, 13) and the oscillation frequency of this mode is 77 GHz.
Experimentally, to visualize the spatial order in the xy plane, one can use either near-field measurements of the ac electric fields near the stack surface 46,47 or low-temperature scanning laser microscopy 48 . The latter technique has already been used to visualize standing electromagnetic waves in BSCCO stacks used for the generation of THz radiation. In these experiments, the sample is usually current biased and a laser beam focused on the sample induces local heating which in turn changes the local electric properties. As a result, the global voltage across the sample changes. These global voltage changes, which are related to the local time average hj 2 z ðx; y; tÞi in our case, can be mapped by scanning. The typical resolution is around 1 μm, allowing to resolve many of the standing wave patterns along x and y. Fourier-transform infrared spectroscopy (FTIR) can be used to verify the predicted spatiotemporal ordered states and resolve the half-harmonic generation of the modulation frequency f mod in time. Figure 5b shows the expected spectrum. Figure 5c shows the calculated parametric modulation band as an indicator of the driving robustness.

Discussion
One may quest whether the spatiotemporal order in BSCCO can be extended to spatially arbitrarily large stacks. With increasing stack size, particularly along the z-direction, the number of resonant modes drastically increases, with a large number of modes that are nearly degenerate. A strong mode competition can occur in this case, eventually leading to a decreased mode stability and/or to chaotic behavior. Our preliminary simulations for larger junction numbers indeed indicate this trend. Therefore, a junction number of order 20 may already be close to the maximum in order not to destroy the space−time order. More indepth studies are needed along this direction.
Besides, we have not specifically mentioned the effect of bath temperature, which is proportional to the thermal fluctuation that may cause mode unlocking. At low temperatures around 4 K, within the timescales of our simulations, we have not observed any mode-unlocking instability. However, at high temperatures, c Fourier transform to j z (x, 0, t) plotted in the magnitude |j z (β, 0, f)|, where β is the wavenumber along x. d Fourier transform to j z (0, z, t) plotted in the magnitude |j z (0, γ, f)|, where γ is the wavenumber along z. e Space−time-averaged correlation function g(δx, δz, δt) in the δz = 0 plane. f Space−time-averaged correlation function in the δx = 0 plane. The corresponding mode indices are l = 6 and q = 6. g Spatiotemporal order parameter Δ xzt versus the modulation frequency f mod when A mod ¼ 0:1 (with fixed f p , temperature T = 4K, and thermal noise). It shows a more continuously filled parametric excitation band (compared with the 1D case) spanning 100−160 GHz for A mod ¼ 0:1. For A mod ¼ 0:1 the correction factor C = 0.91. e.g., above 30 K, we have indeed found that increased thermal fluctuations cause mode unlocking.
A fundamentally intriguing question is how the system will behave in the true thermodynamic limit, i.e., when the degrees of freedom, volume, and time all go to infinity. According to our simulations for a finite system, when the system degrees of freedom increase, the long-wavelength soft modes tend to be suppressed, likely due to the random destruction interference under strong drive. We thus speculate that in the true thermodynamic limit, only the modes with a spatial period on the order of the Josephson length λ J may survive.
In summary, we have shown theoretically that the high T c superconductor BSCCO is a natural candidate for a classical STC, owing to its property to intrinsically form a stack of long (and wide) Josephson junctions. Under a temporally periodic and spatially uniform modulation of the Josephson critical current density, BSCCO can display space−time crystalline order of the supercurrents that breaks the continuous translational symmetry in both space and time. If the size of BSCCO is comparable or after being averaged over (z, t), (x, t), and (y, t) and plotted in the xy, yz, and xz planes, respectively. The modulation frequency is f mod ¼ 154 GHz, the laser power is P = 0.5 mW, and the bath temperature is Expected Fourier spectrum of j z (x, y, z, t) on the edge of the sample (0, 0, 0) where it can convert to electromagnetic radiation and be detected. The spectrum shows the highest spark exactly at the half-harmonic frequency 77 GHz. c Expected parametric excitation band for the spatiotemporal order parameter for various modulation frequencies f mod at the fixed laser power 0.5 mW. This spectrum is further broadened (compared with the 2D case). smaller than the Josephson length in all directions, the system has effectively zero dimension and merely exhibits the well-known parametric oscillations in time and demands an infinite time to evolve towards a (quasi-)Floquet steady state. In contrast, finite (one to three) spatial dimensions ensure a rapid convergence of the system towards the real STC state. The modulation frequency and amplitude span a nonequilibrium phase diagram for a sodefined spatiotemporal order parameter. The phase diagram shows clear boundaries, outside of which the system is either disordered or chaotic. With increasing spatial dimensions, the system shows increasing stability manifested by a broadened modulation bandwidth.
Our calculations are all based on actual properties of BSCCO and so can be used to predict experiments. We have envisioned an experimental scheme to realize the BSCCO STC by using a near-infrared laser with a repetition rate of 100-200 GHz to modulate the Josephson current density. The signature of STC can be indirectly confirmed by measuring the half-harmonic emission spectrum via FTIR or directly visualized by the welldeveloped low-temperature scanning laser microscopy. If experimentally evidenced, this new kind of condensed-matter STCs not only enriches the classical and quantum many-body physics in nonequilibrium, but also holds promise for unprecedented applications in, e.g., tunable emitters and parametric amplifiers at terahertz frequencies.

Methods
Calculation scheme. The basic model describing the electromagnetic and thermal properties of IJJ stacks is given in ref. 31,32 . Here we give a short summary, with an emphasis on the geometry used for the present paper.
We consider a rectangular stack of N IJJs, each having a thickness (layer period) s = 1.5 nm. The thickness of the superconducting layers (CuO 2 double layers) is d s = 0.3 nm and the thickness of the barrier layers is d i = 1.2 nm. Like in the sketch of Fig. 1 the stack of length L and width W shall be a stand-alone stack, i.e., without a BSCCO base crystal underneath. The critical temperature of the stack is T c = 85 K. The electromagnetic part of the circuit is formulated in terms of in-plane and outof-plane current densities flowing, respectively, along the nth CuO 2 layer or across the nth junction in the stack.
The out-of-plane current densities j z,n across the nth IJJ between the CuO 2 planes n and n − 1 consist of Josephson currents with critical current density j c , (ohmic) quasiparticle currents with resistivity ρ c and displacement currents with dielectric constant ε. We also add Nyquist noise created by the quasiparticle currents. The in-plane currents along the nth CuO 2 plane consist of a superconducting part, characterized by a Cooper pair density n s and a quasiparticle component with resistivity ρ ab . The model parameters depend on temperature, as b Proposed setup to achieve 154 GHz laser modulation from 30.8 GHz using the rational harmonic generation technique and to detect the half-harmonic generation either with a Fourier-transform infrared spectroscopy or autocorrelation setup. Inset shows the loop gain profile with a sinusoidal line shape originating from the EOM modulation. The amplitude condition for lasing is satisfied when the gain profile equals unity (G L = 1).
Equations (10) and (11) are discretized in space, with X = 50 points in xdirection and Y points in y-direction. For 3D simulations Y = 10 and for 2D simulations Y = 1. The equations are propagated in time using a 5th order Runge-Kutta method. For each pixel, the normalized spectral density of η z,n is 4ΓXYρ c0 /ρ c , where Γ = 2πk B T/I c0 Φ 0 is the noise parameter.
In the 0D−2D simulations discussed in the main paper, we consider a minimal model where we assume that the stack is at a fixed temperature T, with a homogeneous temperature distribution inside the stack. For these simulations, we assume that the Josephson critical current density oscillates in time with an amplitude A mod , j c ðtÞ ¼ j c ½1 þ A mod cosð2πf mod tÞ=ð1 þ A mod Þ and j c (t) is homogeneous in space.
In the 3D simulations discussed in the main paper, we assume that some power P ac ¼ P ac0 cos 2 ð2πf mod t=2Þ is deposited in the stack and solve Eqs. (10) and (11) in combination with the heat diffusion equation cdT/dt = ∇ (κ ∇ T) + q s , with the specific heat capacity c and the (anisotropic and layer dependent) thermal conductivity κ. The power density q s for heat generation in the stack results from Joule heating due to the electrical part of the circuit and from the ac power deposited by the laser. For the simulations, we assume that the stack is glued to a 30 μm thick substrate with lateral dimensions 2L × 2W, whose lower surface is kept at the bath temperature T b . The thickness of the glue layer is 30 μm. To mimick oscillations in the electronic temperature we use extremely low values for the heat capacities of 5 J/m 3 K. The thermal equations are discretized in space, using 2X points along x and 2Y points along y and are propagated in time using a 5th order Runge-Kutta method.
It is important to note the difference between 0D and finite spatial dimensions. For 0D, within our calculated time scale, the system has never reached a true steady-state, i.e., a desired time-crystal phase. As shown in Fig. 6a, after the system passes the transient state and enters the Floquet (quasi-)steady state, the temporal correlation function persistently decays even by a time difference of 32,768 cycles of modulation period T mod ¼ 1=f mod , with f mod ¼ 85 GHz. In contrast, for 1D, as shown in Fig. 6b, the temporal correlation (evaluated at a specific point with zero space difference) establishes quickly and remains the same up to 32768 cycles of T mod ¼ 1=f mod , with f mod ¼ 140 GHz.
Experimental proposal. A schematic experimental design to generate periodic modulations of the junction parameters is given in Fig. 7a. The major challenge is to achieve a very high repetition rate 100−200 GHz and a strong power so as to induce about 10% change of j c , i.e., A mod ¼ 0:1. One may consider raising the bath temperature T b close to T c to reduce f p and thus alleviate the frequency requirement. However, the simultaneously increased thermal noise significantly suppresses the formation of space−time crystalline order. In fact, in our simulations, we could not find a substantial spatiotemporal order above~40 K even in the 3D case.
Laser pulses with a high repetition rate can be achieved by adopting the rational harmonic mode-locking (RHML) laser technique 49,50 . As shown in Fig. 7b, the construction of a RHML laser consists of a closed-loop fiber ring as the optical cavity. An erbium doped-fiber amplifier (EDFA) is plugged inside the fiber ring to provide gain for lasing in the telecommunication C-band. In general, the lasing frequency (wavelength) of a ring laser is determined by two conditions: (1) amplitude condition: loop gain G L = 1; (2) phase condition: loop phase φ L = 2πM with M being an arbitrary integer. The amplitude condition is determined by a lithium niobate electro-optical modulator (EOM) inserted into the ring, which modulates the loop gain profile and consequently determines the lasing frequency, while the phase condition is controlled by carefully selecting the modulation frequency. When the EOM operates at a frequency f m ¼ s þ 1=p À Á f c , where s and p are integers, and f c is the fundamental frequency of the fiber ring cavity which is determined by the loop length of the fiber ring, the repetition rate of the laser pulses can be obtained as f r = (sp + 1)f c . As an estimation, a 5 m long single mode (SM) fiber loop gives f c = 40 MHz. Choosing s = 770 and p = 5, a high repetition rate of f r = 154 GHz can be obtained, while the required modulation frequency f m = 30.8 GHz is readily accessible using standard EOM and microwave sources. In practice, the polarization of the laser is determined by the polarization-maintaining elements (such as the EOM) inside the ring. The polarization controller (PC) is used to optimize the polarization condition for these elements, and an optical isolator is used to reinforce unidirectional propagation of the laser light inside the ring. Most importantly, such a configuration ensures that the modulation signal, laser repetition rate, and half-harmonic generation are all at different frequencies, and therefore potential crosstalk in the characterization of the generated halfharmonic signal can be drastically reduced.