Skip to main content

Thank you for visiting You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

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.


In recent years, the notion of a space−time crystal (STC) has attracted a great deal of attention1,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 time11. 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 environment7,10,12,14. These extraordinary behaviors have been theoretically conceived and experimentally observed in atomic-molecular-optical systems3,15,16,17,18,19. But few candidates from condensed matter systems have been explored20,21.

In this paper, we present theoretically that the high-critical-temperature (high-Tc) cuprate superconductor Bi2Sr2CaCu2O8+x (BSCCO) is a natural candidate of a classical discrete STC that was recently envisioned22. 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 period23,24,25. Each junction is formed by the insulating BiO and SrO planes sandwiched between the superconducting CuO2 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 CuO2 planes. The number N of junctions in a stack can vary from a few to thousands.

Fig. 1: Schematic of the superconductor BSCCO with intrinsic Josephson junctions.
figure 1

a Atomic structure in which the CuO2 planes are superconducting. b Equivalent circuit. Each discretized element is characterized by the Josephson critical current Ic, c-axis resistance Rc, c-axis capacitance Cc, ab-plane resistance Rab, and ab-plane inductance Lab, which is associated with the in-plane supercurrents and dominantly kinetic in origin. Additional geometric inductances could be added in series to Rab 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.

Our calculations indicate that when the critical current density jc 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 symmetry11. 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 previously26,27,28, their conceptual connection to the space−time crystalline order was only recently proposed22. 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.


Josephson plasma and photon modes in BSCCO

The dynamics of a BSCCO stack follows the inductively coupled sine-Gordon equations29,30,31,32. These are nonlinear differential equations for the gauge-invariant Josephson phase differences, parametrically dependent on the Josephson critical current density jc(T), the in-plane and out-of-plane resistivities ρab(T) and ρc(T), and the Cooper pair density ns(T). All of these quantities are T-dependent. The BSCCO stack supports collective oscillations of the supercurrents—the Josephson-plasma oscillations—across the junctions33,34. These oscillations can couple with the electromagnetic waves in the stack in the frequency range from 0.1 THz to 1 THz24,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 form38,39,

$${E}_{z},{j}_{z}\propto \cos \left(\frac{\pi lx}{L}\right)\cos \left(\frac{\pi my}{W}\right)\sin \left(\frac{\pi qn}{N+1}\right),$$

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

$${f}_{{{{{{{{\rm{res}}}}}}}}}=\sqrt{{\left(C{f}_{{{{{{{{\rm{p}}}}}}}}}\right)}^{2}+{c}_{q}^{2}\left[{\left(\frac{\pi l}{L}\right)}^{2}+{\left(\frac{\pi m}{W}\right)}^{2}\right]},$$

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 jc, fp and \(\bar{c}\) depend strongly on the charge carrier concentration. For convenience, we use jc ≈ 250 A cm−2, fp ≈ 47 GHz, and \(\bar{c}\approx 3.3\times 1{0}^{5}\) m/s at 4 K, which are typical for slightly underdoped BSCCO31,32. These quantities drop to zero at the critical temperature Tc ≈ 85 K.

If l = m = 0, Ez and jz 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 jc. While the case of a 0D single-layer short junction is merely a replica of the children’s swing problem40,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 Lab would be absent (Lab → ) 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 pendulum42,43. If the critical current density jc 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 jz across the junction starts to oscillate, provided that the modulation frequency is about twice the intrinsic Josephson plasma frequency fp and the modulation amplitude is large enough40. 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 μm2. The parametric modulation takes the form

$$j_{{{{{{{\rm{c}}}}}}}}(T,t)=j_{{{{{{{\rm{c}}}}}}}}(T)\frac{1+{A}_{{{{{{{{\rm{mod}}}}}}}}}\cos (2\pi {f}_{{{{{{{{\rm{mod}}}}}}}}}t)}{1+{A}_{{{{{{{{\rm{mod}}}}}}}}}},$$

