The amplitude coherent oscillations detected by THz two-pulse coherent spectroscopy

The equilibrium complex conductivity spectra, i.e., real and imaginary parts, σ₁(ω) and σ₂(ω), of our epitaxial Ba(Fe_1−xCo_x)₂As₂ (x = 0.08) film²⁷ (Methods) measure the low-frequency quasi-particle (QP) electrodynamics and condensate coherence, respectively (Fig. 1b, c)²⁸. The normal state (black circles) displays Drude-like behavior, while the QP spectral weight in σ₁(ω) is depleted in the SC state due to SC gap openings, seen, e.g., in the 4.2 K trace (red circles). The lowest SC gap value 2Δ₁ ~ 6.8 meV obtained is in agreement with the literature values 6.2–7 meV^29,30 (Methods). Such σ₁(ω) spectral weight depletion is accompanied by an increase of condensate fraction n_s/n₀ (inset, Fig. 1c), extracted from a diverging 1/ω condensate inductive response, marked by blue arrow, e.g., in the 4.2 K lineshape of σ₂(ω) (Fig. 1c). Note that the superfluid density n_s vanishes above T_c ~ 23  K (inset Fig. 1c).

We characterize the THz quantum quench coherent dynamics directly in the time domain^31–33 (Methods) by measuring the responses to two phase-locked THz pulses as differential field transmission of the weak THz probe field ΔE/E₀ (blue circles, Fig. 1d) for THz pump field, E_pump  =  56 kV/cm and as a function of pump–probe time delay Δt_pp. The central pump energy ℏω_pump  =  5.4 meV (gray shade, Fig. 1e) is chosen slightly below the 2Δ₁ gap. Intriguingly, the ΔE/E₀ dynamics reveals a pronounced coherent oscillation, superimposed on the overall amplitude change, which persists much longer than the THz photoexcitation (pink shade). This mode is excited by the quadratic coupling of the pump vector potential, $$ {\mathbf{A}}^2\left(t\right)\propto E_{\mathrm{pump}}^2/{\omega }^2$$ due to the SC equilibrium symmetry⁷. Such coherent responses yield information within the general framework of coherent nonlinear spectroscopy³⁴ (Methods). The origin of the observed coherent ΔE/E₀ oscillation is better illustrated by its Fourier transformation (FT), shown in Fig. 1e. The FT spectrum of the coherent nonlinear signals (blue solid line) displays a pronounced resonance at 6.9 meV, indicative of 2Δ₁ coherent amplitude mode oscillations. This FT spectrum strongly differs from the spectra of both THz pump E_pump(ω) centered at ω_pump ~ 5.4 meV (gray shade) and second harmonic, Anderson pseudo-spin (APS) precession at 2ω_pump from $$ E_{\mathrm{pump}}^2\left(\omega \right)$$ (pink shade). The broadband spectrum of the few-cycle pump pulse used in the experiment overlaps with the mode resonances such that the ΔE/E₀ oscillates with the collective mode frequencies⁹. After the oscillation, time-dependent complex conductivity spectra, σ₁(ω, Δt_pp) and σ₂(ω, Δt_pp), can be measured (Supplementary Figs. 5–6). They show that ΔE/E₀ closely follows the pump-induced change in condensate density, Δn_s/n₀¹³. The THz excitation at E_pump  =  56 kV/cm only reduces n_s slightly, Δn_s/n₀ ~ ΔE/E₀ ~ −3% at Δt_pp  =  5 ps. Furthermore, the pump-induced peak amplitude (blue diamond), marked by the red dashed line in Fig. 1d, diminishes above T_c (inset). These evidence indicate that the measured coherent oscillations reflect the emergence of a hybrid Higgs multi-band collective mode between two Coulomb-coupled lower and higher modes, ω_H,1 and ω_H,2.

Figure 2 reveals a strong temperature dependence of the Higgs mode oscillations. The coherent dynamics of ΔE/E₀ is shown in Fig. 2a for temperatures 4.2–30 K. Approaching T_c from below, the coherent oscillations quickly diminish, as seen by comparing the 4.2 K (black line) and 16 K (gray) traces versus 22 K (cyan) and 24 K (pink) traces. Figure 2b shows the temperature-dependent Fourier spectra of ΔE/E₀, in the range 4–12 meV. Figure 2c plots the integrated spectral weight SW_5→14meV of the amplitude mode (Fig. 2b). The strong temperature dependence correlates the mode with SC coherence. Importantly, while the mode frequency is only slightly red-shifted, by less than 10% before full SW depletion close to T_c, SW is strongly suppressed, by  ~55% at 16 K (T/T_c ~ 0.7). Such a spectral weight reduction in FeSCs is much larger than expected in one-band BCS superconductors or for U = 0 shown later. We also note that the temperature dependence of Higgs oscillations, observed in FeSCs here, has not been measured experimentally in conventional BCS systems, and could represent a key signature of quantum quench dynamics of unconventional SCs. The observed behavior is consistent with our simulations of the hybrid Higgs mode in multi-band SCs with dominant interband U, shown later.

Fig. 2

Temperature dependence of hybrid Higgs mode coherent oscillations in FeSCs.

a  Temporal profiles of  −ΔE/E₀ at various temperatures and for a peak THz pump E-field of E_pump = 56 kV/cm. Traces are offset for clarity. b  Fourier spectra of the Higgs mode oscillations derived from the two-pulse coherent pump–probe signals at different temperatures. Dashed red lines indicate the slight redshift of the Higgs mode frequency with the drastic reduction of mode SW as a function of temperature. c  Temperature dependence of the integrated SW_5→14meV of the Higgs mode in b.

Nonlinear THz field dependence of hybrid Higgs mode in FeSCs

Figure 3 presents distinguishing experimental evidence for the hybrid Higgs mode in FeSCs, which is different from one-band SCs – a highly nonlinear THz electric field dependence of coherent 2Δ₁ oscillations that manifests as a huge SW change, yet with persisting mode frequency, i.e., with only very small redshift. Figure 3a shows the detailed pump-fluence dependence of ΔE/E₀ as a function of time delay, presented as a false-color plot at T = 4.2  K for up to E_pump ~ 600 kV/cm. It is clearly seen that amplitude mode oscillations depend nonlinearly on the E_pump field strength. 1 + ΔE/E₀ (red solid line) at Δt_pp  = 5 ps is shown in Fig. 3b. This, together with the measured  ~1/ω divergence in σ₂(ω, Δt_pp), allows the determination of the condensate fraction n_s(E_pump)/n₀ (blue circles) in the driven state (Supplementary Fig. 6). There shows three different excitation regimes, marked by black dash lines in Fig. 3a–b: (1) in regime #1, the condensate quench is minimal, e.g., n_s/n₀ ≥ 97% below the field E_#1  =  56 kV/cm; (2) Regime #2 displays partial SC quench, where n_s/n₀ is still significant, e.g., condensate fraction  ≈60% at E_#2  =  276 kV/cm; (3) A saturation regime #3 is observed  ~E_#3  = 600 kV/cm, which leads to a slowly changing n_s/n₀ approaching a saturation ~25%. The saturation is expected since below gap THz pump is used, especially ℏω_pump ≪ 2Δ₂ ~ 15–19 meV at the e-like pockets^29,30.

Fig. 3

Nonlinear THz field dependence of hybrid Higgs mode in FeSCs.

a  A 2D false-color plot of ΔE/E₀ as a function of pump E-field strength E_pump and pump–probe delay Δt_pp at 4.2 K. b  THz pump E-field E_pump dependence of 1 + ΔE/E₀ (red line) overlaid with the superfluid density fraction n_s/n₀ (blue circles) after THz pump at Δt_pp = 5 ps and T = 4.2 K. Dashed arrows (black) mark the three pump E-field regimes, i.e., weak, partial, and saturation, identified in the main text. c, d Spectra of coherent Higgs mode oscillations show a distinct non-monotonic dependence as a function of THz pump field, i.e., a rapid increase in the mode amplitude for low pump E-field strengths up to 173 kV/cm, saturation up to 276 kV/cm and significant reduction at higher fields. The blue dashed line marks the resonance of the mode and the redshift of the Higgs mode peak, much smaller than the mode SW change. e  Integrated spectral weight SW_5→14meV of the Higgs mode at various pump E-field strengths, indicating the SW reduction of the Higgs mode from dominant ω_H,1 at low driving fields to ω_H,2 due to the interband interaction and coherent coupling, illustrated in the inset, at high driving fields above E_pump = 276 kV/cm. f, g Contrasting the thermal and THz-driven states of coherent hybrid Higgs mode spectral weight by comparing the mode spectra at the similar superfluid density f n_s/n₀  =  60% and g n_s/n₀  =  25% achieved by THz pump (black solid lines) and temperature (red shades).

Quantum quenching of the single-band BCS pairing interaction has been well established to induce Higgs oscillations with amplitude scaling as 1/$$ \sqrt{{\mathrm{\Delta }}_{\mathrm{SC},\infty }}$$, determined by the long-time asymptotic nonthermal order parameter Δ_SC,∞. The latter decreases with pump field^35–37. Both model and experimental results establish that the Higgs mode amplitude increases with THz pumping until full depletion of the condensate, concurrent with a continuous Higgs resonance redshift to zero^11,35–37. In contrast to this expected behavior for conventional SCs, Fig. 3a and the Fourier spectra of the coherent oscillations, Fig. 3c, d, show a distinct, non-monotonic pump-field dependence of the Higgs mode amplitude that is unique here. Specifically, the Fourier spectra exhibit a clear resonance at low pump fluences, which coincides with the frequency of the lower Higgs mode ω_H,1. This resonance grows quickly up to a field of E_pump  =  173 kV/cm (Fig. 3c), saturates up to 276 kV/cm (Fig. 3d) and then exhibits a significant reduction in pump regime #3, e.g., by more than 50% at 629 kV/cm (black line, Fig. 3d). Therefore, the coherent oscillations in Fig. 3a quickly increase in pump regime #1 and saturate in regime #2, prior to any significant mode resonance redshift (blue dashed line, Fig. 3c). Above this relatively low field regime, the oscillation amplitude starts to decrease at 276 kV/cm, even though there is still more than 60% of condensate, marked in Fig. 3b: the driven state is still far from full SC depletion. This striking SW reduction is also seen in the integrated spectral weight analysis, SW_5→14meV, summarized in Fig. 3e. Most intriguingly, while there is a large reduction of the Higgs mode SW ~ 50% at 629 kV/cm (regime #3), the resonance peak position is nearly persistent, with ≤10% redshift (blue dashed lines, Fig. 3d). These observations of the hybrid Higgs mode differ from any known behavior in one-band SCs, but are consistent with expectations from light-induced nonlinear coherent coupling of Higgs modes in multi-band SCs by a dominant interband interaction U.

The distinct mode amplitude and position variation with pump field extracted from the coherent oscillations clearly show the transition from SW growth and saturation to reduction, marked by the black dashed line at E_pump  =  276 kV/cm (Fig. 3e). The saturation and reduction of SW in the amplitude oscillation, yet with a persisting mode frequency, cannot be explained by any known mechanism. This can arise from the coupling of the two amplitude modes ω_H,1 and ω_H,2 expected in iron pnictides due to the strong inter-pocket interaction U. The coherent coupling process can be controlled and detected nonlinearly by THz two-pulse coherent spectroscopy, as clearly shown in Fig. 3e. We argue that (I) At low driving fields, ω_H,1 dominates the hybrid collective mode due to less damping than ω_H,2 arising from the asymmetry between the electron and hole pockets; (II) For higher fields, SW of ω_H,1 mode decreases due to the coherent transfer to ω_H,2 mode expected for the strong interband interaction in iron pnictides. Moreover, it is critical to note that the THz driving is of highly nonthermal nature, which is distinctly different from that obtained by temperature tuning in Fig. 2. Specifically, Fig. 3f, g compare the hybrid Higgs mode spectra for similar condensate faction n_s/n₀, i.e., ≈60% (f) and 25% (g), induced by tuning either the temperature (red shade) or THz pump (black line). The mode amplitude is clearly much stronger in the THz-driven coherent states than in the temperature tuned ones.

Quantum kinetic calculation of the THz-driven hybrid Higgs dynamics

To put the above hybrid Higgs mode findings on a rigorous footing, we calculate the THz coherent nonlinear spectra (Methods) by extending the gauge-invariant density matrix equations of motion theory of ref. ³⁸, as outlined in the Supplementary Note 6. Using the results of these calculations, we propose a physical mechanism that explains the distinct differences of the Higgs mode resonance in the four-wave mixing spectra between the strong and weak interband interaction limits. For this, we calculate the APS and quantum transport nonlinearities³⁸ driven by two intense phase-locked THz E-field pulses for an effective 3-pocket BCS model of FeSCs³⁹. This model includes both intraband and interband pairing interactions, as well as asymmetry between electron and hole pockets. We thus calculate the nonlinear differential field transmission ΔE/E₀ for two phase-locked THz pulses, which allows for a direct comparison of our theory with the experiment (Supplementary Note 6).

The inset of Fig. 4a presents the calculated ΔE/E₀ (black line), shown together with $$ E_{\mathrm{pump}}^2\left(t\right)$$ of the applied experimental pump pulse (pink shade). The calculated Higgs mode spectra, Fig. 4a (regime I) and Fig. 4b (regime II), are dominated by a resonance close to 6.8 meV for low pump fluences, which corresponds to ω_H,1. This resonance grows up to pump fields E_THz ≈ 320 kV/cm for the parameters used here, with minimal redshift and without any significant SW at ω_H,2 (Fig. 4a). Interestingly, for higher fields (Fig. 4b), we obtain both a redshift and a decrease of the oscillation amplitude. In this regime II, SW emerges close to 15.0 meV, outside of our experimental bandwidth, in the frequency regime of the ω_H,2 Higgs mode. The latter mode is strongly suppressed due to damping induced via electron–hole asymmetry (Supplementary Note 6). Specifically, the ellipticity of the e pockets increases the DOS along the pump field direction and thus increases the damping of mainly the ω_H,2 resonance, which leads to a transfer of oscillator strength to the continuum. This damping has a much smaller influence on the ω_H,1 resonance that arises largely from the hole pockets. Most importantly, the strong interband coupling expected in FeSCs leads to a decrease in the ω_H,1 resonance amplitude, with SW reduction accompanied by a persisting mode frequency. This behavior of the multiband model with strong U is clearly seen in the raw experimental data. Note that, while in regime I we observe an increase in the mode amplitude without any significant redshift, in regime II, the decrease in ω_H,1 resonance is accompanied by a small redshift. This behavior of the hybrid Higgs mode contradicts the one-band behavior, recovered by setting U  =  0, and is in excellent agreement with our experimental observations in the FeSC system (Fig. 3c–e).

Fig. 4

Gauge-invariant quantum kinetic calculation of the THz-driven hybrid Higgs dynamics.

a  Calculated Higgs mode spectra for low pump E-field strengths. Inset: Calculated ΔE/E₀ for peak E_pump  =  380 kV/cm (black) and its comparison with the waveform $$ E_{\mathrm{pump}}^2$$ (pink, shaded) of the applied THz pump pulse in the experiment. b  Calculated Higgs mode spectra for higher field strengths show a decrease in amplitude and redshift of the Higgs mode ω_H,1, consistent with experimentally measured coherent responses in Fig. 3c, d. Note that a second ω_H,2 appears at higher E_pump-field strengths and gets stronger at elevated E_pump-field strengths. Although the ω_H,2 mode is outside the experimental sampling width, it is revealed by the distinct nonlinear THz field dependence of spectral weight controlled by THz pump (Fig. 3d). c E_pump-field dependence of the calculated Higgs mode ω_H,1 spectral weight without interband coupling, i.e., U  =  0 (red circles) and with strong interband coupling U  ≠  0 (blue circles). d  Plot of the Higgs mode frequency 2ω_H,1 as a function of E_pump-field strength without interband coupling (red circles) and for strong interband coupling (blue circles). Our simulations for strong U (blue circles) are in full agreement with the experimental results in Fig. 3c-d and in sharp contrast with the one-band SC results obtained for U  =  0 (red circles). e, f The calculated spectral weight SW_0→14meV e and resonance position f are plotted as a function of temperature for a fixed pump field of 56.0 kV/cm for U  =  0 (red circles) and U  ≠  0 (blue circles). With interband coupling, the SW is strongly suppressed, by about 60% up to a temperature of 0.6T_c, while at the same time the mode frequency is only slightly redshifted, by about 15%, before a full spectral weight depletion is observable towards T_c. These simulations are in agreement with the hybrid Higgs behavior in FeSCs and differ from one-band superconductors showing comparable change of SW and position of the Higgs mode with increasing temperature (red circles).

To scrutinize further the critical role of the strong interband interaction U, we show the fluence dependence of the coherent Higgs SW close to ω_H,1 in Fig. 4c. SW_0→14meV differs markedly between the calculation with strong U  ≠  0 (blue circles) and that without inter-pocket interaction U  =  0 (red circles), which resembles the one-band BCS quench results. Importantly, the Higgs mode SW_0→14meV for strong U grows at low pump fluences (regime I), followed by a saturation and then decrease at elevated E_pump (regime II), consistent with the experiment. Meanwhile, Fig. 4d demonstrates that the resonance frequency remains constant in regime I, despite the strong increase of SW, and then redshifts in regime II, yet by much less than in the one-band system (compare U  ≠  0 (blue circles) vs. U  =  0 (red circles)). Without inter-pocket U (U  =  0, red circles in Fig. 4c–d), the SW of Higgs mode ω_H,1 grows monotonically up to a quench of roughly 90%. A further increase of the pump field leads to a complete quench of the order parameter Δ₁ and a decrease of the SW of Higgs mode ω_H,1 to zero, due to transition from a damped oscillating Higgs phase to an exponential decay. Based on the calculations in Fig. 4c–d, the decrease of the spectral weight with interband coupling appears at a Δ₁ quench close to 15%, while without interband coupling the decrease of spectral weight is only observable close to the complete quench (~90%) of the SC order parameter Δ₁. We conclude from this that the spectral weight decrease at the lower Higgs resonance with low redshift is a direct consequence of the strong coupling between the electron and hole pockets due to large U. This E_pump dependence is the hallmark signature of the Higgs mode in FeSCs and is fully consistent with the Higgs mode behaviors observed experimentally.

Finally, the temperature dependence of the hybrid Higgs mode predicted by our model is shown in Fig. 4e, f for U  =  0 (red circles) and U  ≠  0 (blue circles). With interband coupling, the SW is strongly suppressed, by about 60% up to a temperature of 0.6T_c, while at the same time the mode frequency is only slightly redshifted, by about 15%, before a full spectral weight depletion is observable towards T_c. The strong suppression results from transfer of SW from mode ω_H,1 to the higher mode ω_H,2 with increasing temperature, since the higher SC gap Δ₂ experiences stronger excitation by the applied pump E² with growing T. These simulations are in agreement with the hybrid Higgs behavior in Fig. 2 and differ from one-band superconductors showing comparable change of both SW and position of the Higgs mode with increasing temperature (red circles, Fig. 4e, f). Moreover, our calculation without light-induced changes in the collective effects (only charge-density fluctuations) produces a significantly smaller ΔE/E₀ signal in the non-perturbative excitation regime (Supplementary Fig. 11). Therefore, we conclude that the hybrid Higgs mode dominates over charge-density fluctuations in two-pulse coherent nonlinear signals in FeSCs, due to the different effects of the strong interband U and multi-pocket bandstructure on QPs and on Higgs collective modes.

In summary, we provide distinguishing features for hybrid Higgs modes and thier coherent excitations in multi-band FeSCs, which differ significantly from any previously observed collective mode in other superconducting materials: 2Δ_SC amplitude oscillations displaying a robust mode resonance frequency position despite a large change of its spectral weight, more than 50%, on the THz electric field. This unusual nonlinear quantum behavior provides evidence for Higgs hybridization from the interband quantum entanglement in FeSCs. Such discovery and light-control of the hybrid Higgs mode  in multi-band unconventional SCs inspire future research of quantum tomography of  many-body correlated states in complex systems.

Sample preparation

We measure optimally Co-doped BaFe₂As₂ epitaxial single-crystal thin films²⁷ which are discussed in Supplementary Figs. 1 and 2. They are 60 nm thick, grown on 40 nm thick SrTiO₃ buffered (001)-oriented (La;Sr)(Al;Ta)O₃ (LSAT) single-crystal substrates. The sample exhibits a SC transition at T_c ~ 23  K (Supplementary Fig. 2). The base pressure is below 3 ×  10⁻⁵ Pa and the films were synthesized by pulsed laser deposition with a KrF (248 nm) ultraviolet excimer laser in a vacuum of 3  ×  10⁻⁴ Pa at 730^∘C (growth rate: 2.4 nm/sec). The Co-doped Ba-122 target was prepared by solid-state reaction with a nominal composition of Ba/Fe/Co/As = 1:1.84:0.16:2.2. The chemical composition of the thin film is found to be Ba(Fe_0.92,Co_0.08)₂As_1.8, which is close to the stoichiometry of Ba122 with 8% (atomic %) optimal Co-doping. The PLD targets were made in the same way using the same nominal composition of Ba(Fe_0.92,Co_0.08)₂As_2.2 (Supplementary Note 1). The epitaxial and crystalline quality of the optimally doped BaFe₂As₂ thin films was measured by four-circle X-ray diffraction (XRD) shown in Supplementary Fig. 1 for the out-of-plane θ–2θ scans of the films. As an example, the full width half maximum (FWHM) of the (004) reflection rocking curve of the films is as narrow as 0.7^∘, which indicates high-quality epitaxial thin films. Furthermore, we have also performed chemical, structural, and electrical characterizations of epitaxial Co-doped BaFe₂As₂ (Ba-122) superconducting thin films. We determined the chemical composition of Ba122 thin films by wavelength dispersive spectroscopy (WDS) analyses. The chemical composition of the thin film is found to be Ba(Fe_0.92,Co_0.08)₂As_1.8, which is close to the stoichiometry of Ba122 with 8% (atomic %) optimal Co-doping. We measured temperature-dependent electrical resistivity for superconducting transitions by four-point method (Supplementary Fig. 2). Onset T_c and T_c at zero resistivity are as high as 23.4 and 22.0 K, respectively, and ΔT_c is as narrow as 1.4 K, which are the highest and narrowest values for Ba-122 thin films. In our prior papers, one can also check a zero-field-cooled magnetization T_c and it clearly shows a diamagnetic signal by superconducting quantum interference device (SQUID) magnetometer measurements.

Properties of the thin film Co-doped BaFe₂As₂ samples used

The lowest SC gap values 2Δ₁ ~ 6.8 meV were obtained in our THz conductivity data using a Mattis-Bardeen approach similar to that used in^28–30. Specifically, Supplementary Fig. 7 plots together the THz conductivity of our single-crystal film sample and prior FTIR conductivity of single-crystal samples at the similar doping²⁹. They show an excellent agreement which indicate similar superconducting energy gaps of both samples. The two superconducting energy gaps obtained from the FTIR data are 2Δ₁ ~ 6.2 meV that agrees with our 2Δ₁ value, and 2Δ₂ ~ 14.8–19 meV^29,30. The further analysis of THz conductivity spectra of our sample also allows the determination of other key equilibrium electrodynamics and transport parameters, consistent with the prior literature^3,6. First, we can obtain the plasma frequency ω_p, obtained by fitting the normal state conductivity spectra with the Drude model, which gives the plasma frequency to be   ~194 THz (~804 meV) and the scattering rate  ~2.5 THz (~10 meV). These measurements are in excellent agreement with bulk single-crystal samples at optimal doping²⁹: a normal state Drude plasma frequency of 972 meV, scattering rate 15 meV was obtained there. Second, the London penetration depth λ_L = 378 nm, corresponding to a condensed fraction of  ~50%. These are again consistent with the prior measurements of high quality bulk samples, which show ~50% of the free carriers participating in superfluidity below T_c and a penetration depth of 300 nm.

Two-pulse THz coherent nonlinear spectroscopy

