## Abstract

We theoretically demonstrate that the high-critical-temperature (high-*T*_{c}) superconductor Bi_{2}Sr_{2}CaCu_{2}O_{8+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.

## Introduction

In 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 atomic-molecular-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-critical-temperature (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 *a**b* 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 in-plane 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* = *N**s* 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 *n*th 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}_{{{{{{{{\rm{p}}}}}}}}}\propto \sqrt{{j}_{{{\mbox{c}}}}(T)}\) being the Josephson plasma frequency and \({c}_{q}=\bar{c}/\sqrt{1-2\bar{s}\cos (\pi q/(N+1))}\) being the effective speed of light, where \(\bar{c}\) is the Swihart velocity and \(\bar{s}\approx 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 \(\bar{c}\) depend strongly on the charge carrier concentration. For convenience, we use *j*_{c} ≈ 250 A cm^{−2}, *f*_{p} ≈ 47 GHz, and \(\bar{c}\approx 3.3\times 1{0}^{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 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 multi-layer 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}_{{{{{{{{\rm{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}_{{{{{{{{\rm{mod}}}}}}}}}\) is the modulation frequency and \({A}_{{{{{{{{\rm{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}_{{{{{{{{\rm{mod}}}}}}}}}=85\ \,{{\mbox{GHz}}}\,\approx 2{f}_{{{{{{{{\rm{p}}}}}}}}}\) but two different modulation amplitudes \({A}_{{{{{{{{\rm{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}_{{{{{{{{\rm{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 \(2{T}_{{{{{{{{\rm{mod}}}}}}}}}=2/{f}_{{{{{{{{\rm{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}_{{{{{{{{\rm{mod}}}}}}}}}\) or incommensurate. Figure 2g gives a “phase diagram” of Δ_{t} versus \({f}_{{{{{{{{\rm{mod}}}}}}}}}\) and \({A}_{{{{{{{{\rm{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 2*f*_{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}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\), showing a parametric excitation band within 80–90 GHz, near and slightly below 2*f*_{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 2*f*_{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}_{{{{{{{{\rm{mod}}}}}}}}}=140\) GHz and \({A}_{{{{{{{{\rm{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}_{{{{{{{{\rm{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}_{{{{{{{{\rm{mod}}}}}}}}}\) and \({A}_{{{{{{{{\rm{mod}}}}}}}}}\). Figure 3g gives a linescan of Fig. 3f when \({A}_{{{{{{{{\rm{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 2*f*_{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}_{{{{{{{{\rm{mod}}}}}}}}}\approx 2{f}_{{{{{{{{\rm{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 *x**t* and *z**t* plane at the fixed points *z* = 0 and *x* = 0, respectively. The modulation frequency is \({f}_{{{{{{{{\rm{mod}}}}}}}}}=150\) GHz and the amplitude is \({A}_{{{{{{{{\rm{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,

by taking averages in all dimensions and noting *z* = *s**n* 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

and plot its dependence on \({f}_{{{{{{{{\rm{mod}}}}}}}}}\) at the fixed \({A}_{{{{{{{{\rm{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}_{{{{{{{{\rm{mod}}}}}}}}}\), is absorbed by the BSCCO stack where it deposits an ac power \({P}_{{{{{{{{\rm{ac}}}}}}}}}={P}_{{{{{{{{\rm{ac,0}}}}}}}}}{\sin }^{2}(2\pi {f}_{{{{{{{{\rm{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)=\langle {j}_{z}^{2}(x,y,z,t)\rangle /{j_{{{\rm{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}_{{{{{{{{\rm{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 *x**y* 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 \(\langle {j}_{z}^{2}(x,y,t)\rangle\) 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}_{{{{{{{{\rm{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 in-depth 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, 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 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 so-defined 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 well-developed 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 out-of-plane current densities flowing, respectively, along the *n*th CuO_{2} layer or across the *n*th junction in the stack.

The out-of-plane current densities *j*_{z,n} across the *n*th 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 *n*th 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 plotted in detail in ref. ^{31}. Then, the electric part of the circuit is described by sine-Gordon-like equations for the Josephson phase differences *γ*_{n}(*x*, *y*) in the *n*th junction of the IJJ stack:

Here, *n* = 1…*N*, ∇ = (∂/∂*x*, ∂/∂*y*) and \({\lambda }_{k}={[{{{\Phi }}}_{0}{d}_{{{{{{{{\rm{s}}}}}}}}}/(2\pi {\mu }_{0}{j}_{{{{{{{{\rm{c0}}}}}}}}}{\lambda }_{ab0}^{2})]}^{1/2}\), with the in-plane London penetration depth *λ*_{ab0} at 4 K and the magnetic permeability *μ*_{0}. \({\lambda }_{{{{{{{{\rm{c}}}}}}}}}={[{{{\Phi }}}_{0}/(2\pi {\mu }_{0}{j}_{{{{{{{{\rm{c0}}}}}}}}}s)]}^{1/2}\) is the out-of-plane penetration depth. *j*_{c0} is the 4 K Josephson critical current density and *ρ*_{c0} is the BSCCO *c* axis (subgap-)resistivity at 4 K. Φ_{0} is the flux quantum. For the out-of-plane current densities,

with \({\beta }_{{{{{{{{\rm{c0}}}}}}}}}=2\pi {j}_{{{{{{{{\rm{c0}}}}}}}}}{\rho }_{c0}^{2}\varepsilon {\varepsilon }_{0}s/{{{\Phi }}}_{0}\); *ε*_{0} is the vacuum permittivity and *η*_{z,n} are the out-of-plane noise current densities. Time is normalized to Φ_{0}/2*π**j*_{c0}*ρ*_{c0}*s*, resistivities to *ρ*_{c0} and current densities to *j*_{c0}.

For the electrical parameters we use the 4 K values *ρ*_{c0} = 1100 Ω cm, *j*_{c0} = 250 A/cm^{2}, *λ*_{ab0} = 260 nm, and further *ρ*_{ab}(*T*_{c}) = 100 *μ*Ωcm and *ε* = 13, yielding *λ*_{c} = 264 *μ*m, *λ*_{k} = 1.52 μm and *β*_{c0} = 1.5 × 10^{5}. We will further make use of the Josephson length *λ*_{J}, obtained via \({{\lambda }_{{{{{{{{\rm{J}}}}}}}}}}^{-1}=2{\lambda }_{k}^{-1}+{\lambda }_{c}^{-1}\). For our parameters *λ*_{J} ~ 1.1 μm. For the temperature dependence of *j*_{c} we use a parabolic profile, \({j}_{{{{{{{{\rm{c}}}}}}}}}\propto 1-{(T/{T}_{{{{{{{{\rm{c}}}}}}}}})}^{2}\); for the temperature dependence of *ρ*_{c} see ref. ^{31}.

Equations (10) and (11) are discretized in space, with *X* = 50 points in *x-*direction 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Γ*X**Y**ρ*_{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}_{{{{{{{{\rm{mod}}}}}}}}}\), \({j}_{{{{{{{{\rm{c}}}}}}}}}(t)={j}_{{{{{{{{\rm{c}}}}}}}}}[1+{A}_{{{{{{{{\rm{mod}}}}}}}}}\cos (2\pi {f}_{{{{{{{{\rm{mod}}}}}}}}}t)]/(1+{A}_{{{{{{{{\rm{mod}}}}}}}}})\) and *j*_{c}(*t*) is homogeneous in space.

In the 3D simulations discussed in the main paper, we assume that some power \({P}_{{{{{{{{\rm{ac}}}}}}}}}={P}_{{{{{{{{\rm{ac0}}}}}}}}}{\cos }^{2}(2\pi {f}_{{{{{{{{\rm{mod}}}}}}}}}t/2)\) is deposited in the stack and solve Eqs. (10) and (11) in combination with the heat diffusion equation *c*d*T*/d*t* = ∇ (*κ* ∇ *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 2*L* × 2*W*, 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 2*X* points along *x* and 2*Y* 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}_{{{{{{{{\rm{mod}}}}}}}}}=1/{f}_{{{{{{{{\rm{mod}}}}}}}}}\), with \({f}_{{{{{{{{\rm{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}_{{{{{{{{\rm{mod}}}}}}}}}=1/{f}_{{{{{{{{\rm{mod}}}}}}}}}\), with \({f}_{{{{{{{{\rm{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}_{{{{{{{{\rm{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}_{{{\rm m}}}=\left(s+1/p\right){f}_{{{\rm 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} = (*s**p* + 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 half-harmonic signal can be drastically reduced.

## Data availability

The authors declare that all data supporting the findings of this study are available within the paper.

## Code availability

Computer codes are available from the corresponding authors upon reasonable request.

## References

- 1.
Wilczek, F. Quantum time crystals.

*Phys. Rev. Lett.***109**, 160401 (2012). - 2.
Shapere, A. & Wilczek, F. Classical time crystals.

*Phys. Rev. Lett.***109**, 160402 (2012). - 3.
Li, T. C. et al. Space−time crystals of trapped ions.

*Phys. Rev. Lett.***109**, 163001 (2012). - 4.
Bruno, P. Impossibility of spontaneously rotating time crystals: a no-go theorem.

*Phys. Rev. Lett.***111**, 070402 (2013). - 5.
Watanabe, H. & Oshikawa, M. Absence of quantum time crystals.

*Phys. Rev. Lett.***114**, 251603 (2015). - 6.
Lazarides, A., Das, A. & Moessner, R. Periodic thermodynamics of isolated quantum systems.

*Phys. Rev. Lett.***112**, 150401 (2014). - 7.
Sacha, K. Modeling spontaneous breaking of time-translation symmetry.

*Phys. Rev. A***91**, 033617 (2015). - 8.
Lazarides, A., Das, A. & Moessner, R. Fate of many-body localization under periodic driving.

*Phys. Rev. Lett.***115**, 030402 (2015). - 9.
Khemani, V., Lazarides, A., Moessner, R. & Sondhi, S. L. Phase structure of driven quantum systems.

*Phys. Rev. Lett.***116**, 250401 (2016). - 10.
Yao, N. Y., Potter, A. C., Potirniche, I.-D. & Vishwanath, A. Discrete time crystals: rigidity, criticality, and realizations.

*Phys. Rev. Lett.***118**, 030401 (2017). - 11.
Xu, S. & Wu, C. Space−time crystal and space−time group.

*Phys. Rev. Lett.***120**, 096401 (2018). - 12.
Sacha, K. & Zakrzewski, J. Time crystals: a review.

*Rep. Prog. Phys.***81**, 016401 (2018). - 13.
Kozin, V. K. & Kyriienko, O. Quantum time crystals from Hamiltonians with long-range interactions.

*Phys. Rev. Lett.***123**, 210602 (2019). - 14.
Else, D. V., Bauer, B. & Nayak, C. Floquet time crystals.

*Phys. Rev. Lett.***117**, 090402 (2016). - 15.
Smits, J., Liao, L., Stoof, H. T. C. & van der Straten, P. Observation of a space−time crystal in a superfluid quantum gas.

*Phys. Rev. Lett.***121**, 185301 (2018). - 16.
Zhang, J. et al. Observation of a discrete time crystal.

*Nature***543**, 217 (2017). - 17.
Choi, S. W. et al. Observation of discrete time-crystalline order in a disordered dipolar many-body system.

*Nature***543**, 221 (2017). - 18.
Rovny, J., Blum, R. L. & Barrett, S. E. Observation of discrete-time-crystal signatures in an ordered dipolar many-body system.

*Phys. Rev. Lett.***120**, 180603 (2018). - 19.
Pal, S., Nishad, N., Mahesh, T. S. & Sreejith, G. J. Temporal order in periodically driven spins in star-shaped clusters.

*Phys. Rev. Lett.***120**, 180602 (2018). - 20.
Autti, S., Eltsov, V. B. & Volovik, G. E. Observation of a time quasicrystal and its transition to a superfluid time crystal.

*Phys. Rev. Lett.***120**, 215301 (2018). - 21.
Homann, G., Cosme, J. G. & Mathey, L. Higgs time crystal in a high-

*T*_{c}superconductor.*Phys. Rev. Res.***2**, 043214 (2020). - 22.
Yao, N. Y., Nayak, C., Balents, L. & Zaletel, M. P. Classical discrete time crystals.

*Nat. Phys.***16**, 438–447 (2020). - 23.
Kleiner, R., Steinmeyer, F., Kunkel, G. & Müller, P. Intrinsic Josephson effects in Bi

_{2}Sr_{2}CaCu_{2}O_{8}single crystals.*Phys. Rev. Lett.***68**, 2394 (1992). - 24.
Hu, X. & Lin, S. Z. Phase dynamics in a stack of inductively coupled intrinsic Josephson junctions and terahertz electromagnetic radiation.

*Supercond. Sci. Technol.***23**, 053001 (2010). - 25.
Savel’ev, S., Yampol’skii, V. A., Rakhmanov, A. L. & Nori, F. Terahertz Josephson plasma waves in layered superconductors: spectrum, generation, nonlinear and quantum phenomena.

*Rep. Prog. Phys.***73**, 026501 (2010). - 26.
Feldman, M. J., Parrish, P. T. & Chiao, R. Y. Parametric amplification by unbiased Josephson junctions.

*J. Appl. Phys.***46**, 4031 (1975). - 27.
Pedersen, N. F., Soerensen, O. H., Dueholm, B. & Mygind, J. Half-harmonic parametric oscillations in Josephson junctions.

*J. Low. Temp. Phys.***38**, 1 (1980). - 28.
Pedersen, N. F. & Sakai, S. Prediction of half harmonic generation in stacked Josephson junctions and Bi

_{2}Sr_{2}CaCu_{2}O_{x}single crystals.*Phys. Rev. B***61**, 11328 (2000). - 29.
Sakai, S., Bodin, P. & Pedersen, N. F. Fluxons in thin-film superconductor−insulator superlattices.

*J. Appl. Phys.***73**, 2411 (1993). - 30.
Bulaevskii, L. N., Zamora, M., Baeriswyl, D., Beck, H. & Clem, J. R. Time-dependent equations for phase differences and a collective mode in Josephson-coupled layered superconductors.

*Phys. Rev. B***50**, 12831 (1994). - 31.
Rudau, F. et al. Thermal and electromagnetic properties of Bi

_{2}Sr_{2}CaCu_{2}O_{8}intrinsic Josephson junction stacks studied via one-dimensional coupled sine-Gordon equations.*Phys. Rev. B***91**, 104513 (2015). - 32.
Rudau, F. et al. Three-dimensional simulations of the electrothermal and terahertz emission properties of Bi

_{2}Sr_{2}CaCu_{2}O_{8}intrinsic Josephson junction stacks.*Phys. Rev. Appl.***5**, 044017 (2016). - 33.
Tsui, O. K., Ong, N. P., Matsuda, Y., Yan, Y. F. & Peterson, J. B. Sharp magnetoabsorption resonances in the vortex state of Bi

_{2}Sr_{2}CaCu_{2}O_{8+δ}.*Phys. Rev. Lett.***73**, 724 (1994). - 34.
Matsuda, Y., Gaifullin, M. B., Kumagai, K., Kadowaki, K. & Mochiku, T. Collective Josephson plasma resonance in the vortex state of Bi

_{2}Sr_{2}CaCu_{2}O_{8+δ}.*Phys. Rev. Lett.***75**, 4512 (1995). - 35.
Ozyuzer, L. et al. Emission of coherent THz radiation from superconductors.

*Science***318**, 1291 (2007). - 36.
Welp, U., Kadowaki, K. & Kleiner, R. Superconducting emitters of THz radiation.

*Nat. Photonics***7**, 702 (2013). - 37.
Kakeya, I. & Wang, H. B. Terahertz-wave emission from Bi2212 intrinsic Josephson junctions: a review on recent progress.

*Supercond. Sci. Technol.***29**, 073001 (2016). - 38.
Kleiner, R. Two-dimensional resonant modes in stacked Josephson junctions.

*Phys. Rev. B***50**, 6919 (1994). - 39.
Sakai, S., Ustinov, A. V., Kohlstedt, H., Petraglia, A. & Pedersen, N. F. Theory and experiment on electromagnetic-wave-propagation velocities in stacked superconducting tunnel structures.

*Phys. Rev. B***50**, 12905 (1994). - 40.
Belyakov, A. O., Seyranian, A. P. & Luongo, A. Dynamics of the pendulum with periodically varying length.

*Physica D***238**, 1589 (2009). - 41.
Glendinning, P. Adaptive resonance and pumping a swing.

*Eur. J. Phys.***41**, 025006 (2020). - 42.
Stewart, W. C. Current−voltage characteristics of Josephson junctions.

*Appl. Phys. Lett.***12**, 277 (1968). - 43.
McCumber, D. E. Effect of ac impedance on dc voltage−current characteristics of superconductor weak-link junctions.

*J. Appl. Phys.***39**, 3113 (1968). - 44.
Oka, T. & Kitamura, S. Floquet engineering of quantum materials.

*Annu. Rev. Condens. Matter Phys.***10**, 387–408 (2019). - 45.
Sobolewski, R. Ultrafast dynamics of nonequilibrium quasi-particles in hight-temperature superconductors.

*Proc. SPIE***3481**, 480 (1998). - 46.
Stewing, F., Brendel, C. & Schilling, M. Three dimensional near-field radiation imaging up to the THz-regime.

*Frequenz***62**, 149 (2008). - 47.
Adam, A. J. L. Review of near-field terahertz measurement methods and their applications—how to achieve sub-wavelength resolution at THz frequencies.

*J. Infrared Milli. Terahz Waves***32**, 976 (2011). - 48.
Wang, H. B. et al. Hot spots and waves in Bi

_{2}Sr_{2}CaCu_{2}O_{8}intrinsic Josephson junction stacks: a study by low temperature scanning laser microscopy.*Phys. Rev. Lett.***102**, 017006 (2009). - 49.
Yoshida, E. & Nakazawa, M. 80~200 GHz erbium doped fibre laser using a rational harmonic mode-locking technique.

*Electron. Lett.***32**, 1370–1372 (1996). - 50.
Das, P. K., Kaechele, W., Theimer, J. P., & Pirich, A. R. Rational harmonic mode-locking fiber laser,

*Photonic Processing Technology and Applications, Proc. SPIE***3075**, 21−32 (1997).

## Acknowledgements

This work was performed in part at the Center for Nanoscale Materials, a U.S. Department of Energy Office of Science User Facility, and supported by the U.S. Department of Energy, Office of Science, under Contract No. DE-AC02-06CH11357. E. D., D.K., and R.K. acknowledge support by the Deutsche Forschungsgemeinschaft via project KL930-13/2, and the EU-FP6-COST Action NANOCOHYBRI (CA16218).

## Author information

### Affiliations

### Contributions

R.K. and D.J. conceived the idea and led the project. R.K. performed the numerical calculation supported by E.D.X. Zhou analyzed the results and developed the experimental design. X.Z. and D.K. provided additional advices and joined discussions. All authors contributed to the paper writing. R.K. and X.Z. contributed equally to this work.

### Corresponding authors

## Ethics declarations

### Competing interests

Authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks Congjun Wu and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Kleiner, R., Zhou, X., Dorsch, E. *et al.* Space-time crystalline order of a high-critical-temperature superconductor with intrinsic Josephson junctions.
*Nat Commun* **12, **6038 (2021). https://doi.org/10.1038/s41467-021-26132-y

Received:

Accepted:

Published:

DOI: https://doi.org/10.1038/s41467-021-26132-y

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.