where jc(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 jz(t), normalized to jc 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 state14,44, where the oscillation amplitude in time turns into an ordered pattern. Figure 2c, d gives the Fourier transform of jz(t) in terms of the amplitude jz(f) for the two modulation cases. A sharp half-harmonic peak around fp (along with a few other peaks at higher harmonics) can be observed in Fig. 2d for large modulation.

Fig. 2: Calculation on an effectively 0D single-layer short Josephson junction in BSCCO at 4 K under parametric modulation.
figure 2

The junction area is L × W = 4 μm2 and the Josephson plasma frequency is fp = 47 GHz. a Time trace of the current density jz(t) under a modulation frequency \({f}_{{{{{{{{\rm{mod}}}}}}}}}=85\) GHz and small modulation amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.01\). The system does not converge to a Floquet steady state. b Time trace of jz(t) under a large modulation amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\). The system has reached the Floquet steady state. c Fourier transform to jz(t) in (a), plotted in the magnitude of jz(f). d Fourier transform to jz(t) in (b), showing a resolution-limited sharp half-harmonic generation around the Josephson plasma frequency. e Time-averaged current−current correlation function g(δt) for \({f}_{{{{{{{{\rm{mod}}}}}}}}}=85\) GHz and small modulation amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.01\). f Time-averaged correlation function g(δt) for a large modulation amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\), showing an ordered temporal pattern. g Phase diagram of temporal order parameter Δt for an ensemble of 50 junctions versus the modulation frequency \({f}_{{{{{{{{\rm{mod}}}}}}}}}\) and amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}\) (with fixed fp, temperature T = 4 K, and thermal noise). Clear phase boundaries can be identified. h Temporal order parameter versus \({f}_{{{{{{{{\rm{mod}}}}}}}}}\) when \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\), showing a parametric excitation band within 80−90 GHz.

To quantify the temporal order, we define a (dimensionless) current−current correlation function,

$$g(\delta t)=\frac{{\langle \overline{{j}_{z}(t){j}_{z}(t+\delta t)}\rangle }_{t}}{\sqrt{{\langle \overline{{j}_{z}^{2}(t)}\rangle }_{t}{\langle \overline{{j}_{z}^{2}(t+\delta t)}\rangle }_{t}}},$$

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,