Our THz pump–THz probe setup, illustrated in Supplementary Figs. 3 and 4, can be understood within the general framework of THz 2D spectroscopy. The experiment is driven by a 1 kHz fs laser amplifier⁴⁰ and performed in the collinear geometry with two pulses E_A and E_B with wave vectors $$ {\overrightarrow{k}}_A$$=$$ {\overrightarrow{k}}_B$$=$$ {\overrightarrow{k}}_{\mathrm{NL}}$$. Measuring the electric fields in time domain through electro-optic sampling (EOS) by a third pulse (red shade) allows for phase-resolved detection of the sample response. In general, such a two pulse experiment can provide a number of nonlinear (NL) responses as shown in the Supplementary Table 1. The signals arise from the third order (χ⁽³⁾) nonlinear pump–probe responses of the superconducting state, which are separated from the linear response background. Two main contributions relevant here are measured along the same phase matching direction $$ {\overrightarrow{k}}_{\mathrm{NL}}$$: (1) pump–probe (PP) signals that access condensate quench and recovery; (2) Four-wave mixing (FWM) signals that access amplitude channel coherence and/or density fluctuations. For the sake of consistency, we follow our prior publications¹³ to label the time delays Δt_pp and t_gate corresponding to pulse E_A and E_B in Supplementary Fig. 3, which can be varied independently. For the PP contribution, a polarization response $$ \tilde{P}\left(\omega ,\mathrm{\Delta }{t}_{\mathrm{pp}}\right)$$ is measured and used to obtain the time- and frequency-resolved response functions by performing deconvolution along the t_gate axis. For the FWM contribution, one usually measures it in terms of the inter-pulse delay τ  =  t_gate  −  Δt_pp³⁴ as ΔE_FWM(τ, t_gate)  =  E_AB(t_g  −  Δt_pp, t_gate)  −  E_A(Δt_pp)  −  E_B(t_gate) where E_AB means E_A and E_B are both present. In our experiment, two optical choppers synchronized by the f = 1 kHz laser repetition rate modulate the pump and probe THz beams at f/2 (500 Hz) and f/4 (250 Hz), respectively. The measured data is then divided into 4 channels CH0 to CH3 as shown in the schematic of the timing sequence in Supplementary Fig. 4. After the data acquisition in each channel, we obtain different background-free nonlinear signals as ΔE  =  E_AB  −  E_A  −  E_B  =  E_CH0  +  E_CH1  −  E_CH2  −  E_CH3. We provide more details in supplementary.

Time-resolved complex THz conductivity spectra

Unlike for the Fourier spectra of coherent FWM signals of main interest here, THz PP signals (Supplementary Table 1 and Supplementary Eq. (1), Supplementary Note 2) directly access the time-dependent complex conductivity spectra, i.e., real and imaginary parts, σ₁(ω, Δt_PP), σ₂(ω, Δt_PP), shown in Supplementary Fig. 5 for a fixed time delay Δt_PP = 5 ps. PP signal measurements extend the conventional, equilibrium complex conductivity σ₁(ω) and σ₂(ω) measurement (black and gray circles, Supplementary Fig. 5) to characterize the non-equilibrium post-THz-quench superconducting states. σ₁(ω, Δt_PP), σ₂(ω, Δt_PP) are shown for E-field = 44 kV/cm (red circles) and E-field = 88 kV/cm (blue circles) in Supplementary Fig. 5. The static conductivity can be fitted with Drude-like behavior of quasi-particles in the normal state (black circles) and condensate gaps of 2Δ =  6.8 meV in the SC state (red circles) (discussed later in Supplementary Fig. 7). The PP conductivity spectra identify the minimal quenching of condensate density of few percent for tens of kV/cm E-field used in observing coherent oscillations. Time-dependence of superfluid density can be directly obtained by the diverging σ₂(ω, Δt_PP) ∝ n_s/ω in Supplementary Fig. 5. Furthermore, we also plot Supplementary Fig. 6a from three representative line-cuts in Fig. 3a (main text) and then compare this with the condensate fraction (Supplementary Fig. 6b) extracted from the usual diverging response in σ₂(ω) (Supplementary Fig. 6c). It is clear that, after the oscillation, ΔE/E₀ closely follows the pump-induced change in condensate density, Δn_s/n₀. For example, for weak THz excitation E_pump  =  56 kV/cm and below T_c, the THz pump pulse only reduces n_s slightly, e.g., ΔE/E₀ ∝ Δn_s/n₀ ~ −3%.

