# Quantum Coherence Tomography of Lightwave–Controlled Superconductivity

L. Luo1∗, M. Mootz1,2∗, J. H. Kang3∗, C. Huang1, K. Eom3, J. W. Lee3, C.
Vaswani1,
Y. G. Collantes4, E. E. Hellstrom4, I. E. Perakis2, C. B. Eom3 and J. Wang1†
1Department of Physics and Astronomy, Iowa State University,
and Ames Laboratory, Ames, IA 50011 USA
2Department of Physics, University of Alabama at Birmingham,
3Department of Materials Science and Engineering, University of Wisconsin-
Madison,
Madison, WI 53706, USA
4Applied Superconductivity Center, National High Magnetic Field Laboratory,
Florida State University, Tallahassee, FL 32310, USA.
∗Equal contribution
†To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> Lightwave periodic driving of nearly dissipation-less currents has recently
> emerged as a universal control concept for both superconducting (SC) and
> topological electronics applications. While exciting progress has been made
> towards THz-driven superconductivity, our understanding of the interactions
> able to drive non-equilibrium pairing is still limited, partially due to the
> lack of direct measurements of high-order correlation functions. Such
> measurements would exceed conventional single-particle spectroscopies and
> perturbative responses to fully characterize quantum states
> far–from–equilibrium. Particularly, sensing of the exotic collective modes
> that would uniquely characterize lightwave-driven SC coherence, in a way
> analogous to the Meissner effect, is very challenging but much needed. Here
> we report the discovery of lightwave-controlled superconductivity via
> parametric time-periodic driving of the strongly–coupled bands in iron-based
> superconductors (FeSCs) by a unique phase–amplitude collective mode assisted
> by broken–symmetry THz supercurrents. We are able to measure
> non–perturbative, high–order correlations in this strongly–driven
> superconductivity by separating the THz multi-dimensional coherent spectra
> (THz–MDCS) into conventional pump–probe, Higgs collective mode, and
> pronounced bi–Higgs frequency sideband peaks with highly nonlinear field
> dependence. We attribute the drastic transition in the coherent spectra to
> parametric excitation of time–dependent pseudo–spin canting states modulated
> by a phase–amplitude collective mode that manifests as a strongly nonlinear
> shift from $\omega_{\mathrm{Higgs}}$ to 2$\omega_{\mathrm{Higgs}}$.
> Remarkably, the latter higher–order sidebands dominate over the lower–order
> pump–probe and Higgs mode peaks above critical field, which indicates the
> breakdown of the susceptibility perturbative expansion in the
> parametrically-driven SC state. Correlation tomography by THz–MDCS provides
> opportunities for sensing of parametric quantum matter and non-equilibrium
> SC pairing that even processes finite center-of-mass momentum, with
> implications for THz supercurrent acceleration to extend gigahertz quantum
> circuits.

Alternating “electromagnetic” bias, in contrast to DC bias, has shown promise
to enable dynamical functionalities by terahertz (THz) modulation and control
of topological-/super-currents and quantum order parameters during timescales
faster than a cycle of lightwave oscillations [9, 1, 8, 2, 3, 4, 5, 6, 7].
THz–lightwave–accelerated SC and topological currents [10, 11, 12, 13, 14, 8,
15, 16, 9, 17] have revealed exotic quantum dynamics, e. g., harmonic modes
[9, 18, 10] and gapless quantum fluid states [19] forbidden by equilibrium SC
pairing symmetry, or light-induced Weyl and Dirac nodes [14, 1]. However,
high–order correlation characteristics far exceeding the known two-photon
light coupling to superconductors are hidden in conventional spectroscopy
signals and perturbative responses, where a mixture of multiple excitation
pathways contribute to the same low-order responses [20, 21]. A compelling
solution to sensing SC coherence of light–driven states is to be able to
identify unambiguously their collective modes [8, 22, 23, 24, 25, 26, 27]. The
dominant collective excitations of the equilibrium SC phase range from
amplitude fluctuations (Higgs mode) to fluctuations of interband phase
differences (Leggett mode) of the SC order parameter. Although amplitude modes
have been observed close to equilibrium when external DC [28, 17, 29] and AC
[9, 18] fields break inversion symmetry (IS), here we show that the phase
coherent dynamics of the SC order parameters can parametrically drive quantum
states, yet-to-be-observed, accessed by strong THz coherent two–pulse
excitations. These states are characterized by distinct phase-amplitude
collective modes arising from strong light-induced couplings between the
amplitude and phase channels.

Fig. 1. Terahertz (THz) multidimensional coherent spectroscopy of lightwave
accelerated non–equilibrium superfluid states in FeSCs. $\bf{a}$, Schematics
of lightwave supercurrent generation, coherent control and detection of the
parametrically-driven SC state via two phase-locked THz pulses in our
experiment. $\bf{b}$, Temporal waveforms of the nearly single-cycle THz
pulse–pair used in the experiment (red and blue lines), and $\bf{c}$, spectra
of the used pulses, centered at $\omega_{0}=4$ meV. $\bf{d}$, Temporal
dynamics of the measured coherent nonlinear transmission
$E_{\mathrm{NL}}(t,\tau)(\mathrm{pink})=E_{\mathrm{AB}}(t,\tau)(\mathrm{black})-E_{\mathrm{A}}(t)(\mathrm{red})-E_{\mathrm{B}}(t,\tau)(\mathrm{blue})$
as a function of gate time $t$ at a fixed delay time between the two pulses,
$\tau=6.5$ ps, under THz driving fields of 229 kV/cm at temperature of 5 K.
$\bf{e}$, Temporal dynamics of the $E_{\mathrm{NL}}(t,\tau)$ amplitude decay
below (red diamond, 5 K) and above (black cross, 40 K) $T_{\mathrm{c}}$, as a
function of pulse–pair time delay $\tau$ under THz driving fields of 333
kV/cm. The correlated nonlinear signal $E_{\mathrm{NL}}(t,\tau)(\mathrm{red})$
decays over timescales much longer than the pulse duration.

THz frequency, multi-dimensional coherent nonlinear spectroscopy (THz-MDCS) )
[11, 30, 31, 32, 33, 34, 35] represents a correlation tomography tool to
distinguish between different many-body response functions and light–induced
collective modes in superconductors under strong two–pulse THz excitation.
Unlike for THz–MDCS studies of semiconductors [11, 31, 32, 36], magnets [33],
and molecular crystals [34], Fig. 1a illustrates three distinct features of
our scheme in superconductors, not explored so far. First, our approach is
based on measuring the phase of the supercurrent coherent nonlinear emission,
in addition to the amplitude, by using phase-resolved coherent measurements
with two intense phase–locked THz pulses of similar field strengths that has
not been applied on superconductors. Taking advantage of both the real time
and the relative phase of the two THz fields, we separate in two-dimensional
(2D) frequency space spectral peaks generated by light-induced correlations
and collective mode interactions from the conventional pump–probe,
four–wave–mixing, and high–harmonic–generation signals [9, 23]. This 2D
separation of spectral peaks arising from high–order nonlinear processes
achieves a “super” resolution of many-body interactions and collective modes
in highly non-perturbative states, which is not possible with the conventional
spectroscopy one-dimensional measurements of prior works [23, 9, 19]. Second,
as a result of lightwave condensate acceleration by the effective local field
inside a thin–film SC induced by two sub-gap THz pulses and electromagnetic
propagation effects, the Cooper pairs ($\mathbf{k}$, $\mathbf{-k}$) of the
equilibrium BCS state experience a highly nonlinear center–of–mass momentum,
$\mathbf{p}_{\mathrm{S}}(t)$, i. e., SC pairing with finite center-of-mass
momentum (Fig. 1a). Precisely, the state persists well after the two strong
pulses and results in
($\mathbf{k}+\mathbf{p}_{\mathrm{S}}(t)/2,\mathbf{-k}+\mathbf{p}_{\mathrm{S}}(t)/2$)
Cooper pairing, due to dynamical symmetry breaking of the centrosymmetric
pairing states [9]. Third, the driven quantum state with current-flow
$\propto\mathbf{p}_{\mathrm{S}}(t)$ controllable by two–pulse interference can
host distinct collective modes that provide parametric excitation of the time-
dependent quantum states (Fig. 1a), whose nonlinear interactions determine the
THz–MDCS spectral profile [37].

In this Article, we reveal a superconducting state parametrically driven by
time–periodic light–induced dynamics of the order parameter phase in a
Ba(Fe1-xCox)2As2 superconductor. Such parametric driving becomes important
when the phase dynamics is amplified by a phase–amplitude collective mode that
develops with increasing THz pulse–pair driving. The change in the character
of the SC state manifests itself via the drastic changes in the THz–MDCS
spectra observed with increased two–pulse driving. In particular, we observe a
transition in THz–MDCS, from wave–mixing peaks centered at the Higgs amplitude
mode frequency $\omega_{\mathrm{H}}$ to bi–Higgs frequency sidebands at
2$\omega_{\mathrm{H}}$. Remarkably, the higher–order 2$\omega_{\mathrm{H}}$
peaks dominate over the lower order single–Higgs peaks and conventional
pump–probe spectra above critical THz field strength. We attribute the
non–perturbative sidebands at bi–Higgs frequencies to the parametrically
driven SC state by the collective temporal oscillations in the phases of the
$s_{\pm}$–symmetry order parameter, consistent with our quantum kinetic
simulations.

We measured optimally Co-doped BaFe2As2 (Ba-122) epitaxial thin film (60 nm)
with $T_{\mathrm{c}}\sim$23 K and lower SC gap 2$\Delta_{1}\sim$6.8 meV
(Methods section 1.1). We used THz-MDCS to measure the responses to two phase-
locked, nearly single-cycle THz pulses A and B of similar field strength (Fig.
1b), with central frequency $\omega_{0}\sim$4 meV (black arrow, Figs. 1c) and
broadband frequency width of $\Delta\omega\sim$6 meV (purple dashed line, Fig.
1c) (Methods section 1.2). Representative time scans of these THz-MDCS
experiments driven by laser fields $E_{\mathrm{THz,A,B}}=$229 kV/cm, are shown
in Fig. 1d. The measured nonlinear differential emission correlated signal,
$E_{\mathrm{NL}}(t,\tau)=E_{\mathrm{AB}}(t,\tau)-E_{\mathrm{A}}(t)-E_{\mathrm{B}}(t,\tau)$,
was recorded as a function of both the gate time $t$ (Fig. 1d) and the delay
time $\tau$ between the two pulses A and B (Fig. 1e). We note three points.
First, as demonstrated by $E_{\mathrm{NL}}(t,\tau)$ shown in Fig. 1d (pink
cross), measured at fixed delay $\tau=6.5$ ps, the electric field in the time
domain allows for simultaneous amplitude-/phase-resolved detection of the
coherent nonlinear responses induced by the pulse–pair and has negligible
contributions from the individual pulses. This is achieved by subtracting the
individual responses, $E_{\mathrm{A}}(t)$ and $E_{\mathrm{B}}(t,\tau)$ (red
and blue solid lines), from the full signal obtained in response to both
phase–locked pulses, $E_{\mathrm{AB}}(t,\tau)$ (black solid line). Second, the
$E_{\mathrm{NL}}(t,\tau)$ in Fig. 1e vanishes above the SC transition
temperature $T_{\text{c}}$, as seen by comparing the 5 K (red diamond) and 40
K traces (black cross). Third, the THz-MDCS signals persist even when the two
pulses do not overlap in time, e. g., at $\tau=6.5$ ps (Figs. 1d–1e). The
long–lived correlated signal $E_{\mathrm{NL}}(t,\tau)$ indicates that the two
sub-gap laser excitations, centered below 2$\Delta_{1}$ (Fig. 1c), have
generated robust supercurrent–carrying macroscopic states persisting well
after the pulse.

Fig. 2. Drastic changes of correlations and collective modes revealed in the
driving electric field dependence of THz-MDCS. $\bf{a}$-$\bf{c}$, Two-
dimensional (2D) false-colour plot of the measured coherent nonlinear
transmission $E_{\mathrm{NL}}(t,\tau)$ of FeSCs superconductors at 5 K induced
by THz pump electric fields of ($\bf{a}$) 229 kV/cm, ($\bf{b}$) 333 kV/cm and
($\bf{c}$) 475 kV/cm. $\bf{d}$-$\bf{f}$, The corresponding THz 2D coherent
spectra $E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ at 5 K for the above three
pump electric fields, respectively. $\bf{g}$-$\bf{i}$, The normalized
$E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ spectra are plotted for the same
pump fields to highlight the pump field-dependent evolution of the correlation
peaks along the 2D frequency vector space. Peaks marked by the dashed lines
are located at frequencies associated with the Higgs (green) mode and bi-Higgs
frequency sideband (red and blue) consistent with the theory. $\bf{j}$, An
example of calculated $E_{\mathrm{NL}}(t,\tau)$ as a function of gate time $t$
and delay time $\tau$ for 250 kV/cm pump field. $\bf{k}$-$\bf{n}$, 2D Fourier
transform of $E_{\mathrm{NL}}(t,\tau)$. As predicted by the theory, dashed
black (blue) lines indicate pump-probe
$\omega_{t}=\omega_{\mathrm{H,1}}-\omega_{0}$ (bi-Higgs frequency sideband
$\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}$) while IS-breaking signals at
Higgs $\omega_{t}=\omega_{\mathrm{H,1}}$ (bi-Higgs frequency sideband
$\omega_{t}=2\omega_{\mathrm{H,1}}-2\omega_{0}$) are marked by vertical dashed
green (red) line; pump–probe peaks at $\omega_{t}=\omega_{0}$ are indicated by
vertical dashed magenta lines. $\bf{o}$, Field-strength dependence of the
dominant peak in the spectrum of the interband phase difference
$\delta\theta(\omega)$. $\bf{p}$, Field-strength dependence of the bi-Higgs
frequency sidebands at $2\omega_{\mathrm{H,1}}-\omega_{0}$ follows the
$\delta\theta(\omega)$ behavior in (o), which identifies the importance of
light–induced time–periodic phase dynamics at the $\omega_{\mathrm{H,1}}$
frequency in driving a non–equilibrium SC state. Three excitation regimes are
marked (main text).

Figure 2 compares the 2D THz temporal profile of the coherent nonlinear signal
$E_{\mathrm{NL}}(t,\tau)$ for relatively weak (Fig. 2a), intermediate (Fig.
2b), and strong (Fig. 2c) driving fields. The $E_{\mathrm{NL}}(t,\tau)$
dynamics reveals that pronounced coherent temporal oscillations last much
longer than the temporal overlap between the two driving pulses (Fig. 1b).
These long–lived coherent responses generate sharp THz-MDCS spectral peaks
visible up to $\sim$8 meV below substrate absorption (Methods). These spectra
were obtained by Fourier transform of $E_{\mathrm{NL}}(t,\tau)$ with respect
to both $t$ (frequency $\omega_{t}$) and $\tau$ (frequency $\omega_{\tau}$)
(Figs. 2d–2f). We observe multiple distinguishing and well-defined resonances
with unique lineshapes that drastically change with increasing field strength.
These $E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ spectra differ strongly from
the conventional ones measured, e. g., in semiconductors [31, 32, 11]. In the
latter uncorrelated systems, peaks are observable at multiples of the THz
driving pulse frequency $\omega_{0}\sim$4 meV (magenta dashed line), as
expected in the case of a rigid excitation energy bandgap. The observed peaks
in FeSCs are much narrower than the excitation pulse width $\Delta\omega$
(Fig. 1c). This result implies that $E_{\mathrm{NL}}(t,\tau)$ oscillates with
the frequencies of SC collective mode excitations that lie within the
$\Delta\omega$ of the few-cycle driving pulses, with the width of the THz–MDCS
spectral peaks determined by the SC mode damping and not by $\Delta\omega$ of
the driving pulses.

The normalized $E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ experimental
spectra shown in Figs. 2g–2i visualize nonlinear couplings of SC collective
mode resonances and their field-dependences. For the weaker pump field of
$E_{0}=229$ kV/cm in Fig. 2g, the THz–MDCS spectrum shows four dominant peaks.
Intriguingly, the two strongest peaks are located at the higher frequencies,
roughly $(6,0)$ meV and $(6,-6)$ meV, with the weaker peaks at the lower
frequencies, slightly below $(2,0)$ meV and $(2,-2)$ meV. This observation is
in strong contrast to the expectation from conventional harmonic generation
that high–order nonlinear signals should be weaker than lower–order ones. Such
a reversal of coherent nonlinear signal strengths indicates a breakdown of a
susceptibility perturbative expansion around the SC equilibrium state. For the
intermediate field of $E_{0}=333$ kV/cm (Fig. 2h), the THz-MDCS spectrum shows
several peaks close to each other (red and green lines), centered at new
frequencies $\sim(5,0)$ meV and $(5,-5)$ meV which exhibit the similar non-
perturbative behavior with dominant high order THz-MDCS spectral peaks. The
spectral profile changes again with increasing THz driving: four peaks are
observable in the THz-MDCS spectrum for the highest studied pump field of
$E_{0}=475$ kV/cm. For such high fields, the two strongest THz–MDCS peaks are
roughly located at $(2.3,0)$ meV and $(2.3,-2.3)$ meV, while two weaker peaks
become detectable at higher frequencies, at $(6.2,0)$ meV and $(6.2,-6.2)$ meV
(Fig. 2i). These high field peaks should be distinguished from the low field
ones at similar frequencies (Fig. 2g), as the latter have red-shifted with
increasing field due to the SC gap reduction. The evolution of the well-
defined MDCS spectral peaks reflect the emergence of different collective
modes with increasing driving field, which characterize the transition to
different non-equilibrium SC states.

We use three principles to classify the observed peaks in
$(\omega_{t},\omega_{\tau})$ space. First, we introduce frequency vectors
characterizing the two pulses A and B,
$\omega_{\mathrm{A}}=(\omega_{0}\pm\Delta\omega,0)$ and
$\omega_{\mathrm{B}}=(\omega_{0}\pm\Delta\omega,-\omega_{0}\mp\Delta\omega)$,
which are centered around $\omega_{0}\sim$4 meV (black arrow, Fig. 1c). The
corresponding “time vectors” $t^{\prime}=(t,\tau)$ allow us to represent the
driving electric fields of Fig. 1b in the form ${\bf
E}_{\mathrm{A}}(t^{\prime})\sin(\omega_{\mathrm{A}}t^{\prime})$ and ${\bf
E}_{\mathrm{B}}(t^{\prime})\sin(\omega_{\mathrm{B}}t^{\prime})$, respectively.
Second, the non-equilibrium SC state driven by the above pulse–pair is
characterized by a quenched asymptotic value of the time–evolved SC order
parameter, which defines the Higgs frequencies
$\omega_{\mathrm{H},i}=2\Delta_{\infty,i}$, where $i=1$ ($i=2$) denotes the
hole (electron) pocket of the FeAs bandstructure. The above Higgs mode
frequencies decrease from their equilibrium values of $2\Delta_{0,i}$ with
increasing field, which leads to a redshift of the THz–MDCS spectral features
observed in Figs. 2g–2i. Note that we only probe the lower Higgs mode,
$\omega_{\mathrm{H},1}\sim 6.8$ meV, while the higher Higgs frequency,
$\omega_{\mathrm{H},2}\sim 19$ meV, lies outside of the measured spectral
range (Fig. 1c). Third, the THz pulses drive the Anderson pseudospin
oscillators [9, 35] at the different momenta ${\bf k}$ (Methods section 1.4),
whose dynamics is dominated by frequencies
$\sim\omega_{\mathrm{H,i;A}}=(\omega_{\mathrm{H,i}},0)$ and
$\sim\omega_{\mathrm{H,i;B}}=(\omega_{\mathrm{H,i}},-\omega_{\mathrm{H,i}})$,
i. e., field-dependent Higgs and quasi-particle pair excitations, or $\sim
2\omega_{\mathrm{A,B}}$, i. e., quasi-particle excitations driven at the laser
frequency.

To identify which nonlinear process generates each peak measured in Figs.
2g–2i, we use our quantum kinetic simulations (Methods section 1.3) and the
above three principles. Light–wave propagation inside a SC thin film geometry
determines the effective driving field
$E(t)=E_{\mathrm{THz}}(t)-\frac{\mu_{0}c}{2n}J(t)$, which is obtained from
Maxwell’s equations [37] and differs from the applied field
$E_{\mathrm{THz}}(t)$ ($n$ is the refractive index). This effective field
drives the nonlinear supercurrent $J(t)$, described self-consistently by
solving the gauge-invariant SC Bloch equations [37, 23, 35] (Methods section
1.3) for a 3-pocket SC model with strong electron–hole pocket interaction $U$
far exceeding the intra–band pairing interaction. Using the above results, we
simulate directly the $E_{\mathrm{NL}}(t,\tau)$ temporal dynamics measured in
the experiment (Fig. 2j as an example) and then obtain the
$E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ spectra (Figs. 2k–2n). These
simulations are fully consistent with the observed drastic change in the
THz–MDCS spectra, where non-perturbative spectral peaks emerging with
increasing field, as shown in Figs. 2g–2i, are indicative of a transition to
light-driven SC states with different, emergent collective modes.

We elaborate the above quantum state transition by using three different
excitation regimes, marked in Figs. 2o (red arrows) as field-strength
dependence of the interband phase difference $\delta\theta(\omega)$ peak: I,
perturbative regime; II, the non-perturbative state with dominant Higgs
amplitude mode; regime III, the parametrically-driven SC state determined by
phase–amplitude collective mode. We first examine the perturbative
susceptibility regime I, where the Higgs frequency $\omega_{\text{H,1}}$
remains close to its equilibrium value, $2\Delta_{1}\sim 6.8$ meV, similar to
the “rigid” excitation energy gap in semiconductors. The simulated THz-MDCS
spectrum (Fig. 2k) then shows several peaks (Table 1, Methods) splitting along
the $\omega_{\tau}$ vertical axis, at $\omega_{t}=\omega_{0}$ (dashed magenta
line) and $\omega_{t}=\omega_{\text{H,1}}$ (dashed green line). The
conventional pump–probe signals are observed at $(\omega_{0},-\omega_{0})$ and
$(\omega_{0},0)$ in Fig. 2k, generated by the familiar third-order processes
$\omega_{\mathrm{A}}-\omega_{\mathrm{A}}+\omega_{\mathrm{B}}$ and
$\omega_{\mathrm{B}}-\omega_{\mathrm{B}}+\omega_{\mathrm{A}}$, respectively.
Four-wave mixing signals are also observed at $(\omega_{0},\omega_{0})$ and
$(\omega_{0},-2\omega_{0})$, generated by the third–order processes
$2\omega_{\mathrm{A}}-\omega_{\mathrm{B}}$ and
$2\omega_{\mathrm{B}}-\omega_{\mathrm{A}}$. However, the perturbative behavior
in this regime are inconsistent with the dominance of higher–order peaks (Fig.
2g) for the much stronger fields used in the experiment to achieve the
necessary signal–to–noise ratio.

By increasing the field strength (Figs. 2l–2n), the calculated signals along
the $\omega_{\tau}$ vertical axis and at $(\omega_{0},-\omega_{0})$,
$(\omega_{0},0)$ diminish. Only peaks along $(\omega_{t},0)$ and
$(\omega_{t},-\omega_{t})$ are then predicted by our calculation, consistent
with the experiment in Figs. 2g–2i. For the lower field strength of 250 kV/cm
(Fig. 2l), our calculated THz-MDCS spectrum shows two weak peaks at
$\omega_{t}\sim 2$ meV (black dashed line) and two strong broken–IS peaks at
$\omega_{t}=\omega_{\mathrm{H,1}}\sim 6$ meV (green dashed line), similar to
the experimental THz-MDCS peaks shown in Fig. 2g. The weak peaks at
$\omega_{t}\sim 2$ meV (black dashed line) arise from high-order difference-
frequency Raman processes (PP, Table 2 in Methods), which generate pump–probe
signals at $\omega_{t}=\omega_{\mathrm{H,1}}-\omega_{0}$, as observed in Figs.
2d and 2g. The strong peaks at the Higgs frequency
$\omega_{t}=\omega_{\mathrm{H,1}}\sim$6 meV (green dashed line) dominate for
intermediate fields up to $\sim$400 kV/cm (regime II, Fig. 2o), but vanish if
we neglect the electromagnetic propagation effects as discussed later. The BCS
ground state evolves into a finite–momentum–pairing SC state, which is
determined by the condensate momentum $\mathbf{p}_{\mathrm{S}}$ generated by
third–order nonlinear processes (Supplementary Fig. 4c, Note 4) and persisting
well after the pulse. Higgs frequency peaks then arise from ninth–order IS
breaking nonlinear processes generated by the coupling between the Higgs mode
and the lightwave accelerated supercurrent $J(t)$ (IS Higgs, Table 2 in
Methods). The superior resolution achieved for sensing the collective modes
using THz–MDCS with 2D coherent excitation is far more than a static IS
symmetry breaking scheme using a DC current (Supplementary Fig. 9, Note 7).

Fig. 3. Temperature dependence of THz-MDCS signals. $\bf{a}$, THz-MDCS spectra
$E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ at 16 K for pump electric field
333 kV/cm. $\bf{b}$, Temporal profiles of two-pulse THz coherent signals
$E_{\mathrm{NL}}(t,\tau)$ at various temperatures from 5 K to 22 K for a peak
THz pump electric field of $E_{\mathrm{pump}}=229$ kV/cm and $\tau=6.5$ ps.
Traces are offset for clarity. $\bf{c}$, The corresponding Fourier spectra of
the coherent dynamics in ($\bf{b}$). $\bf{d}$-$\bf{e}$, A 2D false-color plot
of THz coherent signals ($\bf{e}$) as a function of temperature and frequency
$\omega_{t}$ with ($\bf{d}$) integrated spectral weight at various
temperatures. Dashed gray line indicates the SC transition temperature.

For an even higher field strength of 350 kV/cm (Fig. 2m) and of 700 kV/cm
(Fig. 2n), the THz–MDCS spectra change above the excitation threshold where
the order parameter phase dynamics becomes significant, as shown in Fig. 2o
(regime III). As discussed in Methods (section 1.4), time–periodic phase
deviations from their equilibrium values, with frequency
$\omega_{\mathrm{H,1}}$, parametrically drive pseudo–spin time–dependent
canting (Supplementary Fig. 4d), which coincides with a highly nonlinear
field–dependence of $\mathbf{p}_{\mathrm{S}}$ (Supplementary Fig. 4c). In this
regime III, new dominant THz–MDCS peaks emerge at
$\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}$ (blue dashed line), referred as
to bi-Higgs frequency sideband, while satellite peaks are also observed at
$\omega_{t}=2\omega_{\mathrm{H,1}}-2\omega_{0}$ (red dashed line). The
spectral position of these emergent sideband peaks in 2D frequency space
indicates that they arise from a second harmonic of the Higgs frequency,
$2\omega_{\mathrm{H,1}}$. Figure 2p demonstrates a threshold nonlinear
behavior of these bi-Higgs frequency sideband peak strengths, which coincides
with the development of strong phase dynamics as seen in Fig. 2o.

These theoretical prediction are fully consistent with our experimental
observations in Figs. 2h and 2i. For the intermediate field in Fig. 2m, the
THz–MDCS peaks at $\omega_{t}=2\omega_{\mathrm{H,1}}-2\omega_{0}\sim$4 meV
(red dashed line) and $\omega_{t}=\omega_{\mathrm{H,1}}\sim$6 meV (green
dashed line) are close to each other, so they merge into a single broad
resonance around $(5,0)$ meV and $(5,-5)$ meV in the calculated 2D spectrum.
This result agrees with the measured broad, overlapping THz–MDCS peaks $\sim$5
meV in Fig. 2h, while the calculated
$\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}$ peak (blue line) is not visible
experimentally due to the substrate absorption. For the highest studied field
strength (Fig. 2n), the calculated THz–MDCS signals are dominated by the
bi–Higgs frequency nonlinear sidebands at
$\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}\sim 6.0$ meV and
$\omega_{t}=2\omega_{\mathrm{H,1}}-2\omega_{0}\sim 2.0$ meV. Both sidebands
peaks now fall into the substrate transparency region and are clearly resolved
in Fig. 2i. The emergence of these new MDCS peaks in Regime III will be
discussed later as manifestations of a phase–amplitude collective mode that
parametrically drive pseudo-spin canting with respect to the
$s_{\pm}$–symmetry equilibrium directions.

Figure 3 demonstrates the strong temperature dependence and redshift of the
observed peaks as we appoach $T_{\mathrm{c}}$. The THz-MDCS spectrum
$E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ at temperature 16 K is shown in
Fig. 3a for the intermediate field of $E_{0}=333$ kV/cm, and is compared in
Fig. 2h with the spectrum at $T$=5 K for same excitation. The broken–IS
signals observed at the Higgs mode frequency
$\omega_{t}=\omega_{\mathrm{H},1}$ red-shift with increasing temperature, from
$(5,0)$ meV and $(5,-5)$ meV peaks at 5 K to broad peaks slightly below
$(2.5,-2.5)$ meV and $(2.5,0)$ meV at 16 K (green line). This redshift arises
from the thermal quench of the SC order parameter $2\Delta_{1}$ with
increasing temperature. However, unlike for the case of THz coherent control
of the order parameter demonstrated in Fig. 2h, a thermal quench does not
produce any obvious bi–Higgs frequency THz–MDCS peaks expected $\sim$1 meV,
which is indicative of the coherent origin of the latter. Figures 3b and 3c
show the temperature dependence of the measured differential coherent emission
$E_{\mathrm{NL}}(t,\tau)$ and the corresponding
$E_{\mathrm{NL}}(\omega_{t},\tau)$ at a fixed pulse separation $\tau=6.5$ ps.
It is clearly seen, by comparing the 5 K (black line) and 22 K (gray) traces,
that, when approaching $T_{\text{c}}$ from below, the coherent nonlinear
emissions quickly diminish and red-shift. Finally, Figs. 3d–3e show a detailed
plot of $E_{\mathrm{NL}}(\omega_{t},\tau)$ up to 100 K. The integrated
spectral weight shown in Fig. 3d correlates with the SC transition at
$T_{\text{c}}$ (gray dashed line).

Fig. 4. Origin of correlation and collective mode peaks in THz-MDCS signals
underpinned by switch-off analysis. $\bf{a}$-$\bf{b}$,,
$E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ for ($\bf{a}$) the full
calculation with lightwave propagation ($E_{0}=250$ kV/cm) and for ($\bf{b}$)
a calculation without propagation ($E_{0}=250$ kV/cm). Dashed black (green)
lines indicate $\omega_{t}=\omega_{\mathrm{H,1}}-\omega_{0}$
($\omega_{t}=\omega_{\mathrm{H,1}}$). The IS peaks at
$\omega_{t}=\omega_{\mathrm{H,1}}$ vanish without persisting IS breaking.
$\bf{c}$-$\bf{e}$, $E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ for ($\bf{c}$)
the full calculation with lightwave propagation ($E_{0}=700$ kV/cm), for
($\bf{d}$) a calculation without phase–amplitude coupling ($E_{0}=400$ kV/cm),
and for ($\bf{e}$) a calculation without interband interaction, $U=0$
($E_{0}=100$ kV/cm). To directly compare the different THz-MDCS spectra, the
field strengths of the different calculations are chosen such that
$\omega_{\mathrm{H,1}}$ are comparable, $\omega_{\mathrm{H,1}}\sim 5.0$ meV.
Dashed blue lines indicate $\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}$,
while IS-breaking signals at Higgs $\omega_{t}=\omega_{\mathrm{H,1}}$ (bi-
Higgs $\omega_{t}=2\omega_{\mathrm{H,1}}-2\omega_{0}$) are marked by vertical
dashed green (red) lines. Note that the bi-Higgs frequency sideband peaks are
strongly suppressed without phase–amplitude coupling or without interband
interaction.

Figure 4 offers more insight into the physical mechanism behind the observed
transition in the THz–MDCS spectra with increasing field. First, we compare
the spectra $E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ for a field strength
of $250$ kV/cm between (i) the full calculation that includes electromagnetic
propagation and interference effects leading to slowly decaying
$\mathbf{p}_{\mathrm{S}}(t)$ after the pulse (Fig. 4a), and (ii) a calculation
where these effects are switched off, in which case
$\mathbf{p}_{\mathrm{S}}(t)$ oscillates during the THz pulse and vanishes
afterwards (Fig. 4b). In the latter case, the
$\omega_{t}=\omega_{\mathrm{H,1}}$ peak vanishes in Fig. 4b (green dashed
line), and the THz-MDCS spectrum is dominated by broad pump–probe (PP) peaks
at $\omega_{t}=\omega_{0}\sim$4 meV, similar to Fig. 2k. This result suggests
that the peaks at $\omega_{\text{H,1}}$ dominating the PP peaks in nonlinear
regime II, as in Fig. 2g (data) and Fig. 4a (theory), provide coherent sensing
of non-perturbative Higgs collective modes underpinning the
finite–momentum–pairing SC phase (Supplementary Note 7) different from the
ground state in perturbative regime I. Next, we turn to the transition from
Higgs to dominant bi-Higgs signals at
$\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}$. We associate this transition
with the development of a time–dependent canted pseudo–spin state, which is
parametrically–driven by amplified phase dynamics at frequency
$\omega_{\mathrm{H,1}}$. The nonlinear enhancement of this phase dynamics in
regime III originates from the development of a phase–amplitude collective
mode at $\omega_{\mathrm{H,1}}$ (Figs. 2o and 2p). This phase–amplitude mode
interacts with a quasi–particle excitation, both with energies
$\sim\omega_{\mathrm{H,1}}$, which amplifies the THz–MDCS sideband peaks at
frequencies $\sim 2\omega_{\mathrm{H,1}}$ (Fig. 2n), at the expense of the
Higgs mode peak at $\omega_{\mathrm{H,1}}$ which dominates in regime II.

We attribute the drastic change in the THz–MDCS spectra with the non-
perturbative emergence of the bi-Higgs frequency sideband in Regime III over
the Higgs peak to parametric excitation of time–dependent pseudo–spin canting
modulation by a phase–amplitude collective mode (Methods sections 1.4 and
1.5). The latter collective mode develops above critical excitation and leads
to pronounced phase dynamics at the Higgs frequency $\omega_{\text{H,1}}$
(Supplementary Figs. 4a and 4b). To further corroborate the transition from
Higgs amplitude to phase–amplitude collective mode, we compare in Figs. 4c and
4d the spectra $E_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ obtained from the
full calculation for 700.0 kV/cm driving with those obtained by turning off
the pseudo–spin canting driven by the $\omega_{\mathrm{H,1}}$ time–periodic
phase fluctuations around the $s_{\pm}$ equilibrium state (Supplementary Figs.
4d and 6, Note 4). Our formulation of the gauge–invariant SC Bloch equations
in terms of two coupled pseudo–spin nonlinear oscillators, presented in the
Methods section 1.4, shows that non–adiabatic pseudo–spin canting is
parametrically driven with time–dependent strength
$\sim|\Delta_{1}|^{2}\sin(2\delta\theta)$. The threshold nonlinear field
dependence of this coupling, shown in Supplementary Fig. 4d, leads to the
strong field dependence of the $\sim 2\omega_{\mathrm{H,1}}$ sideband (Fig.
2p). By comparing Figs. 4c and. 4d, we see that the signals at frequencies
$\omega_{t}=2\omega_{\mathrm{H,1}}-\omega_{0}$ (blue dashed line) and
$\omega_{t}=2\omega_{\mathrm{H,1}}-2\omega_{0}$ (red dashed line) are absent
when the order parameter phase can be approximated by its equilibrium value.
We also compare the full result with a calculation without interband Coulomb
interaction between the electron and hole pockets (Fig. 4e), which diminishes
the bi-Higgs frequency signals. If the inter–band Coulomb coupling exceeds the
intra–band pairing interaction, the Leggett mode phase oscillations lie well
within the quasi–particle continuum (regime I), so they are overdamped (Fig.
2o). Above critical THz driving (regime III), however, the THz–modulated
superfluid density of strongly–Coulomb–coupled electron and hole pockets
(Methods section 1.5 and Supplementary Note 4) enhances the nonlinear coupling
between the phase–amplitude collective mode ($\omega_{\mathrm{H,1}}$) and
quasi–particle excitations ($\omega_{\mathrm{H,1}}$) in the parametrically-
driven SC state that processes supercurrent $J(t)$–oscillations at
$2\omega_{\mathrm{H,1}}$ visualized as the unique bi-Higgs sidebands in the
THz-MDCS spectra (Fig. 2i).