$${{{\Delta }}}_{t}=\; \frac{1}{2}\left[\frac{{\langle {j}_{z}(t){j}_{z}(t+2\nu {T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}{\sqrt{{\langle {j}_{z}^{2}(t)\rangle }_{t,\nu }{\langle {j}_{z}^{2}(t+2\nu {T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}}\right.\\ \left.-\frac{{\langle {j}_{z}(t){j}_{z}(t+(2\nu +1){T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}{\sqrt{{\langle {j}_{z}^{2}(t)\rangle }_{t,\nu }{\langle {j}_{z}^{2}(t+(2\nu +1){T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}}\right],$$

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 jz 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 swing40, only if the modulation frequency falls in a narrow band near 2fp 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 2fp = 94 GHz40.

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 2fp. 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 C2mxmt in ref. 11. To our knowledge, BSCCO provides the first material realization of this new symmetry class. Figure 3a gives the space−time trace jz(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 jz(β, 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,

$$g(\delta x,\delta t)=\frac{{\langle {j}_{z}(x,t){j}_{z}(x+\delta x,t+\delta t)\rangle }_{x,t}}{\sqrt{{\langle {j}_{z}^{2}(x,t)\rangle }_{x,t}{\langle {j}_{z}^{2}(x+\delta x,t+\delta t)\rangle }_{x,t}}},$$

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,

$${{{\Delta }}}_{xt}= \frac{1}{2}\left|\frac{{\langle {j}_{z}(0,t){j}_{z}(L,t+2\nu {T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}{\sqrt{{\langle {j}_{z}^{2}(0,t)\rangle }_{t,\nu }{\langle {j}_{z}^{2}(L,t+2\nu {T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}}\right.\\ \left.-\frac{{\langle {j}_{z}(0,t){j}_{z}(L,t+(2\nu +1){T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}{\sqrt{{\langle {j}_{z}^{2}(0,t)\rangle }_{t,\nu }{\langle {j}_{z}^{2}(L,t+(2\nu +1){T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}}\right|,$$

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.

Fig. 3: Calculation on an effectively 1D single-layer long Josephson junction in BSCCO under parametric modulation at 4 K.
figure 3

The junction length is L = 25 μm and width is W = 2 μm. a Space−time trace of the current density jz(t) in the Floquet steady state under the modulation frequency \({f}_{{{{{{{{\rm{mod}}}}}}}}}=140\) GHz and modulation amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\). b Linescan along the frequency axis in (c) at the peak wavenumber, showing a peak frequency around the half-harmonic frequency \({f}_{{{{{{{{\rm{mod}}}}}}}}}/2\approx 70\) GHz. c Fourier transform to jz(x, t) in (a), plotted in the magnitude jz(β, f), where β is the wavenumber along x. The graph shows a strong half-harmonic peak in the wavenumber-frequency plane. d Linescan along the frequency axis in (c) at the peak frequency, showing a peak wavenumber corresponding to the mode index l = 9. e Space−time-averaged current−current correlation function g(δx, δt) for \({f}_{{{{{{{{\rm{mod}}}}}}}}}=140\) GHz and \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\), showing an ordered spatiotemporal pattern. f Phase diagram of the spatiotemporal order parameter Δxt versus modulation frequency \({f}_{{{{{{{{\rm{mod}}}}}}}}}\) and amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}\) (with fixed fp, temperature T = 4K, and thermal noise). Clear phase boundaries can be identified. g Spatiotemporal order parameter Δxt versus \({f}_{{{{{{{{\rm{mod}}}}}}}}}\) when \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\), showing a broadened parametric excitation band (compared with the 0D case) spanning 100−160 GHz, a range much higher than 2fp = 94 GHz. For \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\) the correction factor C = 0.95.

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 2fp, 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 arrays24,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}_{{{{{{{{\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,

$$g(\delta x,\delta z,\delta t)=\\ \frac{{\langle {j}_{z}(x,z,t){j}_{z}(x+\delta x,z+\delta z,t+\delta t)\rangle }_{x,z,t}}{\sqrt{{\langle {j}_{z}^{2}(x,z,t)\rangle }_{x,z,t}{\langle {j}_{z}^{2}(x+\delta x,z+\delta z,t+\delta t)\rangle }_{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

$${{{\Delta }}}_{xzt}= \frac{1}{2}\left|\frac{{\langle {j}_{z}(0,0,t){j}_{z}(L,H,t+2\nu {T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}{\sqrt{{\langle {j}_{z}^{2}(0,0,t)\rangle }_{t,\nu }{\langle {j}_{z}^{2}(L,H,t+2\nu {T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}}\right.\\ \left.-\frac{{\langle {j}_{z}(0,0,t){j}_{z}(L,H,t+(2\nu +1){T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}{\sqrt{{\langle {j}_{z}^{2}(0,0,t)\rangle }_{t,\nu }{\langle {j}_{z}^{2}(L,H,t+(2\nu +1){T}_{{{{{{{{\rm{mod}}}}}}}}})\rangle }_{t,\nu }}}\right|,$$

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.

Fig. 4: Calculation on an effectively 2D multi-layer long Josephson junction stack in BSCCO at 4 K under parametric modulation.
figure 4

The junction length is L = 25 μm, width is W = 2 μm, and layer number is N = 20. a Space−time trace of jz(x, 0, t) in the xt plane when z = 0 in the Floquet steady state under the modulation frequency \({f}_{{{{{{{{\rm{mod}}}}}}}}}=148\) GHz and modulation amplitude \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\). b Space−time trace of jz(0, z, t) in the xt plane when x = 0. c Fourier transform to jz(x, 0, t) plotted in the magnitude jz(β, 0, f), where β is the wavenumber along x. d Fourier transform to jz(0, z, t) plotted in the magnitude jz(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}_{{{{{{{{\rm{mod}}}}}}}}}\) when \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\) (with fixed fp, 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}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\). For \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\) the correction factor C = 0.91.

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 experiments24. 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 temperature45, 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 Tb, 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/m3K. We find, when starting at a bath temperature of 4 K and using Pac,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 jc, a 2% modulation of ns, 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 jz 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.

Fig. 5: Calculation on a 3D multi-layer long and wide Josephson junction stack in BSCCO under periodic laser illumination.
figure 5

The junction length is L = 25 μm, width is W = 5 μm, and layer number is N = 20. a An example of expected spatial profiles of \(\langle {j}_{z}^{2}(x,y,z,t)\rangle /{j_{{{\rm{c}}}}}^{2}\) in the Floquet steady state 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}_{{{{{{{{\rm{mod}}}}}}}}}=154\) GHz, the laser power is P = 0.5 mW, and the bath temperature is Tb = 4 K. b Expected Fourier spectrum of jz(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}_{{{{{{{{\rm{mod}}}}}}}}}\) at the fixed laser power 0.5 mW. This spectrum is further broadened (compared with the 2D case).

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 surface46,47 or low-temperature scanning laser microscopy48. 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.


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 Tc 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.


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 (CuO2 double layers) is ds = 0.3 nm and the thickness of the barrier layers is di = 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 Tc = 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 nth CuO2 layer or across the nth junction in the stack.

The out-of-plane current densities jz,n across the nth IJJ between the CuO2 planes n and n − 1 consist of Josephson currents with critical current density jc, (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 CuO2 plane consist of a superconducting part, characterized by a Cooper pair density ns 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 nth junction of the IJJ stack:

$$s{d}_{{{{{{{{\rm{s}}}}}}}}}\nabla \left(\frac{\nabla {\dot{\gamma }}_{n}}{{\rho }_{ab}}\right)+{\lambda }_{k}^{2}\nabla ({n}_{{{\mbox{s}}}}\nabla {\gamma }_{n})=\\ \left(2+\frac{{\lambda }_{k}^{2}{n}_{s}}{{\lambda }_{{{{{{{{\rm{c}}}}}}}}}^{2}}\right){j}_{z,n}-{j}_{z,n+1}-{j}_{z,n-1}.$$

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. jc0 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,

$${j}_{z,n}={\beta }_{{{{{{{{\rm{c0}}}}}}}}}{\ddot{\gamma }}_{n}+\frac{{\dot{\gamma }}_{n}}{{\rho }_{c,n}}+{j}_{{{{{{{{\rm{c}}}}}}}}}\sin ({\gamma }_{n})+{\eta }_{z,n},$$

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πjc0ρc0s, resistivities to ρc0 and current densities to jc0.

For the electrical parameters we use the 4 K values ρc0 = 1100 Ω cm, jc0 = 250 A/cm2, λab0 = 260 nm, and further ρab(Tc) = 100 μΩcm and ε = 13, yielding λc = 264 μm, λk = 1.52 μm and βc0 = 1.5 × 105. 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 jc 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ΓXYρc0/ρc, where Γ = 2πkBT/Ic0Φ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 jc(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 cdT/dt =  (κT) + qs, with the specific heat capacity c and the (anisotropic and layer dependent) thermal conductivity κ. The power density qs 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 Tb. 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/m3 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}_{{{{{{{{\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.

Fig. 6: Comparison between the calculated temporal correlation function for 0D and 1D Josephson junctions in BSCCO.
figure 6

The time difference δt takes up to 32768 cycles of modulation period \({T}_{{{{{{{{\rm{mod}}}}}}}}}=1/{f}_{{{{{{{{\rm{mod}}}}}}}}}\) and the modulation amplitude is \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\). a 0D case with \({f}_{{{{{{{{\rm{mod}}}}}}}}}=85\) GHz. The insets show slightly decayed oscillation amplitudes for the first and last 10 modulation periods. b 1D case (evaluated at a specific point with the space difference δx = 0) with \({f}_{{{{{{{{\rm{mod}}}}}}}}}=140\) GHz. The insets show nearly the same oscillation amplitudes for the first and last 10 modulation periods.

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 jc, i.e., \({A}_{{{{{{{{\rm{mod}}}}}}}}}=0.1\). One may consider raising the bath temperature Tb close to Tc to reduce fp 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.

Fig. 7: Experimental scheme of a parametrically modulated stack of multi-layer long and wide Josephson junctions in BSCCO.
figure 7

a Schematics of a laser-pumping induced space−time crystal in BSCCO. 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 (GL = 1).

Laser pulses with a high repetition rate can be achieved by adopting the rational harmonic mode-locking (RHML) laser technique49,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 GL = 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 fc 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 fr = (sp + 1)fc. As an estimation, a 5 m long single mode (SM) fiber loop gives fc = 40 MHz. Choosing s = 770 and p = 5, a high repetition rate of fr = 154 GHz can be obtained, while the required modulation frequency fm = 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.


  1. 1.

    Wilczek, F. Quantum time crystals. Phys. Rev. Lett. 109, 160401 (2012).

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  2. 2.

    Shapere, A. & Wilczek, F. Classical time crystals. Phys. Rev. Lett. 109, 160402 (2012).

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  3. 3.

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

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  4. 4.

    Bruno, P. Impossibility of spontaneously rotating time crystals: a no-go theorem. Phys. Rev. Lett. 111, 070402 (2013).

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  5. 5.

    Watanabe, H. & Oshikawa, M. Absence of quantum time crystals. Phys. Rev. Lett. 114, 251603 (2015).

    ADS  MathSciNet  PubMed  Article  CAS  PubMed Central  Google Scholar 

  6. 6.

    Lazarides, A., Das, A. & Moessner, R. Periodic thermodynamics of isolated quantum systems. Phys. Rev. Lett. 112, 150401 (2014).

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  7. 7.

    Sacha, K. Modeling spontaneous breaking of time-translation symmetry. Phys. Rev. A 91, 033617 (2015).

    ADS  Article  CAS  Google Scholar 

  8. 8.

    Lazarides, A., Das, A. & Moessner, R. Fate of many-body localization under periodic driving. Phys. Rev. Lett. 115, 030402 (2015).

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  9. 9.

    Khemani, V., Lazarides, A., Moessner, R. & Sondhi, S. L. Phase structure of driven quantum systems. Phys. Rev. Lett. 116, 250401 (2016).

    ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

  10. 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).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  11. 11.

    Xu, S. & Wu, C. Space−time crystal and space−time group. Phys. Rev. Lett. 120, 096401 (2018).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  12. 12.

    Sacha, K. & Zakrzewski, J. Time crystals: a review. Rep. Prog. Phys. 81, 016401 (2018).

    ADS  MathSciNet  PubMed  Article  CAS  Google Scholar 

  13. 13.

    Kozin, V. K. & Kyriienko, O. Quantum time crystals from Hamiltonians with long-range interactions. Phys. Rev. Lett. 123, 210602 (2019).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  14. 14.

    Else, D. V., Bauer, B. & Nayak, C. Floquet time crystals. Phys. Rev. Lett. 117, 090402 (2016).

    ADS  PubMed  Article  Google Scholar 

  15. 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).

    ADS  CAS  PubMed  Article  Google Scholar 

  16. 16.

    Zhang, J. et al. Observation of a discrete time crystal. Nature 543, 217 (2017).

    ADS  CAS  PubMed  Article  Google Scholar 

  17. 17.

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

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. 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).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  19. 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).

    ADS  CAS  PubMed  Article  Google Scholar 

  20. 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).

    ADS  CAS  PubMed  Article  Google Scholar 

  21. 21.

    Homann, G., Cosme, J. G. & Mathey, L. Higgs time crystal in a high-Tc superconductor. Phys. Rev. Res. 2, 043214 (2020).

    CAS  Article  Google Scholar 

  22. 22.

    Yao, N. Y., Nayak, C., Balents, L. & Zaletel, M. P. Classical discrete time crystals. Nat. Phys. 16, 438–447 (2020).

    CAS  Article  Google Scholar 

  23. 23.

    Kleiner, R., Steinmeyer, F., Kunkel, G. & Müller, P. Intrinsic Josephson effects in Bi2Sr2CaCu2O8 single crystals. Phys. Rev. Lett. 68, 2394 (1992).

    ADS  CAS  PubMed  Article  Google Scholar 

  24. 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).

    ADS  Article  CAS  Google Scholar 

  25. 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).

    ADS  Article  CAS  Google Scholar 

  26. 26.

    Feldman, M. J., Parrish, P. T. & Chiao, R. Y. Parametric amplification by unbiased Josephson junctions. J. Appl. Phys. 46, 4031 (1975).

    ADS  Article  Google Scholar 

  27. 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).

    ADS  CAS  Article  Google Scholar 

  28. 28.

    Pedersen, N. F. & Sakai, S. Prediction of half harmonic generation in stacked Josephson junctions and Bi2Sr2CaCu2Ox single crystals. Phys. Rev. B 61, 11328 (2000).

    ADS  CAS  Article  Google Scholar 

  29. 29.

    Sakai, S., Bodin, P. & Pedersen, N. F. Fluxons in thin-film superconductor−insulator superlattices. J. Appl. Phys. 73, 2411 (1993).

    ADS  Article  Google Scholar 

  30. 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).

    ADS  CAS  Article  Google Scholar 

  31. 31.

    Rudau, F. et al. Thermal and electromagnetic properties of Bi2Sr2CaCu2O8 intrinsic Josephson junction stacks studied via one-dimensional coupled sine-Gordon equations. Phys. Rev. B 91, 104513 (2015).

    ADS  Article  CAS  Google Scholar 

  32. 32.

    Rudau, F. et al. Three-dimensional simulations of the electrothermal and terahertz emission properties of Bi2Sr2CaCu2O8 intrinsic Josephson junction stacks. Phys. Rev. Appl. 5, 044017 (2016).

    ADS  Article  CAS  Google Scholar 

  33. 33.

    Tsui, O. K., Ong, N. P., Matsuda, Y., Yan, Y. F. & Peterson, J. B. Sharp magnetoabsorption resonances in the vortex state of Bi2Sr2CaCu2O8+δ. Phys. Rev. Lett. 73, 724 (1994).

    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

  34. 34.

    Matsuda, Y., Gaifullin, M. B., Kumagai, K., Kadowaki, K. & Mochiku, T. Collective Josephson plasma resonance in the vortex state of Bi2Sr2CaCu2O8+δ. Phys. Rev. Lett. 75, 4512 (1995).

    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

  35. 35.

    Ozyuzer, L. et al. Emission of coherent THz radiation from superconductors. Science 318, 1291 (2007).

    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

  36. 36.

    Welp, U., Kadowaki, K. & Kleiner, R. Superconducting emitters of THz radiation. Nat. Photonics 7, 702 (2013).

    ADS  CAS  Article  Google Scholar 

  37. 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).

    ADS  Article  CAS  Google Scholar 

  38. 38.

    Kleiner, R. Two-dimensional resonant modes in stacked Josephson junctions. Phys. Rev. B 50, 6919 (1994).

    ADS  CAS  Article  Google Scholar 

  39. 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).

    ADS  CAS  Article  Google Scholar 

  40. 40.

    Belyakov, A. O., Seyranian, A. P. & Luongo, A. Dynamics of the pendulum with periodically varying length. Physica D 238, 1589 (2009).

    ADS  MATH  Article  Google Scholar 

  41. 41.

    Glendinning, P. Adaptive resonance and pumping a swing. Eur. J. Phys. 41, 025006 (2020).

    Article  Google Scholar 

  42. 42.

    Stewart, W. C. Current−voltage characteristics of Josephson junctions. Appl. Phys. Lett. 12, 277 (1968).

    ADS  Article  Google Scholar 

  43. 43.

    McCumber, D. E. Effect of ac impedance on dc voltage−current characteristics of superconductor weak-link junctions. J. Appl. Phys. 39, 3113 (1968).

    ADS  Article  Google Scholar 

  44. 44.

    Oka, T. & Kitamura, S. Floquet engineering of quantum materials. Annu. Rev. Condens. Matter Phys. 10, 387–408 (2019).

    ADS  Article  Google Scholar 

  45. 45.

    Sobolewski, R. Ultrafast dynamics of nonequilibrium quasi-particles in hight-temperature superconductors. Proc. SPIE 3481, 480 (1998).

    ADS  CAS  Article  Google Scholar 

  46. 46.

    Stewing, F., Brendel, C. & Schilling, M. Three dimensional near-field radiation imaging up to the THz-regime. Frequenz 62, 149 (2008).

    ADS  Article  Google Scholar 

  47. 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).

    Article  Google Scholar 

  48. 48.

    Wang, H. B. et al. Hot spots and waves in Bi2Sr2CaCu2O8 intrinsic Josephson junction stacks: a study by low temperature scanning laser microscopy. Phys. Rev. Lett. 102, 017006 (2009).

    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

  49. 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).

    ADS  CAS  Article  Google Scholar 

  50. 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).

Download references


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




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

Correspondence to Reinhold Kleiner or Dafei Jin.

Ethics declarations

Competing interests

Authors declare no competing interests.

Additional information

Peer review informationNature 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

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

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).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


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.


Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing