Visualization and quantum control of light-accelerated condensates by terahertz multi-dimensional coherent spectroscopy

Abstract

Characterizing and controlling high-order correlation of quantum systems is key for developing quantum devices and switching technologies. Although conventional static and ultrafast spectroscopy gives access to collective excitations characterizing quantum states, more exotic correlations cannot be easily separated from other contributions. Here we develop density matrix simulations to show that seventh-order-wave-mixing peaks with distinct temperature and field dependences in two-dimensional terahertz nonlinear spectra reveal light-induced correlations in non-equilibrium superconducting states. Above critical terahertz driving, these emerging peaks split from conventional peaks along the second axis introduced by pump-probe relative phase in two-dimensional frequency space. They are photo-generated by correlations between two-photon fluctuations and interacting quasi-particle and quasi-particle/Higgs superconductor excitations. By photo-inducing persistent symmetry breaking via light-wave propagation, we also demonstrate seventh-order-wave-mixing sensing of Higgs collective modes. Our theory suggests to use multi-dimensional spectroscopy for quantum sensing of light-driven superconductivity and paves a path for quantum operations by few-cycle-THz-periodic photocurrent modulation.

Introduction

Recent works have established ultrafast pump-probe (PP) terahertz (THz) spectroscopy (Fig. 1a) as a powerful tool for probing and manipulating the properties of complex quantum materials. Intense THz fields can now be used to explore the competition between pairing, excitonic, spin, and lattice degrees of freedom far from equilibrium^{1,2,3,4,5,6,7,8,9,10,11,12,13}. Such light-wave fields can also break symmetries of condensed matter systems during a cycle of electric field oscillations, i.e., well before the establishment of quasi-stationary states or thermalization^{9,11,14,15,16,17}. In superconductors, quantum quench of the superconducting order parameter Δ by THz laser fields gives access to tunable long-living (100’s ps to 10’s ns) non-equilibrium states. For example, single-cycle THz pulses can drive a clean superconductor into quasi-particle quantum fluid states¹⁸. Few-cycle THz pulses give access to gapless superconducting states with broken inversion-symmetry, infinite conductivity, and nearly intact phase coherence¹⁶. Strong THz fields can also be used to induce or enhance superconductivity^1,19. However, unambiguous identification and characterization of non-equilibrium superconducting states requires quantum sensing of their collective modes and underlying order parameters. Such characterization is challenging with ultrafast PP THz spectroscopy^20,21,22,23. First, the collective excitations of superconducting states couple linearly to electromagnetic fields only in the presence of a finite Cooper-pair center-of-mass momentum p_S. Such a condensate momentum can be induced, e.g., by coexisting charge-density order^24,25, by direct supercurrent injection²⁶, or by THz-light induced breaking of inversion symmetry via electromagnetic light-wave propagation inside the superconductor^9,17. Second, the lack of experimental features uniquely attributed to collective modes has been an obstacle in the field, as several different processes contribute to the same nonlinear signals²⁷. For example, excitation of the Higgs mode, which corresponds to amplitude fluctuations of the superconducting order parameter, is only one of the processes contributing to the third-order nonlinear response of superconductors²⁸. The Higgs mode has frequency ω_H = 2Δ_∞, where Δ_∞ is the quenched order parameter of the non-equilibrium steady state. This mode is damped, since it is located at the edge of the quasi-particle continuum (Fig. 1b). Higgs-mode excitation thus competes with the excitation of quasi-particles. Depending on the strength of electron-phonon, electron-impurity, or interband Coulomb coupling in multi-band superconductors, either the Higgs-mode or the quasi-particle excitation process dominates the third-order nonlinear response^21,28,29,30.

Fig. 1: Phase-coherent nonlinear terahertz (THz) spectroscopy.