## References

  * [1] Vaswani, C. _et al._ Light-driven Raman coherence as a nonthermal route to ultrafast topology switching in a Dirac semimetal.  _Phys. Rev. X_ 10, 021013 (2020).
  * [2] Fausti, D. _et al._ Light-induced superconductivity in a stripe-ordered cuprate.  _Science_ 331, 189–191 (2011).
  * [3] Knap, M., Babadi, M., Refael, G., Martin, I. & Demler, E.  Dynamical Cooper pairing in nonequilibrium electron-phonon systems.  _Phys. Rev. B_ 94, 214504 (2016).
  * [4] Mitrano, M. _et al._ Possible light-induced superconductivity in K3C${}_{6}0$ at high temperature.  _Nature_ 530, 461–464 (2016).
  * [5] Budden, M. _et al._ Evidence for metastable photo-induced superconductivity in 3C${}_{6}0$.  _Nature Physics_ 17, 611–618 (2021).
  * [6] Isoyama, K. _et al._ Light-induced enhancement of superconductivity in iron-based superconductor FeSe${}_{0}.5$Te${}_{0}.5$.  _Communications Physics_ 4, 160 (2021).
  * [7] Buzzi, M. _et al._ Higgs-mediated optical amplification in a nonequilibrium superconductor.  _Phys. Rev. X_ 11, 011055 (2021).
  * [8] Matsunaga, R. _et al._ Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor.  _Science_ 345, 1145–1149 (2014).
  * [9] Yang, X. _et al._ Lightwave-driven gapless superconductivity and forbidden quantum beats by terahertz symmetry breaking.  _Nat. Photon._ 13, 707–713 (2019).
  * [10] Linder, J. & Robinson, J. W.  Superconducting spintronics.  _Nat. Phys._ 11, 307 (2015).
  * [11] Maag, T. _et al._ Coherent cyclotron motion beyond Kohn’s theorem.  _Nature Physics_ 12, 119–123 (2016).
  * [12] Reimann, J. _et al._ Subcycle observation of lightwave-driven Dirac currents in a topological surface band.  _Nature_ 562, 396–400 (2018).
  * [13] Lingos, P. C., Kapetanakis, M. D., Wang, J. & Perakis, I. E.  Light-wave control of correlated materials using quantum magnetism during time-periodic modulation of coherent transport.  _Communications Physics_ 4, 60 (2021).
  * [14] Luo, L. _et al._ A light-induced phononic symmetry switch and giant dissipationless topological photocurrent in ZrTe5.  _Nat. Mater._ 20, 329–334 (2021).
  * [15] Rajasekaran, S. _et al._ Parametric amplification of a superconducting plasma wave.  _Nat. Phys._ 12, 1012–1016 (2016).
  * [16] Dienst, A. _et al._ Bi-directional ultrafast electric-field gating of interlayer charge transport in a cuprate superconductor.  _Nature Photon._ 5, 485–488 (2011).
  * [17] Nakamura, S. _et al._ Infrared activation of the Higgs mode by supercurrent injection in superconducting NbN.  _Phys. Rev. Lett._ 122, 257001 (2019).
  * [18] Vaswani, C. _et al._ Terahertz second-harmonic generation from lightwave acceleration of symmetry-breaking nonlinear supercurrents.  _Phys. Rev. Lett._ 124, 207003 (2020).
  * [19] Yang, X. _et al._ Terahertz-light quantum tuning of a metastable emergent phase hidden by superconductivity.  _Nat. Mater._ 17, 586 (2018).
  * [20] Cea, T., Castellani, C. & Benfatto, L.  Nonlinear optical effects and third-harmonic generation in superconductors: Cooper pairs versus Higgs mode contribution.  _Phys. Rev. B_ 93, 180507 (2016).
  * [21] Murotani, Y. & Shimano, R.  Nonlinear optical response of collective modes in multiband superconductors assisted by nonmagnetic impurities.  _Phys. Rev. B_ 99, 224510 (2019).
  * [22] Krull, H., Bittner, N., Uhrig, G., Manske, D. & Schnyder, A.  Coupling of Higgs and Leggett modes in non-equilibrium superconductors.  _Nat. Commun._ 7, 11921 (2016).
  * [23] Vaswani, C. _et al._ Light quantum control of persisting Higgs modes in iron-based superconductors.  _Nat. Commun._ 12, 258 (2021).
  * [24] Giorgianni, F. _et al._ Leggett mode controlled by light pulses.  _Nat. Phys._ 15, 341–346 (2019).
  * [25] Udina, M., Cea, T. & Benfatto, L.  Theory of coherent-oscillations generation in terahertz pump-probe spectroscopy: From phonons to electronic collective modes.  _Phys. Rev. B_ 100, 165131 (2019).
  * [26] Chu, H. _et al._ Phase-resolved Higgs response in superconducting cuprates.  _Nature Commun._ 11, 1793 (2020).
  * [27] Schwarz, L., Fauseweh, B. & Tsuji, N. e. a.  Classification and characterization of nonequilibrium Higgs modes in unconventional superconductors.  _Nat. Commun._ 11, 287 (2020).
  * [28] Moor, A., Volkov, A. F. & Efetov, K. B.  Amplitude Higgs mode and admittance in superconductors with a moving condensate.  _Phys. Rev. Lett._ 118, 047001 (2017).
  * [29] Puviani, M., Schwarz, L., Zhang, X.-X., Kaiser, S. & Manske, D.  Current-assisted Raman activation of the Higgs mode in superconductors.  _Phys. Rev. B_ 101, 220507 (2020).
  * [30] Cundiff, S. T. & Mukamel, S.  Optical multidimensional coherent spectroscopy.  _Physics Today_ 66, 44–49 (2013).
  * [31] Kuehn, W., Reimann, K., Woerner, M., Elsaesser, T. & Hey, R.  Two-dimensional terahertz correlation spectra of electronic excitations in semiconductor quantum wells.  _J. Phys. Chem. B_ 115, 5448–5455 (2011).
  * [32] Junginger, F. _et al._ Nonperturbative interband response of a bulk InSb semiconductor driven off resonantly by terahertz electromagnetic few-cycle pulses.  _Phys. Rev. Lett._ 109, 147403 (2012).
  * [33] Lu, J. _et al._ Coherent two-dimensional terahertz magnetic resonance spectroscopy of collective spin waves.  _Phys. Rev. Lett._ 118, 207204 (2017).
  * [34] Johnson, C. L., Knighton, B. E. & Johnson, J. A.  Distinguishing nonlinear terahertz excitation pathways with two-dimensional spectroscopy.  _Phys. Rev. Lett._ 122, 073901 (2019).
  * [35] Mootz, M., Luo, L., Wang, J. & Perakis, l. E.  Visualization and quantum control of light-accelerated condensates by terahertz multi-dimensional coherent spectroscopy.  _Communications Physics_ 5, 47 (2022).
  * [36] Mahmood, F., Chaudhuri, D., Gopalakrishnan, S., Nandkishore, R. & Armitage, N. P.  Observation of a marginal Fermi glass.  _Nature Physics_ 17, 627–631 (2021).
  * [37] Mootz, M., Wang, J. & Perakis, I. E.  Lightwave terahertz quantum manipulation of nonequilibrium superconductor phases and their collective modes.  _Phys. Rev. B_ 102, 054517 (2020).
  * [38] Lee, S. _et al._ Template engineering of Co-doped BaFe2As2 single-crystal thin films.  _Nat. Mater._ 9, 397–402 (2010).
  * [39] Tu, J. J. _et al._ Optical properties of the iron arsenic superconductor ${\text{BaFe}}_{1.85}{\text{Cco}}_{0.15}{\text{As}}_{2}$.  _Phys. Rev. B_ 82, 174509 (2010).
  * [40] Charnukha, A.  Optical conductivity of iron-based superconductors.  _J. Phys.: Condens. Matter._ 26, 253203 (2014).
  * [41] Stephen, M. J.  Transport equations for superconductors.  _Phys. Rev._ 139, A197–A205 (1965).
  * [42] Yang, F. & Wu, M. W.  Gauge-invariant microscopic kinetic theory of superconductivity: Application to the optical response of Nambu-Goldstone and Higgs modes.  _Phys. Rev. B_ 100, 104513 (2019).
  * [43] Anderson, P. W.  Random-phase approximation in the theory of superconductivity.  _Phys. Rev._ 112, 1900–1916 (1958).
  * [44] Nambu, Y.  Quasi-particles and gauge invariance in the theory of superconductivity.  _Phys. Rev._ 117, 648–663 (1960).
  * [45] Liu, C. _et al._ Evidence for a Lifshitz transition in electron-doped iron arsenic superconductors at the onset of superconductivity.  _Nature Physics_ 6, 419–423 (2010).
  * [46] Fernandes, R. M. & Schmalian, J.  Competing order and nature of the pairing state in the iron pnictides.  _Phys. Rev. B_ 82, 014521 (2010).
  * [47] Yang, X. _et al._ Nonequilibrium pair breaking in $\mathrm{Ba}({\mathrm{Fe}}_{1-x}{\mathrm{Co}}_{x}{)}_{2}{\mathrm{As}}_{2}$ superconductors: Evidence for formation of a photoinduced excitonic state.  _Phys. Rev. Lett._ 121, 267001 (2018).

## Acknowledgments

THz spectroscopy work was supported by National Science Foundation 1905981
(M.M.). THz MDCS instrument was supported by the W.M. Keck Foundation (initial
design and commission) and by the U.S. Department of Energy, Office of
Science, National Quantum Information Science Research Centers,
Superconducting Quantum Materials and Systems Center (SQMS) under the contract
No. DE-AC02-07CH11359 (upgrade for improved cryogenic operation). THz spectral
analysis was supported by the Ames Laboratory, the US Department of Energy,
Office of Science, Basic Energy Sciences, Materials Science and Engineering
Division under contract No. DEAC0207CH11358 (L.L.). The work at UW-Madison
(synthesis and characterizations of epitaxial thin films) was supported by the
US Department of Energy (DOE), Office of Science, Office of Basic Energy
Sciences (BES), under award number DE-FG02-06ER46327. Modeling work at the
University of Alabama, Birmingham was supported by the US Department of Energy
under contract # DE-SC0019137 and was made possible in part by a grant for
high performance computing resources and technical support from the Alabama
Supercomputer Authority.