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

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. At ultrafast time scales, for example, through light-matter interactions, exotic phases in quantum systems can be realised but identifying and characterising these phases from other excitations is complex. Here, the authors theoretically demonstrate how to identify the experimental signatures of a non-equilibrium superconducting state using multi-dimensional coherent spectroscopy in the terahertz regime.

R ecent 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 quasistationary 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 18 . Few-cycle THz pulses give access to gapless superconducting states with broken inversion-symmetry, infinite conductivity, and nearly intact phase coherence 16 . Strong THz fields can also be used to induce or enhance superconductivity 1,19 . However, unambiguous identification and characterization of non-equilibrium 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. E 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 ω 0 = 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 E NL ðt; ω τ Þ 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 ω τ = ω 0 , which results from inversion-symmetry breaking by the probe-induced condensate momentum. Inset of d: The t-dependence of the ω τ = ω 0 inversion-symmetry breaking peak strength (black line) follows that of ½p pro S ðtÞ 2 (red line) where p pro S ðtÞ is probe-induced condensate center-of-mass momentum.
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 26 , 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 27 . 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 28 . 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 electronphonon, electron-impurity, or interband Coulomb coupling in multi-band superconductors, either the Higgs-mode or the quasiparticle excitation process dominates the third-order nonlinear response 21,[28][29][30] .
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 40 , 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 21 .
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 17 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 lightwave 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 nonequilibrium 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 lightinduced 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 symmetryforbidden 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 ω 0 = 1 THz. The experimentally measured nonlinear signal, E NL ðt; τÞ ¼ E pp ðt; τÞ À E p ðt; τÞ À E pro ðtÞ, is obtained by subtracting from the transmitted E-field excited by both pump and probe, E pp ðt; τÞ, the individual transmitted field transients, excited by probe (E pro ðtÞ) and pump (E p ðt; τÞ). 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 E NL ðt; τÞ 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 0 ¼ ðt; τÞ. The pump and probe electric fields used in the calculation have the form E p ðt 0 Þ sinðω p t 0 Þ and E pro ðt 0 Þ sinðω pro t 0 Þ, respectively, with Gaussian envelope functions E p ðt 0 Þ and E pro ðt 0 Þ. The pump and probe frequency vectors are ω p = (ω 0 , −ω 0 ) and ω pro = (ω 0 + Δω 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 onedimensional (1D) phase-coherent nonlinear spectrum used to observe the hybrid-Higgs mode 21 . 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ω 0 = 2 THz) oscillations during the laser pulse (vertical dashed line in Fig. 1d), followed by damped Higgs mode oscillations with frequency ω H = 2Δ ∞ < 2Δ 0~2 ω 0 (vertical solid line in Fig. 1d). For weak fields, the ω H and second harmonic resonances overlap, due to resonant Higgs mode excitation 28 . 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Δ 0 → 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 lowfrequency 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 17 , 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 E NL ðt; τÞ on the coherence time τ between weak and strong driving, for fixed real time t. Figure 1f shows the corresponding 1D spectrum, E NL ðt; ω τ Þ, for t ≈ 0 fixed at the probe Efield maximum. This ω τ spectrum displays two (three) distinct peaks under weak (strong) excitation. Two peaks are centered at frequencies ω τ = 2ω 0 and ω τ = ω H < 2ω 0 ; they merge into one for weak pump excitation with ω 0~Δ0 28 . The signal at ω τ = 2ω 0 is generated by four-wave mixing, 2ω p − ω pro = (ω 0 ± Δω pro , −2ω 0 ), and harmonic generation, 2ω p + ω pro = (3ω 0 ± Δω pro , −2ω 0 ), third-order nonlinear processes, discussed in the literature before. The ω τ = ω H signal is generated by four-wave mixing, ω H − ω pro = (ω H − ω 0 ± Δω pro , −ω H ), 3HG, ω H + ω pro = (ω H + ω 0 ± Δω 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ω 0 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 = (ω 0 ± 2Δω pro , ω 0 ), 2ω pro + ω p = (3ω 0 ± 2Δω pro , −ω 0 ) and ω pro − ω pro + ω p = (ω 0 , −ω 0 ) arise when pump and probe pulses overlap in time. These processes contribute to the ω τ = ω 0 peak in Fig. 1f, which is due to inversion-symmetry breaking by the probe-induced condensate center-of-mass momentum p pro S ðt ¼ 0Þ at the fixed time t ≈ 0. This peak corresponds to the symmetry-forbidden fundamental harmonic observed experimentally 16 . The broken inversionsymmetry 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 ω τ = ω 0 broken symmetry peak strength (black line) to ½p pro S ðtÞ 2 (red line). The ω τ = ω 0 peak follows the t-dependence of ½p pro S ð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 17 , and thus the symmetrybreaking p S (t) vanishes after the pulses, both Cooper-pair and amplitude Higgs mode excitation processes still contribute to E NL , as in previous works 28 . We isolate the Cooper-pair excitation contributions by comparing in Fig. 1e and f the full calculation of E 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 E NL ðt; τÞ 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, E NL ðt; τÞ shows oscillation frequencies determined by the laser driving frequency ω 0 , with increasing E p , new interactiondependent 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 E NL ðω t ; ω τ Þ 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 36 and consists of PP, fourwave 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 , 0) and (ω 0 , − ω 0 ) are generated by the PP processes ω p + ω pro − ω p = (ω 0 ± Δω pro , 0) and ω pro + ω p − ω pro = (ω 0 , −ω 0 ), respectively. The 3HG processes 2ω p + ω pro and 2ω pro + ω p produce peaks at (3ω 0 ± Δω pro , −2ω 0 ) and (3ω 0 ± 2Δω pro , −ω 0 ), respectively. Finally, two four-wave mixing peaks arise from 2ω p − ω pro and 2ω pro − ω p processes, at (ω 0 ± Δω pro , −2ω 0 ) and (ω 0 ± 2Δω pro , ω 0 ), respectively. 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 ω 0 (magenta dashed lines), as the excitation energy~ω H shifts from the second harmonic frequency 2ω 0 . 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 − ω 0 ± Δω pro , −ω H + 2ω 0 ) 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 − ω 0 ± Δω, −ω H ) and C3HG = ω H + ω pro + 2ω p − 2ω p = (ω H + ω 0 ± Δω, −ω 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 + ω 0 ± Δω, −ω H − 2ω 0 ), 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 + ω 0 peak along the vertical axis, to ω τ = −ω H − 2ω 0 , 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 45 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ρ i ð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 45 . 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 pseudospins, 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 (2)   , 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 = ω 0 and ω t = 3ω 0 (ω t = ω H ± ω 0 , ω 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 E NL ðt; τÞ 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 lightinduced higher-order correlation in THz-time periodic-driven superconducting quantum states.
COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-022-00822-5 ARTICLE COMMUNICATIONS PHYSICS | (2022) 5:47 | https://doi.org/10.1038/s42005-022-00822-5 | www.nature.com/commsphys differential transmission (third term in Eq. (S20)). δρ pp 2 describes probe-induced deviations from the pump-driven state Δρ p 2 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 Δρ 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 46 . 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 46 . 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 47 . Here, while the conventional Raman processes, δS (2) 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 $ Δρ p 2 ðkÞ δjΔ pp ðt; τÞj, 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 quasiparticle 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 lightinduced 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 fourwave-mixing or PP peaks and the superconductor order parameter, Δ ∞ or Δ 0 . Figure 3a and b compare the ω τ -frequency dependence of E NL slices obtained by fixing ω t = ω 0 (magenta dashed line in Fig. 2e-h) and ω t = ω H + ω 0 (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 + ω 0 . 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 = ω 0 (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 pdependences of the PP, (ω 0 , 0), and CWM, (ω H + ω 0 , −ω H − 2 ω 0 ), 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 + ω 0 slice differs distinctly from that of the PP signal seen in the ω t = ω 0 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 Tdependence 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 wavemixing peaks in MDC-THz spectra could be used as experimental signatures of light-induced superconducting states in hightemperature superconductors.
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 E NL ðt; τÞ 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 17 , 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 ± ω 0 ) and (ω H , −ω H ± 2ω 0 ). 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 E NL for ω t = ω H (vertical magenta dashed line in Fig. 4c, d) and ω t = ω H + ω 0 (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 + ω 0 (dashed vertical line in Fig. 4f) is not much affected by the lightwave electromagnetic propagation and the persisting inversionsymmetry 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.

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 : Here, ψ y α ð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 μ α F ðxÞ ensures charge conservation, while the Hartree energy μ H (x) moves the phase mode of the order parameter into the quasiparticle continuum 17 . 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 17 . 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ρðk; RÞ which, unlike for the original density matrix, is invariant under gauge transformation 17 : Here, ΨðxÞ ¼ ðψ " ðxÞ; ψ y # ðxÞÞ T is the field operator in Nambu space 56 and σ 3 is the Pauli spin matrix, while R ¼ ðx þ x 0 Þ=2 and r ¼ x À x 0 are the Cooper pair centerof-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 gaugeinvariant density matrix is expressed in terms of Anderson pseudo-spins defined at Fig. 3 Pump field and temperature dependences of the different multi-dimensional spectral peaks. a, b Correlated nonlinear signal E NL ðω t ; ω τ Þ at fixed ω t = ω 0 (magenta dashed line in Fig. 2e-h) and ω t = ω H + ω 0 (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 0 . 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.
where σ n , n = 1 ⋯ 3, are the Pauli spin matrices, σ 0 is the unit matrix, andρ n ðk; 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 17 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 17 , 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): where 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), 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 E NL ðt; τÞ 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 + ω 0 (vertical dashed black line). e, f Slices of E NL for ω t = ω H (vertical magenta dashed line in c and d) and ω t = ω H + ω 0 (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 + ω 0 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.
which also depends on the scalar potential ϕ(t). For a homogeneous system, μ eff ðtÞ ¼ e ϕðtÞ þ 1 2 ∂ ∂t θðtÞ À μ; ð6Þ while the superfluid momentum p S (t) is given by 17 ∂ t p S ðtÞ ¼ 2 e EðtÞ ; p S ðtÞ ¼ À2 e AðtÞ: ð7Þ Compared to previously studied Anderson pseudo-spin models 45,50-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 E Á ∇ kρ3 ð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ρ 0 ðkÞ andρ 3 ð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 E NL ðt; τÞ is where E THz (t) is the applied THz electric field, n is the refractive index of the SC system, and 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Δ 0 = 8.4 meV, using a square lattice nearestneighbor tight-binding band dispersion εðkÞ ¼ À2 ½J x cosðk x Þ þ J y cosðk y Þ þ μ, 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.

Data availability
Source data are available for this paper. All other data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Code availability
All computer codes are available from the corresponding author upon reasonable request.