a Schematic representation of the two-dimensional terahertz configuration used here. The superconducting (SC) system is excited by collinear phase-locked pump and probe pulses. The probe electric field E_pro is centered at t_pro = 0 ps while the pump electric field E_p is centered at t_p = τ such that τ = t_p − t_pro corresponds to the pump-probe delay. ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{pp}}}}}}}}}$ is the transmitted electric field after excitation of the SC film with pump and probe pulses. b Spectra E(ω) of narrowband pump (gray shaded area) and broadband probe (solid line) laser pulses used here, with central frequency ω₀ = 1 THz. The amplitude Higgs mode frequency ω_H (quasi-particle continuum) for a pump-field strength of 160.0 kV cm^−1 is indicated by a vertical red line (purple shaded area). c, d Dynamics and spectra of the order parameter Δ for two different pump field strengths and fixed probe field strength. The order parameter spectra Δ(ω) show a second harmonic generation peak (vertical dashed line) and a Higgs mode peak at the Higgs mode frequency ω_H (vertical solid line), which separate in frequency for high pump fields (red line), but merge into one peak for lower fields (purple line, third-order susceptibility regime). e, f The corresponding dynamics and spectra of the nonlinear differential transmission ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,{\omega }_{\tau })$ at a fixed sampling time t = 0 ps. Traces are offset for clarity. The result of the calculation without pump-probe coherent modulation of the order parameter, δΔ_pp = 0 (shaded area), is close to that of the full calculation (solid line), which indicates that collective effects do not play a major role here. In addition to the two peaks of the order parameter spectrum in d, the one-dimensional pump-probe spectrum shows a third peak at ω_τ = ω₀, which results from inversion-symmetry breaking by the probe-induced condensate momentum. Inset of d: The t-dependence of the ω_τ = ω₀ inversion-symmetry breaking peak strength (black line) follows that of ${[{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}^{{{{{{{{\rm{pro}}}}}}}}}(t)]}^{2}$ (red line) where ${{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}^{{{{{{{{\rm{pro}}}}}}}}}(t)$ is probe-induced condensate center-of-mass momentum.

In semiconductors, multi-dimensional coherent nonlinear spectroscopy^31,32 at optical frequencies has been used to distinguish between different nonlinear and correlation processes^33,34,35. At THz frequencies, such spectroscopy has been applied, e.g., to study electronic excitations in semiconductors and graphene^36,37,38,39, investigate spin waves⁴⁰, and analyze vibrational modes in liquids and solids^27,41. In superconductors, two-dimensional (2D) spectroscopy measurements and the necessary theories to interpret and guide such experiments are mostly lacking. The first 2D spectroscopy experimental results have been obtained in cuprate superconductors at optical frequencies^22,42. We recently showed that purely electronic processes in iron-based superconductors with strong interband Coulomb coupling lead to a controllable hybrid-Higgs collective mode. This mode was observed as THz four-wave-mixing signal coherent oscillations, whose field and temperature dependences are clearly distinguished from those observed in single-band or weakly-coupled multi-band superconductors²¹.

In this article, we provide a comprehensive theoretical description and interpretation of the Multi-Dimensional Coherent THz (MDC-THz) spectra of superconductors, with applicability to other quantum systems. Our theory is based on gauge-invariant density matrix equations of motion in the time domain¹⁷ which are summarized in Supplementary Notes 1–4 and “Methods”. With this gauge-invariant theory, we can consistently describe superconductor amplitude, phase, and spatial dynamics driven by the phase and amplitude of THz electromagnetic field pulses. Additionally to conventional Anderson pseudo-spin precession nonlinearities, our gauge-invariant superconductor Bloch equations describe Cooper pair nonlinear quantum transport, driven by light-wave acceleration during cycles of THz electric field oscillations. Such back-and-forth acceleration leads to a moving condensate non-equilibrium state, with a finite photogenerated Cooper pair center-of-mass momentum p_S(t). This finite-momentum Cooper pairing state persists well after the THz pulse when electromagnetic propagation effects inside the superconducting system are included. By using the phase-locked pulse configuration of Fig. 1a, we control the nonlinear responses via both the phase and the amplitude of the THz fields, as well as via light-wave propagation effects dictated by Maxwell’s equations. We show that, in this way, one can use MDC-THz to disentangle the different multi-photon, quantum transport, and quantum interference nonlinear processes and characterize a THz-periodically driven moving condensate state. We focus on the non-perturbative regime of intense ultrashort excitation, where we show that susceptibility expansions are not sufficient to capture the essential signatures of a non-equilibrium light-induced superconducting state. We find that time modulation of the superconductor energy gap, controlled by the relative phase accumulation of two phase-locked pump and probe pulses, generates correlated wave-mixing spectral peaks with distinct pump field and temperature dependences. These peaks arise from seventh-order or higher wave-mixing nonlinear processes and are separated from the conventional peaks generated by previously studied PP processes. Above critical photoexcitation, they split from the known PP, four-wave mixing, and third harmonic generation (3HG) peaks along the second axis introduced by the PP relative phase in 2D frequency space. We show that such higher-order wave-mixing nonlinear peaks result from light-induced correlations between two superconductor excitations (quasi-particle pair and Higgs/quasi-particle pair excitation) and a two-photon excitation. We also identify MDC-THz peaks at equilibrium–symmetry–forbidden frequencies. These symmetry-forbidden peaks signify the Higgs collective modes characteristic of the superconducting state. These Higgs MDC-THz peaks arise from seventh-order wave mixing processes when THz dynamical symmetry breaking persists after the pulse and allow for dynamical quantum sensing of the non-equilibrium state. Finally, we show that MDC-THz can be used to both sense and coherently control PP nonlinear photogeneration of long-lived dissipationless photocurrents, which are tunable via the light-wave field phase and its electromagnetic propagation inside the superconductor^17,43.

Results

We demonstrate MDC-THz visualization of THz-light-driven moving condensate states by considering a one-band BCS superconductor. We excite the system with two phase-locked pump and probe pulses, which arrive at times t_p = τ and t_pro = 0 ps, respectively (Fig. 1a). We use a narrowband THz pump field to resolve correlated wave-mixing signals, which provide signatures of light-induced correlations in 2D frequency space. To simplify the interpretation of these signals, we consider a broadband THz probe field which is weaker than the narrowband THz pump field. The spectra of the pump and probe pulses used in the calculations are shown in Fig. 1b. Both are centered at ω₀ = 1 THz. The experimentally measured nonlinear signal, ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )={{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{pp}}}}}}}}}(t,\tau )-{{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{p}}}}}}}}}(t,\tau )-{{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{pro}}}}}}}}}(t)$, is obtained by subtracting from the transmitted E-field excited by both pump and probe, ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{pp}}}}}}}}}(t,\tau )$, the individual transmitted field transients, excited by probe (${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{pro}}}}}}}}}(t)$) and pump (${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{p}}}}}}}}}(t,\tau )$). The calculated signal therefore vanishes in the case of single-field nonlinear excitation. We study the dynamics of this correlated signal as a function of both sampling (real) time t and PP time delay τ = t_p − t_pro (the coherence time controlling the relative PP electromagnetic field phase). The corresponding MDC-THz spectra are obtained by Fourier transform of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ with respect to both t (frequency ω_t) and τ (frequency ω_τ). Similar to previous studies in semiconductors^36,37,38,39, we introduce “frequency vectors”, (ω_t, ω_τ), and “time vectors”, ${t}^{\prime}=(t,\tau )$. The pump and probe electric fields used in the calculation have the form ${{{{{{{{\bf{E}}}}}}}}}_{{{{{{{{\rm{p}}}}}}}}}({t}^{\prime})\sin ({\omega }_{{{{{{{{\rm{p}}}}}}}}}{t}^{\prime})$ and ${{{{{{{{\bf{E}}}}}}}}}_{{{{{{{{\rm{pro}}}}}}}}}({t}^{\prime})\sin ({\omega }_{{{{{{{{\rm{pro}}}}}}}}}{t}^{\prime})$, respectively, with Gaussian envelope functions ${{{{{{{{\bf{E}}}}}}}}}_{{{{{{{{\rm{p}}}}}}}}}({t}^{\prime})$ and ${{{{{{{{\bf{E}}}}}}}}}_{{{{{{{{\rm{pro}}}}}}}}}({t}^{\prime})$. The pump and probe frequency vectors are ω_p = (ω₀, −ω₀) and ω_pro = (ω₀ + Δω_pro, 0), where we use Δω_pro to denote the broadband probe frequency width. Important for obtaining the MDC-THz spectral features of main interest here is that we calculate the density matrix time evolution exactly, i.e., without any expansions.

One-dimensional coherent nonlinear spectra: role of THz dynamical symmetry breaking

Before turning to MDC-THz spectroscopy for full characterization of the non-equilibrium condensate state, we first consider the pump-induced order-parameter dynamics and the one-dimensional (1D) phase-coherent nonlinear spectrum used to observe the hybrid-Higgs mode²¹. In this way, we make the connection with previous works and also identify the effects of THz dynamical breaking of the equilibrium inversion symmetry. Figure 1c compares the time evolution of the order parameter amplitude, 2Δ(t), in response to weak (purple line) or strong (red line) THz-time-periodic fields. Oscillations at different frequencies, already visible in 2Δ(t), are characterized in Fig. 1d by the order parameter spectrum Δ(ω). The latter spectrum demonstrates second harmonic (2ω₀ = 2 THz) oscillations during the laser pulse (vertical dashed line in Fig. 1d), followed by damped Higgs mode oscillations with frequency ω_H = 2Δ_∞ < 2Δ₀ ~ 2ω₀ (vertical solid line in Fig. 1d). For weak fields, the ω_H and second harmonic resonances overlap, due to resonant Higgs mode excitation²⁸. For strong fields, the excitation energy of the superconductor is quenched non-instantaneously, during few cycles of THz-time-periodic driving. For the highest fields used here, 2Δ₀ → 2Δ_∞ is only quenched by ~20–30%. A stationary state with order parameter Δ_∞ is reached after few oscillations. This quantum quench coherent dynamics results in the low-frequency enhancement of Δ(ω) seen in Fig. 1d. On the other hand, in the extreme nonlinear regime close to complete quench of the order parameter, THz-driven Rabi flopping excites different classes of Rabi–Higgs collective modes¹⁷, which modify the MDC-THz spectra as will be shown elsewhere.