Gauge-invariant theory and simulations of THz coherent nonlinear spectroscopy in FeSCs

We start from the microscopic spatial-dependent Boguliobov–de Gennes Hamiltonian for multi-band superconductors

Here the Fermionic field operators $$ {\psi }_{\alpha ,\nu }^{†}\left(\mathbf{x}\right)$$ and ψ_α,ν(x) create and annihilate an electron with spin α in pocket ν; ξ_ν(p + eA(x, t)) corresponds to the dispersion of the pocket with momentum operator p = −i∇_x (ℏ = 1), vector potential A(x, t), and electron charge  −e; μ is the chemical potential while ϕ(x, t) denotes the scalar potential.

The SC complex order parameter components arising from the condensates in the different Fermi sea pockets are

Hamiltonian (1) is gauge-invariant under the gauge transformation

The equation of motion for $$ {\tilde{\rho }}^{\left(\nu \right)}\left(\mathbf{r},\mathbf{R}\right)$$ is now derived by using the Heisenberg equation of motion technique. To simplify the equations of motion, we Fourier transform the obtained exact results with respect to the relative coordinate r and then apply a gradient expansion which is valid for SC systems with condensate center-of-mass spatial fluctuations smoother than the spatial dependence of Cooper pair relative motion. To simplify the problem further, we eliminate the phase of the order parameter $$ {\mathrm{\Delta }}_{{\nu }_0}\left(\mathbf{R}\right)$$ by applying the gauge transformation

In our calculations we solve the gauge-invariant optical Bloch equations (9) for a 3-pocket model with a hole (h) pocket centered at the Γ-point and two electron (e) pockets located at (π, 0) and (0, π). We include the inter e–h pocket interactions (U  =  g_e,h  =  g_h,e) as well as intra-pocket interactions (V_λ  =  g_λ,λ) while inter e–e pocket interactions are neglected for simplicity. The dominance of interband coupling U between e–h pockets over intraband in Fe-based SCs is taken into account by using an interband-to-intraband interaction ratio of r  =  U/V  =  10. The pockets are modeled using the square lattice nearest-neighbor tight-binding dispersion $$ {\xi }_{\nu }\left(\mathbf{k}\right)=-2\left[J_{\nu ,x}\cos \left(k_x\right)+J_{\nu ,y}\cos \left(k_y\right)\right]+{\mu }_{\nu }$$ with hopping parameter J_ν,i and band-offset μ_ν. We choose a circular hole pocket with J_1,x  =  J_1,y  =  10.0 meV and μ₁  = −37.5   meV. We introduce the known particle–hole asymmetry between electron and hole pockets in our system by considering elliptical electron pockets with J_2,x  =  J_3,y  =  −10.0 meV, J_2,y  =  J_3,x  = −50.0 meV, and μ₂  =  μ₃  =  57.5 meV. The latter can lead to coexistence of superconductivity and spin-density wave by changing the doping level. For the doping levels considered here, such asymmetry in our calculation strongly suppresses the ω_H,2 mode in the spectra of coherent ΔE/E dynamics presented in Fig. 4a–b in the main text as discussed in more detail below. We assume s± -pairing symmetry with equilibrium SC order parameters Δ₁  =  3.4 meV for the hole pocket and Δ₂  =  Δ₃  =  9.7 meV for the electron pockets. To directly model our phase-coherent nonlinear pump–probe spectroscopy experiments analogous to 2D THz phase-coherent nonlinear spectroscopy in semiconductors, we calculate the nonlinear differential transmission, ΔE, and not just the order parameter dynamics. The nonlinear differential transmission ΔE/E₀ is obtained by computing the transmitted E-field of both pump and probe pulse, E_pp(t, τ), as a function of gate time t and pump–probe delay τ, as well as the transmitted electric field resulting from the probe pulse, E_probe(t), and the pump pulse, E_pump(t, τ) separately, following the experimental protocols discussed above. Here, the calculated transmitted E-field is given by