To make the connection with previous 1D phase-coherent experiments, Fig. 1e compares the dependence of the differential transmission signal ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ on the coherence time τ between weak and strong driving, for fixed real time t. Figure 1f shows the corresponding 1D spectrum, ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,{\omega }_{\tau })$, for t ≈ 0 fixed at the probe E-field maximum. This ω_τ spectrum displays two (three) distinct peaks under weak (strong) excitation. Two peaks are centered at frequencies ω_τ = 2ω₀ and ω_τ = ω_H < 2ω₀; they merge into one for weak pump excitation with ω₀ ~ Δ₀²⁸. The signal at ω_τ = 2ω₀ is generated by four-wave mixing, 2ω_p − ω_pro = (ω₀ ± Δω_pro, −2ω₀), and harmonic generation, 2ω_p + ω_pro = (3ω₀ ± Δω_pro, −2ω₀), third-order nonlinear processes, discussed in the literature before. The ω_τ = ω_H signal is generated by four-wave mixing, ω_H − ω_pro = (ω_H − ω₀ ± Δω_pro, −ω_H), 3HG, ω_H + ω_pro = (ω_H + ω₀ ± Δω_pro, −ω_H), and by the seventh-order wave-mixing processes discussed later. All of the above processes contribute at the same frequency. They are therefore hard to distinguish in 1D nonlinear spectra. Similar to previous works^20,28,44, the ω_τ = 2ω₀ and ω_H peaks in Fig. 1d predominantly arise from processes where the pump comes first and creates two-photon excitations, later sensed by the probe pulse (Supplementary Note 3). Additional processes 2ω_pro − ω_p = (ω₀ ± 2Δω_pro, ω₀), 2ω_pro + ω_p = (3ω₀ ± 2Δω_pro, −ω₀) and ω_pro − ω_pro + ω_p = (ω₀, −ω₀) arise when pump and probe pulses overlap in time. These processes contribute to the ω_τ = ω₀ peak in Fig. 1f, which is due to inversion-symmetry breaking by the probe-induced condensate center-of-mass momentum ${{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}^{{{{{{{{\rm{pro}}}}}}}}}(t=0)$ at the fixed time t ≈ 0. This peak corresponds to the symmetry-forbidden fundamental harmonic observed experimentally¹⁶. The broken inversion-symmetry non-equilibrium state can be controlled by varying the sampling time t. This is seen in the inset of Fig. 1d, which compares the t-dependence of the ω_τ = ω₀ broken symmetry peak strength (black line) to ${[{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}^{{{{{{{{\rm{pro}}}}}}}}}(t)]}^{2}$ (red line). The ω_τ = ω₀ peak follows the t-dependence of ${[{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}^{{{{{{{{\rm{pro}}}}}}}}}(t)]}^{2}$, with additional contributions for t after the probe coming from the correlated wave-mixing higher order processes discussed later. While Fig. 1e and f do not include electromagnetic propagation effects¹⁷, and thus the symmetry-breaking p_S(t) vanishes after the pulses, both Cooper-pair and amplitude Higgs mode excitation processes still contribute to ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}$, as in previous works²⁸. We isolate the Cooper-pair excitation contributions by comparing in Fig. 1e and f the full calculation of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}$ (solid lines) with the calculation for δ∣Δ_pp∣ = 0 (shaded area), where δ∣Δ_pp∣ = ∣Δ_pp∣ − ∣Δ_p∣ is the probe-induced order parameter amplitude fluctuation in the pump-driven state. Here, ∣Δ_pp∣ (∣Δ_p∣) is the order parameter amplitude that is driven by both pump and probe fields (pump field alone) (Supplementary Note 3). The results of the full calculation agree with the calculation for δ∣Δ_pp∣ = 0, which indicates that the 1D PP spectra are dominated by the light-induced Cooper-pair excitation process in one-band superconductors.

Multi-dimensional coherent nonlinear spectra

The 1D spectra for given t, discussed in the previous section and in previous works, do not provide the necessary resolution to observe higher-order correlations and corresponding nonlinear responses, since multiple different nonlinear processes contribute simultaneously at the same frequencies. For this reason, we use MDC-THz coherent spectroscopy and the 2D frequency space (ω_t, ω_τ) to separate the contributions of different nonlinear processes along the ω_τ vertical axis corresponding to the coherence time τ. Figure 2a–c shows the calculated ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ as a function of both t and τ. They demonstrate qualitative changes in the temporal profile between weak (Fig. 2a), intermediate (Fig. 2b), and strong (Fig. 2c) multi-cycle pump fields, E_p, under fixed weak probe field. In particular, while at the lowest studied field, ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ shows oscillation frequencies determined by the laser driving frequency ω₀, with increasing E_p, new interaction-dependent t-oscillation frequencies emerge, controlled by the phase difference between pump and probe fields (Fig. 2a–c). Figure 2e–g characterize these new frequency sidebands, arising from quantum interference and nonlinearity, by using the (ω_t, ω_τ)-space. We compare the calculated 2D Fourier transform spectra ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}({\omega }_{t},{\omega }_{\tau })$ for the three different E_p in Fig. 2a–c. We start with the weaker E_p = 40 kV cm^−1, where the quench of 2Δ is small in Fig. 1c. In this case, the MDC-THz signal spectral profile is similar to semiconductors³⁶ and consists of PP, four-wave mixing (4WM), and 3HG peaks described by third-order responses and Raman processes similar to previous works (Supplementary Note 3). In particular, in Fig. 2e, the peaks at (ω_t, ω_τ) = (ω₀, 0) and (ω₀, − ω₀) are generated by the PP processes ω_p + ω_pro − ω_p = (ω₀ ± Δω_pro, 0) and ω_pro + ω_p − ω_pro = (ω₀, −ω₀), respectively. The 3HG processes 2ω_p + ω_pro and 2ω_pro + ω_p produce peaks at (3ω₀ ± Δω_pro, −2ω₀) and (3ω₀ ± 2Δω_pro, −ω₀), respectively. Finally, two four-wave mixing peaks arise from 2ω_p − ω_pro and 2ω_pro − ω_p processes, at (ω₀ ± Δω_pro, −2ω₀) and (ω₀ ± 2Δω_pro, ω₀), respectively.

Fig. 2: Multi-dimensional coherent teraherz nonlinear spectroscopy.

a–c Nonlinear differential transmission ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ as a function of sampling time t (real time) and pump-probe delay τ (coherence time) for a low (40.0 kV cm^−1), intermediate (100.0 kV cm^−1), and high (160.0 kV cm^−1) pump driving field. d ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ for a calculation without pump-probe coherent modulation of the order parameter, δ∣Δ_pp∣ = 0. The result is shown for the highest pump field strength of 160.0 kV cm^−1 used here. e–g Two-dimensional (2D) Fourier transform of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ from a–c. Labels and circles indicate spectral positions of coherent pump-probe (PP, black circles), four-wave mixing (FWM, red circles), third-harmonic generation (3HG, yellow circles), and correlated high-order wave-mixing (CPP, CFWM, C3HG, and CWM, cyan circles) peaks arising from different nonlinear responses separated in 2D frequency space (ω_t, ω_τ) as discussed in the main text. Vertical magenta (black) dashed lines indicate ω_t = ω₀ and ω_t = 3ω₀ (ω_t = ω_H ± ω₀, ω_H is the Higgs mode frequency) peaks, which are split in 2D frequency space along the vertical ω_τ axis corresponding to the phase coherence time τ. Strong CPP, CFWM, C3HG, and CWM sharp peaks emerge above critical field strength and dominate the conventional peaks (PP, FWM, and 3HG). Magenta arrows indicate the 2D frequency vectors of pump (ω_p) and probe (ω_pro) pulses in e. h Two-dimensional Fourier transform of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ from d. Despite the same strong pump excitation as g, the correlated high-order wave-mixing peaks are absent for δ∣Δ_pp∣ = 0 even for the highest fields. Without δ∣Δ_pp∣, the results agree with previous results in semiconductors or the results obtained from the third-order nonlinear response (compare e and h). The distinct correlated wave-mixing peak (CWM) in the right bottom corner of the 2D space, seen in f and g, is generated only when δ∣Δ_pp∣ ≠ 0 and presents the most dramatic experimental signature of light-induced higher-order correlation in THz-time periodic-driven superconducting quantum states.

Figure 2d, h demonstrates that the above dynamics and the corresponding MDC-THz spectral profile are maintained for high E_p if we neglect the order parameter PP coherent modulation, i.e., if we set δ∣Δ_pp∣ = 0 in the gauge-invariant density matrix equations of motion. In contrast, the full calculation with δ∣Δ_pp∣ ≠ 0, Fig. 2f and g, shows a drastic change in the MDC-THz profile with increasing E_p. First, along the ω_t axis (real time), strong MDC-THz peaks (black dashed lines) split away from integer multiples of the fixed laser frequency ω₀ (magenta dashed lines), as the excitation energy ~ ω_H shifts from the second harmonic frequency 2ω₀. This shift is a result of THz-light induced quantum quench of the order parameter determining the excitation energies, an effect absent in semiconductors. More important, however, are the changes along the vertical ω_τ axis. New signals are observable when δ∣Δ_pp∣ ≠ 0, which are due to interference of pump and probe excitations that differs from previously discussed Raman processes (Supplementary Note 3). Four dominant peaks emerge with significant coherent modulation δ∣Δ_pp∣ ≠ 0 (cyan circles), all generated by seventh-order nonlinear processes not discussed before. More specifically, CPP = (ω_H − ω₀ ± Δω_pro, −ω_H + 2ω₀) is generated by ω_H − 2ω_p + ω_pro + ω_p − ω_p processes and splits from the third-order PP peak (black circle); CFWM = ω_H − ω_pro + 2ω_p − 2ω_p = (ω_H − ω₀ ± Δω, −ω_H) and C3HG = ω_H + ω_pro + 2ω_p − 2ω_p = (ω_H + ω₀ ± Δω, −ω_H) peaks also separate from third-order four-wave mixing (red circle) and 3HG (yellow circle) peaks. Importantly, a fourth strong peak, on the right bottom corner of Fig. 2f and g, is not observable for δ∣Δ_pp∣ = 0 (Fig. 2h), or for weak E_p (Fig. 2e). This correlated wave-mixing signal, CWM = 2ω_p + ω_H − (ω_pro + ω_p) + ω_p = (ω_H + ω₀ ± Δω, −ω_H − 2ω₀), at a 2D frequency where no signal is expected, is the most direct experimental evidence of light-induced correlation between two superconductor excitations (2ω_p and ω_H) and a two-photon fluctuation ω_p + ω_pro, which is controlled by the coherence time τ (Supplementary Note 3). Such coherent control by tuning τ manifests itself by a nontrivial shift of the ω_t = ω_H + ω₀ peak along the vertical axis, to ω_τ = −ω_H − 2ω₀, where no signal is expected, or observed below critical THz driving. As we discuss below, this correlated wave-mixing signal, which arises from seventh-order or higher nonlinear processes, can be used to sense light-induced superconducting states and correlations. All observed signals in the MDC-THz spectra in Fig. 2 and the corresponding nonlinear processes are summarized in Supplementary Tables 1–4.

Discussion

To understand the origin of the new correlated wave-mixing peaks characterizing the THz-light-driven moving condensate non-equilibrium state, we have derived, starting from the full gauge-invariant theory, equations of motion for pseudo-spin oscillators at each wavevector k. In this way, we generalize Anderson’s⁴⁵ original Random Phase Approximation analysis of the collective modes characterizing equilibrium superconductivity to the case of non-perturbative light-induced correlation and THz light-wave condensate acceleration (Eq. (S23) in Supplementary Note 3). Anderson pseudo-spins, with components ${\tilde{\rho }}_{i}({{{{{{{\bf{k}}}}}}}})$, i = 0 ⋯ 3, are defined by expanding the gauge-invariant density matrix in terms of the Pauli matrices, as in a standard qubit analysis (Eq. (3) in Methods). Here, up and down pseudo-spins correspond to filled and empty electronic k-states, while canted (tilted) pseudo-spins correspond to quantum superpositions of up and down states⁴⁵. As demonstrated in the Supplementary Note 3, the pseudo-spins are driven by three main coherent nonlinear processes: (i) sum-frequency Raman processes involving pump and probe excitations whose relative phase is controlled by the coherence time τ (Eq. (S24) in Supplementary Note 3, third-order response), (ii) difference-frequency Raman processes assisted by a superconductor excitation (Eq. (S25) in Supplementary Note 3, fifth-order response), and (iii) parametric driving by PP coherent modulation of the order parameter amplitude, δ∣Δ_pp(t, τ)∣ (interaction between different pseudo-spins, seventh-order response). While the first two Raman processes have been considered before, the nonlinear signal of main interest here comes from the third process, which parametrically drives seventh-order or higher responses not discussed before. The driving terms δS⁽²⁾ and δS^R in Eq. (S23) (Supplementary Note 3) are proportional to p_S(t) and the THz electric fields. They determine the third- and fifth-order responses, respectively (Supplementary Note 3). However, here we are mainly interested in the light-induced correlations driven by the source term $| {{{\Delta }}}_{{{{{{{{\rm{p}}}}}}}}}(t)| \,\delta {\tilde{\rho }}_{2}^{{{{{{{{\rm{pp}}}}}}}}}(t,\tau )$ in the equation of motion for the nonlinear differential transmission (third term in Eq. (S20)). $\delta {\tilde{\rho }}_{2}^{{{{{{{{\rm{pp}}}}}}}}}$ describes probe-induced deviations from the pump-driven state ${{\Delta }}{\tilde{\rho }}_{2}^{{{{{{{{\rm{p}}}}}}}}}$ and is excited by two phase-coherent pulses. In contrast, the nonlinear processes when the probe arrives after the pump are fully determined by the single-pulse ${{\Delta }}{\tilde{\rho }}_{2}^{{{{{{{{\rm{p}}}}}}}}}$ (first term in Eq. (S20), Supplementary Note 3). The pump-driven ${{\Delta }}{\tilde{\rho }}_{2}^{{{{{{{{\rm{p}}}}}}}}}({{{{{{{\bf{k}}}}}}}})$ describes pseudo-spin canting away from the equilibrium x–z plane orientations analogous to Anderson’s original description of the collective effects⁴⁵. On the other hand, the correlated wave-mixing peaks discussed below arise from PP fluctuations $\delta {\tilde{\rho }}_{2}^{{{{{{{{\rm{pp}}}}}}}}}(t,\tau )$ and are driven by $\sim {{\Delta }}{\tilde{\rho }}_{2}^{{{{{{{{\rm{p}}}}}}}}}({{{{{{{\bf{k}}}}}}}})\,\delta | {{{\Delta }}}_{{{{{{{{\rm{pp}}}}}}}}}(t,\tau )|^{2}$ in Eq. (S23) (Supplementary Note 3).

To clarify the origin of the proposed coherent non-perturbative signals, we make an analogy to the known four-wave mixing signals used to identify exciton–exciton interactions in semiconductors⁴⁶. There, ω_p − ω_pro second-order processes create a coherent population grating δn, which interacts with the exciton polarization P to generate an exciton–exciton interaction four-wave mixing signal driven by ~Pδn⁴⁶. Similarly, magneto–plasmon fluctuations by difference-frequency Raman processes couple to magnetoexciton coherences and generate the four-wave mixing signal of the 2D electron gas in the quantum Hall regime⁴⁷. Here, while the conventional Raman processes, δS⁽²⁾ and δS^R in Eq. (S23), determine the third- and fifth-order responses dominating in the perturbative low-intensity regime, photogenerated time-dependent interactions between different Anderson pseudo-spins, described by $\sim {{\Delta }}{\tilde{\rho }}_{2}^{{{{{{{{\rm{p}}}}}}}}}({{{{{{{\bf{k}}}}}}}})\,\delta | {{{\Delta }}}_{{{{{{{{\rm{pp}}}}}}}}}(t,\tau )|$, dominate the 2D spectra above critical driving, where the order parameter PP coherent modulation δ∣Δ_pp(t, τ)∣ becomes significant. As discussed in the Supplementary Note 3, the interaction-dependent signal discovered here is coherently driven by difference-frequency processes ω_H − ω_p − ω_pro and 2ω_p − ω_p − ω_pro, as described by Eq. (S22) in Supplementary Note 3. The correlated wave-mixing signals of main interest here thus originate from light-induced correlation between a superconductor excitation ω_H (Higgs or quasi-particle pair), a quasi-particle pair excitation 2ω_p, and a phase-coherent sum-frequency excitation ω_p + ω_pro (Supplementary Note 3). These MDC-THz peaks emerge above critical driving and are signatures of light-induced higher-order correlation.

Experimental evidence for light-induced correlation is provided by studying the pump field and temperature dependences of the correlated wave-mixing peaks in the MDC-THz spectra. These dependences are distinct from those of conventional four-wave-mixing or PP peaks and the superconductor order parameter, Δ_∞ or Δ₀. Figure 3a and b compare the ω_τ-frequency dependence of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}$ slices obtained by fixing ω_t = ω₀ (magenta dashed line in Fig. 2e–h) and ω_t = ω_H + ω₀ (black dashed line in Fig. 2e–h). They also compare the full calculation (shaded area) and that with δ∣Δ_pp∣ = 0 (solid line). The correlated wave-mixing peak, CWM, is observed in the slice ω_t = ω_H + ω₀. With increasing field, it emerges by shifting along the ω_τ-axis (dashed vertical line in Fig. 3b). This peak is completely suppressed if we set δ∣Δ_pp∣ = 0 in the equations of motion. The PP peak, on the other hand, is observed in the slice ω_t = ω₀ (vertical dashed line in Fig. 3a). In contrast to the CWM peak, it is not affected significantly when setting δ∣Δ_pp∣ = 0. Figure 3c, d compares the E_p-dependences of the PP, (ω₀, 0), and CWM, (ω_H + ω₀, −ω_H − 2ω₀), peaks. They also compare the slices from the full calculation (shaded area) with the calculation for δ∣Δ_pp∣ = 0 (solid line). Above critical field, the field-dependence of δ∣Δ_pp∣ dominates the CWM signal. As a result, the CWM field-dependence seen in the ω_t = ω_H + ω₀ slice differs distinctly from that of the PP signal seen in the ω_t = ω₀ slice. More intriguingly, Fig. 3e, f compare the temperature dependence of CWM and PP peaks up to the transition temperature T_c = 28 K, for fixed E_p = 240 kV cm^−1. For the full calculation (shaded area), the conventional PP signal mostly follows the T-dependence of the Higgs mode energy (inset Fig. 3f), since δ∣Δ_pp∣ plays a minor role (compare red and black curves). In contrast, the CWM signal follows the distinct T-dependence of the coherent modulation δ∣Δ_pp∣ of the off-diagonal coherence, which differs from that of the pump-driven order parameter ∣Δ_p∣. The T-dependence of δ∣Δ_pp∣ is also expected to reflect the coupling of the superconductor order parameter to additional degrees of freedom. As a result, correlated wave-mixing peaks in MDC-THz spectra could be used as experimental signatures of light-induced superconducting states in high-temperature superconductors.

Fig. 3: Pump field and temperature dependences of the different multi-dimensional spectral peaks.

a, b Correlated nonlinear signal ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}({\omega }_{t},{\omega }_{\tau })$ at fixed ω_t = ω₀ (magenta dashed line in Fig. 2e–h) and ω_t = ω_H + ω₀ (black dashed line in Fig. 2e–h) for the full calculation (shaded area) and the calculation without pump-probe coherent modulation of the order parameter, δ∣Δ_pp(t, τ)∣ = 0 (solid line). Vertical dashed lines indicate the pump-probe (a) and correlated wave-mixing (b) signals. c, d Pump-probe (PP) and correlated wave mixing (CWM) signals as a function of pump field strength E₀. The result of the full calculation (shaded area) is compared with the calculation where δ∣Δ_pp∣ = 0 (solid line). With increasing pump field, δ∣Δ_pp∣ ≠ 0 is responsible for the entire correlated wave mixing signal, but only slightly modifies the pump-probe signal. Inset d: Pump-field dependence of the Higgs mode energy 2Δ_∞ (non-equilibrium state order parameter). e, f Temperature dependence of the pump-probe and correlated wave-mixing peaks for the full calculation (shaded area) and the calculation with δ∣Δ_pp∣ = 0 (solid line) up to the transition temperature T_c. The pump-probe signal follows the T-dependence of the non-equilibrium order parameter 2Δ_∞ (inset f) and is only slightly affected by δ∣Δ_pp∣. In contrast, the correlated wave mixing signal is dominated by the time-dependence of δ∣Δ_pp∣ and, for this simple one-band BCS model, is suppressed well below T_c. Inset f Temperature dependence of the non-equilibrium order parameter 2Δ_∞(T). This differs distinctly from the T-dependence of the correlated wave mixing peak in f, unlike for the pump-probe peak in e.

Finally, MDC-THz spectroscopy provides direct experimental evidence for our recently proposed THz-light induced dynamical inversion symmetry breaking mechanism^16,17,43 and allows for quantum sensing of the collective mode. Figure 4a, b compares ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ with or without electromagnetic propagation effects, respectively. With propagation, the equilibrium inversion symmetry remains broken after the pump pulse, and the MDC-THz signals persist even when the pump and probe pulses do not overlap in time. This is a consequence of the predicted coherent photogeneration of a persistent nonlinear supercurrent¹⁷, which was observed experimentally^16,43. This dissipationless photocurrent forms via dynamical interplay of pump-induced nonlinear response with lightwave propagation inside the superconducting system. It persists up to many 100’s of ps due to the robustness of the condensate electronic state. The MDC-THz spectra (Fig. 4c, d) then show a series of new seventh–order wave-mixing peaks (ISWM, green circles) along the ω_τ-axis if we fix ω_t = ω_H (magenta dashed line). These inversion symmetry-forbidden peaks (Supplementary Table 5) emerge along the ω_τ axis at the Higgs mode frequency (ω_H, −ω_H) due to dynamical symmetry breaking during light-wave propagation inside the superconductor. Higgs sidebands are also seen at (ω_H, −ω_H ± ω₀) and (ω_H, −ω_H ± 2ω₀). To clarify such experimental evidence for a THz-driven non-equilibrium superconducting state with dissipationless flowing current, Fig. 4e and f compare the slices of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}$ for ω_t = ω_H (vertical magenta dashed line in Fig. 4c, d) and ω_t = ω_H + ω₀ (vertical black dashed line in Fig. 4c and d) between (i) full calculation with electromagnetic propagation effects (shaded area, inversion-symmetry breaking persists after the pump pulse), (ii) calculation without propagation (blue line, inversion symmetry only broken during the pump pulse), and (iii) calculation with propagation, but with δ∣Δ_pp∣ = 0 (red line, no order parameter coherent modulation). The CWM peak at ω_t = ω_H + ω₀ (dashed vertical line in Fig. 4f) is not much affected by the light-wave electromagnetic propagation and the persisting inversion-symmetry breaking. On the other hand, the signals at ω_t = ω_H vanish when propagation effects or δ∣Δ_pp∣ are switched off: they provide direct experimental evidence of Higgs collective modes (Supplementary Note 3) and can be used to sense non-equilibrium superconductivity.

Fig. 4: Quantum sensing of the Higgs collective mode and terahertz (THz) dynamical inversion symmetry breaking in the multi-dimensional THz coherent nonlinear spectra.

a, b Correlated differential transmission signal ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ for calculations with and without electromagnetic propagation effects leading to persistent inversion-symmetry breaking. The signal with propagation effects persists even when the pump and probe pulses do not overlap in time. c, d The corresponding two-dimensional spectra. Green circles indicate the dominant inversion-symmetry breaking signals (ISWM). Lightwave propagation generates a series of inversion-symmetry breaking wave-mixing signals along the ω_τ-axis (coherence time) for fixed ω_t = ω_H (vertical dashed magenta line) at the Higgs mode frequency in addition to the correlated wave-mixing peaks at ω_t = ω_H + ω₀ (vertical dashed black line). e, f Slices of ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}$ for ω_t = ω_H (vertical magenta dashed line in c and d) and ω_t = ω_H + ω₀ (vertical black dashed line in c and d). The full calculation with lightwave propagation (shaded area) is compared with the calculation without propagation (blue line), and the calculation with propagation but without pump-probe coherent modulation of the order parameter, δ∣Δ_pp∣ = 0 (red line). Vertical dashed lines indicate dominant inversion-symmetry breaking wave-mixing (e) and correlated wave mixing (CWM) (f) signals. The correlated high-order wave mixing peak at ω_t = ω_H + ω₀ is not affected by propagation-induced inversion symmetry-breaking effects, while the inversion-symmetry breaking wave-mixing signals at the Higgs mode frequency ω_t = ω_H vanish when propagation effects and/or δ∣Δ_pp∣ are switched off.

Methods

Gauge-invariant non-equilibrium density matrix theory

We develop quantum kinetic modeling to simulate the experimentally measurable MDC-THz spectra in quantum systems^{36,37,38,39,41}. Different nonlinear processes determine the nonlinear response to two phase-locked laser pulses and generate MDC-THz spectral peaks that are separated in 2D frequency space. In order to (i) analyze the MDC-THz spectra obtained in the non-perturbative regime of intense THz excitation, and (ii) include nonlinear quantum transport effects important for describing light-wave acceleration of the condensate^16,18,43, we use gauge-invariant density matrix equations of motion derived from the real space Bogolubov–de Gennes Hamiltonian for s-wave superconductors^17,48,49:

$$H= \mathop{\sum}\limits_{\alpha }\int {{{{{{{{\rm{d}}}}}}}}}^{3}{{{{{{{\bf{x}}}}}}}}\,{\psi }_{\alpha }^{{{{\dagger}}} }({{{{{{{\bf{x}}}}}}}})\left[\varepsilon ({{{{{{{\bf{p}}}}}}}}+e{{{{{{{\bf{A}}}}}}}}({{{{{{{\bf{x}}}}}}}},t))-\mu -e\phi ({{{{{{{\bf{x}}}}}}}},t)+{\mu }_{{{{{{{{\rm{H}}}}}}}}}({{{{{{{\bf{x}}}}}}}})+{\mu }_{{{{{{{{\rm{F}}}}}}}}}^{\alpha }({{{{{{{\bf{x}}}}}}}})\right]{\psi }_{\alpha }({{{{{{{\bf{x}}}}}}}})\\ -\int {{{{{{{{\rm{d}}}}}}}}}^{3}{{{{{{{\bf{x}}}}}}}}\left[{{\Delta }}({{{{{{{\bf{x}}}}}}}}){\psi }_{\uparrow }^{{{{\dagger}}} }({{{{{{{\bf{x}}}}}}}}){\psi }_{\downarrow }^{{{{\dagger}}} }({{{{{{{\bf{x}}}}}}}})+{{{{{{{\rm{h.c.}}}}}}}}\right].$$

(1)

Here, ${\psi }_{\alpha }^{{{{\dagger}}} }({{{{{{{\bf{x}}}}}}}})$ and ψ_α(x) are the electron creation and annihilation operators with spin index α. The superconductor complex order parameter is Δ(x) = − 2 g〈ψ_↓(x)ψ_↑(x)〉 = ∣Δ(x)∣e^iθ(x) and depends on the electron–electron pairing interaction g and the phase θ(x). The energy band dispersion is ε(p), where p = − i∇_x (ℏ=1) is the momentum operator; − e is the electron charge and μ denotes the equilibrium chemical potential. The Fock energy ${\mu }_{{{{{{{{\rm{F}}}}}}}}}^{\alpha }({{{{{{{\bf{x}}}}}}}})$ ensures charge conservation, while the Hartree energy μ_H(x) moves the phase mode of the order parameter into the quasi-particle continuum¹⁷. The driving of the superconducting system by the electromagnetic fields is described by the coupling of the vector potential A(x, t) and the scalar potential ϕ(x, t).

We numerically model the coherent dynamics driven by two THz laser pulses by using a gauge-invariant density-matrix theory¹⁷. Compared to BCS models^{45,50,51,52,53,54,55}, this gauge-invariant theory includes lightwave acceleration of condensate electrons (nonlinear quantum transport) and spatial variations. The nonlinear dynamics is described non-perturbatively by using the transformed Wigner function $\tilde{\rho }({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{R}}}}}}}})$ which, unlike for the original density matrix, is invariant under gauge transformation¹⁷:

$$\tilde{\rho }({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{R}}}}}}}})= \int {{{{{{{{\rm{d}}}}}}}}}^{3}{{{{{{{\bf{r}}}}}}}}\,\exp \left[-{{{{{{{\rm{i}}}}}}}}e\int\nolimits_{0}^{\frac{1}{2}}{{{{{{{\rm{d}}}}}}}}\lambda \,{{{{{{{\bf{A}}}}}}}}({{{{{{{\bf{R}}}}}}}}+\lambda \,{{{{{{{\bf{r}}}}}}}},t)\cdot {{{{{{{\bf{r}}}}}}}}\,{\sigma }_{3}\right]\left\langle {{{\Psi }}}^{{{{\dagger}}} }\left({{{{{{{\bf{R}}}}}}}}+\frac{{{{{{{{\bf{r}}}}}}}}}{2}\right){{\Psi }}\left.\left({{{{{{{\bf{R}}}}}}}}-\frac{{{{{{{{\bf{r}}}}}}}}}{2}\right)\right)\right\rangle \\ \times \exp \left[-{{{{{{{\rm{i}}}}}}}}e\int\nolimits_{-\frac{1}{2}}^{0}{{{{{{{\rm{d}}}}}}}}\lambda \,{{{{{{{\bf{A}}}}}}}}({{{{{{{\bf{R}}}}}}}}+\lambda \,{{{{{{{\bf{r}}}}}}}},t)\cdot {{{{{{{\bf{r}}}}}}}}\,{\sigma }_{3}\right]\,{{{{{{{{\rm{e}}}}}}}}}^{-{{{{{{{\rm{i}}}}}}}}{{{{{{{\bf{k}}}}}}}}\cdot {{{{{{{\bf{r}}}}}}}}}\,.$$

(2)

Here, ${{\Psi }}({{{{{{{\bf{x}}}}}}}})={({\psi }_{\uparrow }({{{{{{{\bf{x}}}}}}}}),{\psi }_{\downarrow }^{{{{\dagger}}} }({{{{{{{\bf{x}}}}}}}}))}^{T}$ is the field operator in Nambu space⁵⁶ and σ₃ is the Pauli spin matrix, while ${{{{{{{\bf{R}}}}}}}}=({{{{{{{\bf{x}}}}}}}}+{{{{{{{\bf{x}}}}}}}}^{\prime} )/2$ and ${{{{{{{\bf{r}}}}}}}}={{{{{{{\bf{x}}}}}}}}-{{{{{{{\bf{x}}}}}}}}^{\prime}$ are the Cooper pair center-of-mass and relative coordinates, respectively. We transform the Cooper pair relative motion to k-space, but include the spatial dependence on the Cooper pair center-of-mass, R, in order to describe the condensate motion. The above gauge-invariant density matrix is expressed in terms of Anderson pseudo-spins defined at each k-point⁴⁵ by the standard qubit expansion

$$\tilde{\rho }({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{R}}}}}}}})=\mathop{\sum }\limits_{n=0}^{3}{\tilde{\rho }}_{n}({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{R}}}}}}}}){\sigma }_{n}\,,$$

(3)

where σ_n, n = 1 ⋯ 3, are the Pauli spin matrices, σ₀ is the unit matrix, and ${\tilde{\rho }}_{n}({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{R}}}}}}}})$ are the pseudo-spin components, which depend on the center-of-mass spatial coordinate R.

We expand the full spatially-dependent equations of motion¹⁷ by assuming a weak R-dependence relative to the Cooper pair size. Below we neglect the Hartree potential for simplicity, as its effect on the nonlinear response is small for weak spatial dependence. After eliminating the order parameter phase via a gauge transformation¹⁷, we obtain gauge-invariant superconductor Bloch equations that describe a nonlinearly driven moving condensate quantum state with time-dependent center-of-mass momentum p_S(t):

$${\partial }_{t}{\tilde{\rho }}_{0}({{{{{{{\bf{k}}}}}}}})= -e{{{{{{{\bf{E}}}}}}}}\cdot {\nabla }_{{{{{{{{\bf{k}}}}}}}}}{\tilde{\rho }}_{3}({{{{{{{\bf{k}}}}}}}})+| {{\Delta }}| \left[{\tilde{\rho }}_{2}({{{{{{{\bf{k}}}}}}}}+{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)-{\tilde{\rho }}_{2}({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)\right]\,,\\ {\partial }_{t}{\tilde{\rho }}_{1}({{{{{{{\bf{k}}}}}}}})= -\left[\varepsilon ({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+\varepsilon ({{{{{{{\bf{k}}}}}}}}+{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+2\,{\mu }_{{{{{{{{\rm{eff}}}}}}}}}(t)+2{\mu }_{{{{{{{{\rm{F}}}}}}}}}(t)\right]{\tilde{\rho }}_{2}({{{{{{{\bf{k}}}}}}}})\,,\\ {\partial }_{t}{\tilde{\rho }}_{2}({{{{{{{\bf{k}}}}}}}})= \left[\varepsilon ({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+\varepsilon ({{{{{{{\bf{k}}}}}}}}+{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+2\,{\mu }_{{{{{{{{\rm{eff}}}}}}}}}(t)+2{\mu }_{{{{{{{{\rm{F}}}}}}}}}(t)\right]{\tilde{\rho }}_{1}({{{{{{{\bf{k}}}}}}}})\\ \ +| {{\Delta }}| \left[{\tilde{\rho }}_{3}({{{{{{{\bf{k}}}}}}}}+{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+{\tilde{\rho }}_{3}({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)-{\tilde{\rho }}_{0}({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+{\tilde{\rho }}_{0}({{{{{{{\bf{k}}}}}}}}+{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)\right]\,,\\ {\partial }_{t}{\tilde{\rho }}_{3}({{{{{{{\bf{k}}}}}}}})= -e{{{{{{{\bf{E}}}}}}}}\cdot {\nabla }_{{{{{{{{\bf{k}}}}}}}}}{\tilde{\rho }}_{0}({{{{{{{\bf{k}}}}}}}})-| {{\Delta }}| \left[{\tilde{\rho }}_{2}({{{{{{{\bf{k}}}}}}}}+{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)+{\tilde{\rho }}_{2}({{{{{{{\bf{k}}}}}}}}-{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}/2)\right]\,,$$

(4)

where

$$| {{\Delta }}| =-2g\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}{\tilde{\rho }}_{1}({{{{{{{\bf{k}}}}}}}})$$

(5)

is the time-dependent order parameter amplitude and θ(t) is the order parameter phase. The latter time-dependent phase determines an effective potential μ_eff(t), which also depends on the scalar potential ϕ(t). For a homogeneous system,

$${\mu }_{{{{{{{{\rm{eff}}}}}}}}}(t)=e\,\phi (t)+\frac{1}{2}\frac{\partial }{\partial t}\theta (t)-\mu ,$$

(6)

while the superfluid momentum p_S(t) is given by¹⁷

$${\partial }_{t}{{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}(t)=2\,e\,{{{{{{{\bf{E}}}}}}}}(t)\,,\quad {{{{{{{{\bf{p}}}}}}}}}_{{{{{{{{\rm{S}}}}}}}}}(t)=-2\,e\,{{{{{{{\bf{A}}}}}}}}(t).\,$$

(7)

Compared to previously studied Anderson pseudo-spin models^{45,50,51,52,53,54,55}, Eq. (4) is gauge invariant. In addition to the real-valued ∣Δ(t)∣, μ_eff(t), and p_S(t) characterizing the time-dependent driven condensate state, Eq. (4) includes quantum transport terms such as $e\,{{{{{{{\bf{E}}}}}}}}\cdot {\nabla }_{{{{{{{{\bf{k}}}}}}}}}{\tilde{\rho }}_{3}({{{{{{{\bf{k}}}}}}}})$. The latter terms are absent in previous pseudo-spin models and displace the electronic populations in k–space as a result of lightwave acceleration. The coupling between ${\tilde{\rho }}_{0}({{{{{{{\bf{k}}}}}}}})$ and ${\tilde{\rho }}_{3}({{{{{{{\bf{k}}}}}}}})$ in Eq. (4) arises from inversion-symmetry breaking induced by the lightwave field E(t) during oscillation cycles^16,18,43. For strong p_S(t) ≠ 0, the dynamical breaking of the equilibrium inversion-symmetry results in symmetry-forbidden harmonics, gapless superconductivity, and quasi-particle quantum states, which can be controlled over long time intervals of 100’s ps by tuning the ultrashort THz pulse shape^16,17,18,43.

Calculation of MDC-THz coherent spectra

The transmitted E-field used for calculating the nonlinear differential transmission ${{{{{{{{\mathcal{E}}}}}}}}}_{{{{{{{{\rm{NL}}}}}}}}}(t,\tau )$ is

$${{{{{{{\mathcal{E}}}}}}}}(t)={E}_{{{{{{{{\rm{THz}}}}}}}}}(t)-\frac{{\mu }_{0}c}{2n}J(t)\,,$$

(8)

where E_THz(t) is the applied THz electric field, n is the refractive index of the SC system, and

$$J=2e\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}{\nabla }_{{{{{{{{\bf{k}}}}}}}}}\xi ({{{{{{{\bf{k}}}}}}}}){\tilde{\rho }}_{0}({{{{{{{\bf{k}}}}}}}})$$

(9)

is the current, obtained by solving the gauge-invariant Bloch equations (4). Electromagnetic propagation inside the thin SC film is included in our calculation by self-consistently solving Eq. (8) and the gauge-invariant SC Bloch equations (4). We computed the MDC-THz spectra presented in the main text by solving Eq. (4) for equilibrium order parameter 2Δ₀ = 8.4 meV, using a square lattice nearest-neighbor tight-binding band dispersion $\varepsilon ({{{{{{{\bf{k}}}}}}}})=-2\,[{J}_{x}\cos ({k}_{x})+{J}_{y}\cos ({k}_{y})]+\mu$, with hopping parameter J_x,y = 40.0 meV and μ = 0. The coupled Bloch equations (4) and order parameter equation (5) were solved in the time domain for a 1200 × 1200 square lattice by using a fourth-order Runge–Kutta method.