* Jacob et al. (2018) Jacob S., Pakmor R., Simpson C. M., Springel V., Pfrommer C., 2018, MNRAS, 475, 570
  * Kazantsev (1968) Kazantsev A. P., 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Kennicutt (1998) Kennicutt Robert C. J., 1998, ApJ, 498, 541
  * Kereš et al. (2005) Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
  * Kolmogorov (1941) Kolmogorov A., 1941, The Proceedings of the USSR Academy of Sciences, 30, 301
  * Kotarba et al. (2010) Kotarba H., Karl S. J., Naab T., Johansson P. H., Dolag K., Lesch H., Stasyszyn F. A., 2010, ApJ, 716, 1438
  * Kotarba et al. (2011) Kotarba H., Lesch H., Dolag K., Naab T., Johansson P. H., Donnert J., Stasyszyn F. A., 2011, MNRAS, 415, 3189
  * Kulsrud (2005) Kulsrud R. M., 2005, Plasma physics for astrophysics
  * Kulsrud & Zweibel (2008) Kulsrud R. M., Zweibel E. G., 2008, Reports on Progress in Physics, 71, 046901
  * Lacki & Beck (2013) Lacki B. C., Beck R., 2013, MNRAS, 430, 3171
  * Libeskind et al. (2020) Libeskind N. I., et al., 2020, MNRAS,
  * Luo et al. (2014) Luo W., Yang X., Zhang Y., 2014, ApJ, 789, L16
  * Lupton et al. (2004) Lupton R., Blanton M. R., Fekete G., Hogg D. W., O’Mullane W., Szalay A., Wherry N., 2004, PASP, 116, 133
  * Man et al. (2016) Man A. W. S., Zirm A. W., Toft S., 2016, ApJ, 830, 89
  * Marinacci & Vogelsberger (2016) Marinacci F., Vogelsberger M., 2016, MNRAS, 456, L69
  * Marinacci et al. (2014) Marinacci F., Pakmor R., Springel V., 2014, MNRAS, 437, 1750
  * Marinacci et al. (2018) Marinacci F., Vogelsberger M., Kannan R., Mocz P., Pakmor R., Springel V., 2018, MNRAS, 476, 2476
  * Marinacci et al. (2019) Marinacci F., Sales L. V., Vogelsberger M., Torrey P., Springel V., 2019, MNRAS, 489, 4233
  * Martin-Alvarez et al. (2018) Martin-Alvarez S., Devriendt J., Slyz A., Teyssier R., 2018, MNRAS, 479, 3343
  * Martin-Alvarez et al. (2020) Martin-Alvarez S., Slyz A., Devriendt J., Gómez-Guijarro C., 2020, MNRAS, 495, 4475
  * Mocz et al. (2014) Mocz P., Vogelsberger M., Sijacki D., Pakmor R., Hernquist L., 2014, MNRAS, 437, 397
  * Monachesi et al. (2019) Monachesi A., et al., 2019, MNRAS, 485, 2589
  * Moreno et al. (2020) Moreno J., et al., 2020, MNRAS, p. 10.1093/mnras/staa2952
  * Moss et al. (2014) Moss D., Sokoloff D., Beck R., Krause M., 2014, A&A, 566, A40
  * Naab & Burkert (2003) Naab T., Burkert A., 2003, ApJ, 597, 893
  * Nulsen & Fabian (2000) Nulsen P. E. J., Fabian A. C., 2000, MNRAS, 311, 346
  * Okamoto et al. (2010) Okamoto T., Frenk C. S., Jenkins A., Theuns T., 2010, MNRAS, 406, 208
  * Pakmor & Springel (2013) Pakmor R., Springel V., 2013, MNRAS, 432, 176
  * Pakmor et al. (2011) Pakmor R., Bauer A., Springel V., 2011, MNRAS, 418, 1392
  * Pakmor et al. (2014) Pakmor R., Marinacci F., Springel V., 2014, ApJ, 783, L20
  * Pakmor et al. (2016a) Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal K., Zhu C., 2016a, MNRAS, 455, 1134
  * Pakmor et al. (2016b) Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016b, ApJ, 824, L30
  * Pakmor et al. (2017) Pakmor R., et al., 2017, MNRAS, 469, 3185
  * Pakmor et al. (2018) Pakmor R., Guillet T., Pfrommer C., Gómez F. A., Grand R. J. J., Marinacci F., Simpson C. M., Springel V., 2018, MNRAS, 481, 4410
  * Pakmor et al. (2020) Pakmor R., et al., 2020, MNRAS,
  * Pfrommer et al. (2017a) Pfrommer C., Pakmor R., Schaal K., Simpson C. M., Springel V., 2017a, MNRAS, 465, 4500
  * Pfrommer et al. (2017b) Pfrommer C., Pakmor R., Simpson C. M., Springel V., 2017b, ApJ, 847, L13
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Powell et al. (1999) Powell K. G., Roe P. L., Linde T. J., Gombosi T. I., De Zeeuw D. L., 1999, Journal of Computational Physics, 154, 284
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
  * Rautiainen & Salo (2000) Rautiainen P., Salo H., 2000, IAU Colloq. 174: Small Galaxy Groups, 209, 330
  * Raymond (1992) Raymond J. C., 1992, ApJ, 384, 502
  * Renaud et al. (2014) Renaud F., Bournaud F., Kraljic K., Duc P. A., 2014, MNRAS, 442, L33
  * Rieder & Teyssier (2017) Rieder M., Teyssier R., 2017, MNRAS, 472, 4368
  * Rodenbeck & Schleicher (2016) Rodenbeck K., Schleicher D. R. G., 2016, A&A, 593, A89
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Rodriguez-Gomez et al. (2017) Rodriguez-Gomez V., et al., 2017, MNRAS, 467, 3083
  * Ruszkowski et al. (2017) Ruszkowski M., Yang H. Y. K., Zweibel E., 2017, ApJ, 834, 208
  * Sales et al. (2012) Sales L. V., Navarro J. F., Theuns T., Schaye J., White S. D. M., Frenk C. S., Crain R. A., Dalla Vecchia C., 2012, MNRAS, 423, 1544
  * Scannapieco et al. (2012) Scannapieco C., et al., 2012, MNRAS, 423, 1726
  * Schleicher et al. (2013) Schleicher D. R. G., Schober J., Federrath C., Bovino S., Schmidt W., 2013, New Journal of Physics, 15, 023017
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Schober et al. (2013) Schober J., Schleicher D. R. G., Klessen R. S., 2013, A&A, 560, A87
  * Scudder et al. (2012) Scudder J. M., Ellison S. L., Torrey P., Patton D. R., Mendel J. T., 2012, MNRAS, 426, 549
  * Sérsic (1963) Sérsic J. L., 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shukurov et al. (2006) Shukurov A., Sokoloff D., Subramanian K., Brand enburg A., 2006, A&A, 448, L33
  * Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  * Sijacki et al. (2012) Sijacki D., Vogelsberger M., Kereš D., Springel V., Hernquist L., 2012, MNRAS, 424, 2999
  * Sparre & Springel (2016) Sparre M., Springel V., 2016, MNRAS, 462, 2418
  * Sparre & Springel (2017) Sparre M., Springel V., 2017, MNRAS, 470, 3946
  * Sparre et al. (2020) Sparre M., Pfrommer C., Ehlert K., 2020, MNRAS, 499, 4261
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel (2015) Springel V., 2015, N-GenIC: Cosmological structure initial conditions (ascl:1502.003)
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Steinwandel et al. (2019) Steinwandel U. P., Beck M. C., Arth A., Dolag K., Moster B. P., Nielaba P., 2019, MNRAS, 483, 1008
  * Su et al. (2018) Su K.-Y., Hayward C. C., Hopkins P. F., Quataert E., Faucher-Giguère C.-A., Kereš D., 2018, MNRAS, 473, L111
  * Subramanian (1998) Subramanian K., 1998, MNRAS, 294, 718
  * Tacchella et al. (2019) Tacchella S., et al., 2019, MNRAS, 487, 5416
  * Teyssier et al. (2010) Teyssier R., Chapon D., Bournaud F., 2010, ApJ, 720, L149
  * Thomas et al. (2020) Thomas T., Pfrommer C., Enßlin T., 2020, ApJ, 890, L18
  * Thorp et al. (2019) Thorp M. D., Ellison S. L., Simard L., Sánchez S. F., Antonio B., 2019, MNRAS, 482, L55
  * Toomre & Toomre (1972) Toomre A., Toomre J., 1972, ApJ, 178, 623
  * Torrey et al. (2012) Torrey P., Cox T. J., Kewley L., Hernquist L., 2012, ApJ, 746, 108
  * Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
  * Vogelsberger et al. (2014a) Vogelsberger M., et al., 2014a, MNRAS, 444, 1518
  * Vogelsberger et al. (2014b) Vogelsberger M., et al., 2014b, Nature, 509, 177
  * Weinberger et al. (2019) Weinberger R., Springel V., Pakmor R., 2019, ApJS,
  * Welker et al. (2017) Welker C., Dubois Y., Devriendt J., Pichon C., Kaviraj S., Peirani S., 2017, MNRAS, 465, 1241
  * Wetzel et al. (2009) Wetzel A. R., Cohn J. D., White M., 2009, MNRAS, 395, 1376
  * Zweibel (2017) Zweibel E. G., 2017, Physics of Plasmas, 24, 055402
  * van de Voort et al. (2021) van de Voort F., Bieri R., Pakmor R., Gómez F. A., Grand R. J. J., Marinacci F., 2021, MNRAS, 501, 4888

##  Appendix A Impact of MHD on the ISM

Figure 16: Top panel: star formation rate surface density as a function of gas
surface density for 1349-3M, as seen at a lookback time of $\sim$4 Gyr. The
dashed line shows a Kennicutt-Schmidt relation (Schmidt, 1959; Kennicutt,
1998) with exponent 1.5. The dotted line indicates the approximate position of
the cut-off in the star formation rate. Bottom panel: as above, but for
1349-3H. Both follow the same relation, despite the strong amplification of
the magnetic field in 1349-3M during this time.

In Section 2.2, we claimed that it was not necessary to recalibrate our ISM
subgrid model when including magnetic fields in the simulation. This statement
is non-trivial: whilst the effective pressure in our ISM subgrid model is a
function of density only in the hydrodynamic simulations (Springel &
Hernquist, 2003), in the MHD simulations there is an additional magnetic
pressure term to consider. We can check the influence of this additional term
on the ISM model by investigating its impact on the relation between the star
formation rate surface density, $\Sigma_{\text{SFR}}$, and the gas surface
density, $\Sigma_{\text{gas}}$. Stars form probabilistically out of our
simulated ISM with a gas consumption time-scale set to match that observed by
Kennicutt (1998) for disc galaxies in the local Universe. This results in the
ISM following the well-known Kennicutt-Schmidt (Schmidt, 1959; Kennicutt,
1998) relation of: $\Sigma_{\text{SFR}}\propto(\Sigma_{\text{gas}})^{n}$. If
our ISM subgrid model required recalibration, the form of this relation would
be dependent on the physics model used.

We show the relation between $\Sigma_{\text{SFR}}$ and $\Sigma_{\text{gas}}$
for 1349-3M and 1349-3H in Fig. 16. The 1349 simulations were chosen due to
the particularly strong amplification in 1349-3M, but naturally the results
apply to all our simulations. We have chosen a time when the magnetic field is
highly amplified in the MHD simulation. At later times we observe a similar
relation, but with the gas density covering a much narrower range. In both
cases, the surface densities are calculated by taking a face-on projection
with depth $\pm$1 kpc from the midplane. We choose this depth to make sure
that we are predominantly considering gas in the disc. It can be seen that for
both physics models, the star formation rate surface density follows the
Kennicutt-Schmidt relation with exponent $n=1.5$ (cf. Springel & Hernquist,
2003). This is true over a broad range of values. At lower gas surface density
values, the relation becomes more scattered as the star formation threshold
density is approached. The peak of the distribution is higher in the MHD
simulation than the hydrodynamic analogue, which is a result of morphological
differences between the two remnants. The strong similarity between both
relations otherwise supports our assertion that the ISM subgrid model must not
be recalibrated when introducing magnetic fields.

##  Appendix B Galaxy Tracking and Black Hole Dynamics

Figure 17: The distance between the galactic centre and the closest black hole
for all high-resolution simulations as a function of time. In general, this
distance stays well below 5 kpc, confirming the reliability of our galaxy
tracking method. The ‘wandering’ nature of the black hole in simulation
1526-3H likely contributes to the unusual morphology displayed by the
corresponding merger remnant.

In Section 2.3, we claimed that tracking a galaxy between snapshots is
frequently akin to tracking the black hole particle that resides in that
galaxy. We support this claim in Fig. 17, where we show the distance between
the most bound gas cell in the galaxy (i.e. the galactic centre) and the
closest black hole for each high-resolution simulation. For the vast majority
of snapshots, this distance is always well under 5 kpc, confirming the
validity of our tracking method. An exception to this rule is simulation
1526-3H, where the black hole is not well-tied to the galactic centre. In this
merger scenario, the secondary progenitor passes directly through the primary
progenitor. For a short time immediately afterwards, the black hole then
‘hitches a ride’, becoming gravitationally bound to the merging galaxy. This
can be confirmed by comparing its distance from the main galaxy between 7.11
and 6.35 Gyr to the distance between the merging galaxies, as seen in Fig. 1.

On the black hole’s return to the main galaxy, it never quite loses its newly-
gained orbital energy, oscillating around the galactic centre instead. During
this time, the black hole continues to accrete gas, and consequently continues
to inject energy into neighbouring gas cells. The subsequent AGN outbursts
disrupt the gas, producing a similar morphology to that of 1605-3M (see Fig.
6). This similarity is unexpected as the black hole in 1526-3H grows only a
quarter as large as that in 1605-3H post-merger, meaning that the AGN
outbursts should be significantly less influential (see Section 2.2). Indeed,
the black hole growth in 1526-3H is similar to that of 1349-3H and 1330-3H,
both of which show no signatures of strong outbursts in their gas morphology.
We therefore argue that it is the unlocalised nature of the feedback in
1526-3H that is behind the strongly disrupted morphology. This, in turn,
produces the unusual stellar morphology seen in this simulation.

##  Appendix C Evidence for a Small-scale Dynamo in Action

In Section 3.3.3, we discussed the development of a small-scale dynamo in our
simulations as a function of increasing resolution. In Fig. 18, we support
this with kinetic and magnetic energy power spectra for each of our high-
resolution MHD simulations. These were created in the same manner as for Fig.
15 over the same radius of 5 kpc within a zero-padded box of 10 kpc across.
The profiles are shown at 0.5 Gyr after the respective time of periapsis in
each simulation (see Section 2.4). In choosing this time, we show the power
spectra at a similar point in the evolution of the magnetic field. It can be
seen that for each simulation, the kinetic energy exhibits a Kolmogorov-like
spectrum, which is consistent with the volume-filling phase of the gas being
both turbulent and subsonic. The difference in normalisation in each plot can
be mostly explained by the difference in mass evolution. At high $k$ values,
the magnetic energy saturates at around 50 per cent of the kinetic energy,
which is consistent with subsonic turbulent box simulations where the forcing
was solenoidal (Federrath, 2016). In each simulation, the magnetic field has
grown similarly quickly. Indeed, the magnetic energy is almost fully saturated
in each case, with the peak magnetic energy occuring at the driving scale of
the turbulence. This means almost the entire magnetic energy spectrum is in
the non-linear dynamo phase, and consequently there is little evidence of the
Kazantsev-like slope at large scales. The decline of the power spectra at
scales larger than 10 kpc is an artefact of the zero-padding we use. Overall,
Fig. 18 supports our understanding of the growth of the magnetic field, as
analysed in Section 3.3.3.

Figure 18: Kinetic and magnetic energy power spectra for the highest-
resolution MHD simulations, calculated for gas within a sphere of 5 kpc
centred on the galactic centre. The power spectra are shown at 0.5 Gyr after
the time of periapsis for each simulation. The black dotted lines show the
slopes of a Kolmogorov spectrum ($\propto k^{-5/3}$) (Kolmogorov, 1941), which
is theoretically expected for a small-scale dynamo resulting from
incompressible turbulence. In each case, the magnetic field is strongly
saturated.


# Light Quantum Control of Persisting Higgs Modes in Iron-Based
Superconductors

C. Vaswani1∗, J. H. Kang2∗, M. Mootz3∗, L. Luo1, X. Yang1, C. Sundahl2, D.
Cheng1, C. Huang1, R. H. J. Kim1, Z. Liu1, Y. G. Collantes4, E. E. Hellstrom4,
I. E. Perakis3, C. B. Eom2 and J. Wang1† 1Department of Physics and Astronomy,
Iowa State University, and Ames Laboratory, Ames, IA 50011 USA.
2Department of Materials Science and Engineering, University of Wisconsin-
Madison, Madison, WI 53706, USA.
3Department of Physics, University of Alabama at Birmingham, Birmingham, AL
35294-1170, USA.
4Applied Superconductivity Center, National High Magnetic Field Laboratory,
Florida State University, Tallahassee, FL 32310, USA.

###### Abstract

The Higgs mechanism, i.e., spontaneous symmetry breaking of the quantum
vacuum, is a cross-disciplinary principle, universal for understanding dark
energy, antimatter and quantum materials, from superconductivity to magnetism.
Yet, Higgs modes in one-band superconductors (SCs) are currently under debate
due to their competition with charge-density fluctuations. A distinct Higgs
mode, controllable by terahertz (THz) laser pulses, can arise in multi-band,
unconventional SCs via strong interband Coulomb interaction, but is yet to be
accessed. Here we both discover and demonstrate quantum control of such
collective mode in iron-based high-temperature superconductors. Using two-
pulse, phase coherent THz spectroscopy, we observe a tunable and coherent
2$\Delta_{\mathrm{SC}}$ amplitude oscillation of the complex order parameter
in such SC with coupled lower and upper bands. The nonlinear dependence of the
amplitude mode oscillations on the THz driving fields is distinct from any
one-band and conventional SC results: we observe a large nonlinear change of
resonance strength, yet with a persisting mode frequency. We argue that this
result provides compelling evidence for a transient coupling between the
electron and hole amplitude modes via strong interband coherent interaction.
To support this scenario, we perform quantum kinetic modeling of a hybrid
Higgs mechanism without invoking extra disorder or phonons. In addition to
distinguishing between collective modes and charge fluctuations, the light
quantum control of multiband SCs can be extended to probe and manipulate many-
body entanglement and hidden symmetries in different quantum materials.

Phase coherence between multiple SC condensates in different strongly
interacting bands is well–established in iron-based superconductors (FeSCs).
As illustrated in Fig. 1a, a dominant Coulomb coupling between the $h$\- and
$e$-like Fermi sea pockets, unlike in other SCs, is manifested by, e.g.,
$s\pm$ pairing symmetry Chubukov:2012 ; Chubukov:2015 , spin density wave
resonant peaks and nesting wave vectors (black arrows) ref5 ; 1 ; 2 ; ref4 .
Experimental evidence for Higgs amplitude coherent excitations in FeSCs has
not been reported yet, despite recent progress in non-equilibrium
superconductivity and collective modes giorgianni2019leggett ; Matsunaga:2013
; matsunaga2014 ; Cea2016 ; Rajasekaran ; n1 ; n2 ; n3 ; n4 ; Wu2019 ;
Udina:2019 ; Kumar:2019 ; Manske:2020 ; Kaiser2020 ; Podolsky2020 .

The condensates in different bands of Ba(Fe1-xCox)2As2 studied here, shown in
Fig. 1a, are coupled by the strong interband $e$–$h$ interaction $U$ (blue
double arrow vector), which is about one order of magnitude stronger than the
intraband interaction $V$ (gray and red thin lines) Fernandes:2010 ;
Akbari:2013 . For $U\gg V$, the formation channels of collective modes are
distinct from one-band SC matsunaga2014 ; Cea2016 ; Cea2018 ; Udina:2019  and
multiband MgB2 with dominant intraband interaction, $V\gg U$ Shimano2019 ;
Aoki2017 ; krull2016 . For the latter, only Leggett modes are observed thus
far giorgianni2019leggett ; ref1 . In contrast, for $U\gg V$, one expects
Higgs amplitude modes arising from the condensates in all Coulomb–coupled
bands, i.e., in the $h$ pocket at the $\Gamma$-point (gray circle, mode
frequency $\omega_{\mathrm{H,1}}$), and in the two $e$ pockets at (0, $\pi$)
and ($\pi$, 0) (red ellipses, $\omega_{\mathrm{H,2}}$). A single-cycle THz
oscillating field (red pulse) can act like a “quantum quench”, with impulsive
non-adiabatic driving of the Mexican-hat-like quantum fields (dark green) and,
yet, with minimum heating of other degrees of freedoms. Consequently, the
multi-band condensates are forced out of the free energy minima, since they
cannot follow the quench adiabatically. Most intriguingly, such coherent
nonlinear driving not only excites amplitude mode oscillations in the
different Fermi sea pockets, but also transiently modifies their coupling,
assisted by the strong interband interaction $U$. Such coherent transient
coupling can be regarded as nonlinear amplitude mode hybridization with a
time-dependent phase coherence. In this way, THz laser fields can manipulate
phase-coherent, hybrid Higgs emerging collective modes.

Figure 1: The 2$\Delta_{\mathrm{SC}}$ oscillations detected by two-pulse THz
coherent spectroscopy of multi-band FeSCs. $\textbf{a},$ Illustration of
coherent excitation of hybrid Higgs mode via THz quantum quench. An effective
three-band model has a $h$ pocket at the $\Gamma$ point and two $e$ pockets at
X/Y points, with strong inter- (blue) and weak intra-band (gray) interactions
marked by arrows. b,c, Real and imaginary parts of the complex THz
conductivity spectra $\sigma_{1}(\omega)$ and $\sigma_{2}(\omega)$ in the
superconducting (4.2 K) and normal states (23 K) of Ba(Fe1-xCox)2As2 (x=0.08)
in equilibrium. Inset of c shows the temperature dependence of the superfluid
density normalized to its value at 4.2K, n${}_{s}(T)$/n0, as determined from
$1/\omega$ divergence of $\sigma_{2}(\omega)$ (blue arrow in c). $\textbf{d},$
Differential THz transmission $\Delta E/E_{0}$ (blue circles) measured by
phase-locked two-THz-pulse pump–probe spectroscopy at 4.2 K shows pronounced
coherent oscillations for a peak THz field strength of Epump=56 kV/cm. The
pink shaded curve denotes the square of the pump THz waveform
E${}_{\mathrm{pump}}^{2}$. Inset: Temperature dependence of the peak amplitude
of $\Delta E/E_{0}$. $\textbf{e},$ Fourier spectrum of the coherent
oscillations in $\Delta E_{\mathrm{}}/E_{0}$ (blue line) exhibits a resonance
peak at $\sim$6.9 meV (blue line) and is distinct from both the pump Epump
(gray) and pump-squared E${}_{\mathrm{pump}}^{2}$ (pink) spectra.

Here we present evidence of Higgs modes that are controlled by THz-field-
driven interband quantum entanglement in a multi-band SC, optimally-doped
Ba(Fe1-xCox)2As2, using two phase-locked near-single–cycle THz laser fields.
We thus reveal a striking nonlinear THz field dependence of coherent amplitude
mode oscillations: quick increase to maximum spectral weight (SW) with
negligible mode frequency shift, followed by a huge SW reduction by more than
50$\%$, yet with robust mode frequency position, with less than 10$\%$
redshift. These distinguishing features of the observed collective mode are
different from any one-band and conventional SC results and predictions so
far. Instead, they are predicted by our quantum kinetic calculation, which
identifies the key role of the interband interaction $U$ for coherently
coupling two amplitude modes, in the $h$\- and $e$-like Fermi sea pockets, and
for controlling the SW of the lower Higgs mode observed in the experiment.

The equilibrium complex conductivity spectra, i.e., real and imaginary parts,
${\sigma}_{1}(\omega)$ and ${\sigma}_{2}(\omega)$, of our epitaxial
Ba(Fe1-xCox)2As2 (x=0.08) film lee2010  (Methods) measure the low-frequency
quasi-particle (QP) electrodynamics and condensate coherence, respectively
(Figs. 1b and 1c) Xu . The normal state (black circles) displays Drude–like
behavior, while the QP spectral weight in ${\sigma}_{1}(\omega)$ 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$\Delta_{1}\sim$6.8 meV obtained is in
agreement with the literature values 6.2-7 meV ref2 ; ref3  (Methods). Such
${\sigma}_{1}(\omega)$ spectral weight depletion is accompanied by an increase
of condensate fraction $n_{\mathrm{s}}/n_{\mathrm{0}}$ (inset, Fig. 1c),
extracted from a diverging $1/\omega$ condensate inductive response, marked by
blue arrow, e. g., in the 4.2 K lineshape of ${\sigma}_{2}(\omega)$ (Fig. 1c).
Note that superfluid density $n_{\mathrm{s}}$ vanishes above T${}_{c}\sim$23 K
(inset Fig. 1c).

We characterize the THz quantum quench coherent dynamics directly in the time
domain (Methods) by measuring the pump–probe responses to two phase-locked THz
pulses as differential field transmission $\Delta E_{\mathrm{}}/E_{0}$ (blue
circles, Fig. 1d) for THz pump field, $E_{\mathrm{pump}}=$56 kV/cm as a
function of pump–probe time delay $\Delta t_{\mathrm{pp}}$. The central pump
energy $\hbar\omega_{\mathrm{pump}}=$5.4 meV (gray shade, Fig. 1e) is chosen
slightly below the 2$\Delta_{1}$ gap. Intriguingly, the $\Delta E/E_{0}$
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}(t)\propto
E^{2}_{\mathrm{pump}}/\omega^{2}$ due to the SC equilibrium symmetry
matsunaga2014 . Such coherent responses yield information within the general
framework of 2D coherent nonlinear spectroscopy rupert  (Methods). The origin
of the observed coherent $\Delta E/E_{0}$ 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.9meV, indicative of 2$\Delta_{1}$ coherent amplitude mode oscillations.
This FT spectrum strongly differs from the spectra of both THz pump
E${}_{\mathrm{pump}}(\omega)$ centered at $\omega_{{\mathrm{pump}}}\sim$5.4
meV (gray shade) and second harmonic, Anderson pseudo-spin (APS) precession at
2$\omega_{\mathrm{pump}}$ from E${}^{2}_{\mathrm{pump}}(\omega)$ (pink shade).
This is a consequence of the broadband spectrum of the few-cycle pump pulse
used in the experiment which overlaps with the mode resonances such that
$\Delta E/E_{0}$ oscillates with the collective mode frequencies Udina:2019 .
After the oscillation, time-dependent complex conductivity spectra,
$\sigma_{1}(\omega,\Delta t_{\mathrm{pp}})$ and $\sigma_{2}(\omega,\Delta
t_{\mathrm{pp}})$, can be measured (Figs.S5-S6, Supplementary). They show that
$\Delta E/E_{0}$ closely follows the pump-induced change in condensate
density, $\Delta n_{\mathrm{s}}/n_{\mathrm{0}}$ n1 . The THz excitation at
$E_{\mathrm{pump}}=$56 kV/cm only reduces $n_{\mathrm{s}}$ slightly, $\Delta
n_{\mathrm{s}}/n_{\mathrm{0}}\sim\Delta E_{\mathrm{}}/E_{0}\sim-3~{}\%$ at
$\Delta t_{\mathrm{pp}}=5$ ps. Furthermore, the pump-induced peak amplitude
(blue diamond), marked by the red dashed line in Fig. 1d, diminishes above Tc
(inset). The measured coherent oscillations reflect the emergence of a hybrid
Higgs multi-band collective mode between two Coulomb–coupled lower and higher
modes, $\omega_{\mathrm{H,1}}$ and $\omega_{\mathrm{H,2}}$, as shown later.

Fig. 1. The 2$\Delta_{\mathrm{SC}}$ oscillations detected by two-pulse THz
coherent spectroscopy of multi-band FeSCs. $\textbf{a},$ Illustration of
coherent excitation of hybrid Higgs mode via THz quantum quench. An effective
three-band model has a $h$ pocket at the $\Gamma$ point and two $e$ pockets at
X/Y points, with strong inter- (blue) and weak intra-band (gray) interactions
marked by arrows. b,c, Real and imaginary parts of the complex THz
conductivity spectra $\sigma_{1}(\omega)$ and $\sigma_{2}(\omega)$ in the
superconducting (4.2 K) and normal states (23 K) of Ba(Fe1-xCox)2As2 (x=0.08)
in equilibrium. Inset of c shows the temperature dependence of the superfluid
density normalized to its value at 4.2K, n${}_{s}(T)$/n0, as determined from
$1/\omega$ divergence of $\sigma_{2}(\omega)$ (blue arrow in c). $\textbf{d},$
Differential THz transmission $\Delta E/E_{0}$ (blue circles) measured by
phase-locked two-THz-pulse pump–probe spectroscopy at 4.2 K shows pronounced
coherent oscillations for a peak THz field strength of Epump=56 kV/cm. The
pink shaded curve denotes the square of the pump THz waveform
E${}_{\mathrm{pump}}^{2}$. Inset: Temperature dependence of the peak amplitude
of $\Delta E/E_{0}$. $\textbf{e},$ Fourier spectrum of the coherent
oscillations in $\Delta E_{\mathrm{}}/E_{0}$ (blue line) exhibits a resonance
peak at $\sim$6.9 meV (blue line) and is distinct from both the pump Epump
(gray) and pump-squared E${}_{\mathrm{pump}}^{2}$ (pink) spectra.

Figure 2: Temperature dependence of hybrid Higgs coherent oscillations in
FeSCs. $\textbf{a},$ Temporal profiles of $-\Delta E_{\mathrm{}}/E_{0}$ at
various temperatures and for a peak THz pump E-field of Epump=56 kV/cm. Traces
are offset for clarity. $\textbf{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. $\textbf{c},$ Temperature dependence of the integrated
SW${}_{5\rightarrow 14~{}\mathrm{meV}}$ of the Higgs mode in b.

Figure 2 reveals a strong temperature dependence of the Higgs oscillations.
The coherent dynamics of $\Delta E/E_{0}$ is shown in Fig. 2a for temperatures
4.2 K–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. Fig. 2b shows the
temperature–dependent Fourier spectra of $\Delta E/E_{0}$, in the range 4–12
meV. Fig. 2c plots the integrated spectral weight SW${}_{5\rightarrow
14~{}\mathrm{meV}}$ 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 $\sim$55 $\%$ at 16
K ($T/T_{c}\sim$0.7). Such a spectral weight reduction in FeSCs is much larger
than that expected in one-band BCS superconductors or for $U$=0 shown later
(Figs. 4e-4f). 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. Note that the observed behavior is
consistent with our simulations of the Higgs mode in multi-band SCs with
dominant interband $U$, shown later in Figs. 4e-4f.

Figure 3 presents the distinguishing experimental evidence for the
unconventional Higgs mode in FeSCs, which is different from one-band SCs – a
highly nonlinear THz electric field dependence of coherent $2\Delta_{1}$
oscillations that manifests as a huge SW change, yet with persisting mode
frequency, with only very small redshift. Fig. 3a shows the detailed pump-
fluence dependence of $\Delta E_{\mathrm{}}/E_{0}$ as a function of time
delay, presented as a false-color plot at $T$=4.2 K for up to
E${}_{\mathrm{pump}}\sim 600$ kV/cm. It is clearly seen that amplitude mode
oscillations depend nonlinearly on the $E_{\mathrm{pump}}$ field strength.
1$+\Delta E_{\mathrm{}}/E_{0}$ (red solid line) at $\Delta t_{\mathrm{pp}}=5$
ps is shown in Fig. 3b. This, together with the measured $\sim 1/\omega$
divergence in $\sigma_{2}(\omega,\Delta t_{\mathrm{pp}})$, allows the
determination of the condensate fraction
$n_{\mathrm{s}}(E_{\mathrm{pump}})/n_{\mathrm{0}}$ (blue circles) in the
driven state (Fig. S6, Supplementary). This shows three different excitation
regimes, marked by black dash lines in Figs. 3a-3b: (1) in regime #1, the
condensate quench is minimal, e. g., $n_{\mathrm{s}}/n_{\mathrm{0}}\geq$97%
below the field $E_{\\#1}=56$ kV/cm; (2) Regime #2 displays partial SC quench,
where $n_{\mathrm{s}}/n_{\mathrm{0}}$ is still significant, e. g., condensate
fraction $\approx$60% at $E_{\\#2}=276$ kV/cm; (3) A saturation regime #3 is
observed $\sim E_{\\#3}=600$ kV/cm, which leads to a slowly changing
$n_{\mathrm{s}}/n_{\mathrm{0}}$ approaching 25 %. The saturation is expected
since the below gap THz pump is used, especially
$\hbar\omega_{pump}\ll$2$\Delta_{2}\sim$15-19 meV at the $e$ pockets ref2 ;
ref3 .

Figure 3: Nonlinear THz field dependence of spectral weight of Higgs mode in
FeSCs. $\textbf{a},$ A 2D false-color plot of $\Delta E_{\mathrm{}}/E_{0}$ as
a function of pump E-field strength Epump and pump–probe delay $\Delta
t_{\mathrm{pp}}$ at 4.2 K. $\textbf{b},$ THz pump $E$-field Epump dependence
of $1+\Delta E_{\mathrm{}}/E_{0}$ (red line) overlaid with the superfluid
density fraction ns/n0 (blue circles) after THz pump at $\Delta
t_{\mathrm{pp}}$=5ps and T=4.2 K. Dashed arrows 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 173kV/cm, saturation up to
276kV/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. $\textbf{e},$ Integrated spectral weight
SW${}_{5\rightarrow 14~{}\mathrm{meV}}$ of the Higgs mode at various pump
$E$-field strengths, indicating the SW reduction of the Higgs mode from
dominant $\omega_{\mathrm{\mathrm{H},1}}$ at low driving fields to
$\omega_{\mathrm{\mathrm{H},2}}$ due to the interband interaction and coherent
coupling at high driving fields above $E_{\mathrm{pump}}=$276 kV/cm (inset).
f, g Contrasting the thermal and THz-driven states of coherent hybrid Higgs
responses by comparing the mode spectra at the similar superfluid density f
$n_{s}/n_{0}=60\%$ and g $n_{s}/n_{0}=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{\Delta_{SC,\infty}}$, determined by the long–time asymptotic
nonthermal order parameter $\Delta_{SC,\infty}$. The latter decreases with
pump field Yuzbashyan:2006 ; Axt:2007 ; Forster:2017 . 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 Yuzbashyan:2006 ; Axt:2007 ;
Forster:2017 ; Matsunaga:2013 . In contrast to this expected behavior for
conventional SCs, Fig. 3a and the Fourier spectra of the coherent
oscillations, Figs. 3c and 3d, show a 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 $\omega_{\mathrm{\mathrm{H},1}}$. This
resonance grows quickly up to a field of $E_{\mathrm{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). This non-monotonic field dependence clearly shows that 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\rightarrow 14~{}\mathrm{meV}}$, summarized
in Fig. 3e. Most intriguingly, while there is a large reduction of the Higgs
mode SW$\sim$50% at 629 kV/cm (regime #3), the resonance peak position is
nearly persistent, with $\leq$10 $\%$ red shift (blue dashed lines, Fig. 3d).
These observations differ from any known behavior in one-band SCs, but are
consistent with expectations from coherent coupling of Higgs modes in multi-
band SCs by a dominant interband interaction $U$, discussed below.

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_{\mathrm{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 $\omega_{\mathrm{\mathrm{H},1}}$ and
$\omega_{\mathrm{\mathrm{H},2}}$ expected in iron pnictides due to the strong
inter-pocket interaction $U$. As we demonstrate theoretically below, the
coherent coupling process shown in Fig. 3e can be controlled and detected
nonlinearly by two phase-locked THz pulses. We argue that (I) At low driving
fields, $\omega_{\mathrm{\mathrm{H},1}}$ dominates the hybrid collective mode
due to less damping than $\omega_{\mathrm{\mathrm{H},2}}$ arising from the
asymmetry between the electron and hole pockets; (II) For higher fields, SW of
$\omega_{\mathrm{\mathrm{H},1}}$ mode decreases due to the strong coupling
with $\omega_{\mathrm{\mathrm{H},2}}$ expected for the strong inter–band
interaction in iron pnictides. Moreover, it is critical to note that the THz
driving is of highly non-thermal nature, which is distinctly different from
that obtained by temperature tuning in Fig. 2. Specifically, Figs. 3f and 3g
compare the hybrid Higgs mode spectra for similar condensate faction
$n_{\mathrm{s}}/n_{\mathrm{0}}$, i.e., $\approx$60% (f) and 25% (g), induced
by tuning either the temperature (red shade) or THz pump (black line). The
mode is clearly much stronger in the THz driven coherent states than in the
temperature tuned ones.

Figure 4: Fig. 4. Gauge-invariant quantum kinetic calculation of the THz-
driven hybrid Higgs dynamics. $\textbf{a},$ Calculated Higgs mode spectra for
low pump $E$-field strengths. Inset: Calculated $\Delta E_{\mathrm{}}/E_{0}$
for peak Epump = 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. $\textbf{b},$ Calculated Higgs mode spectra for higher field
strengths show a decrease in amplitude and redshift of the Higgs mode
$\omega_{H,1}$, consistent with experimentally measured coherent responses in
Figs. 3c and 3d. Note that a second $\omega_{H,2}$ appears at higher
$E_{\mathrm{pump}}$-field strengths and gets stronger at elevated
$E_{\mathrm{pump}}$-field strengths. Although the $\omega_{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). $\textbf{c},$ $E_{\mathrm{pump}}$-field dependence of the calculated
Higgs mode $\omega_{H,1}$ spectral weight without inter-band coupling, i.e.,
$U=0$ (red circles) and with strong inter-band coupling $U\neq 0$ (blue
circles). $\textbf{d},$ Plot of the Higgs mode frequency $2\omega_{H,1}$ as a
function of $E_{\mathrm{pump}}$-field strength without inter-band coupling
(red circles) and for strong inter-band coupling (blue circles). Our
simulations for strong $U$ (blue circles) are in full agreement with the
experimental results in Figs. 3c-3d and in sharp contrast with the one-band SC
results obtained for $U=0$ (red circles). e, f, The calculated spectral weight
SW${}_{0\rightarrow 14~{}\mathrm{meV}}$ 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\neq 0$ (blue circles). With inter-band coupling, the SW
is strongly suppressed, by about 60$\%$ up to a temperature of $0.6~{}T_{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 put the above Higgs mode findings on a rigorous footing, we calculate
two–dimensional THz coherent nonlinear spectra by extending the gauge-
invariant density matrix equations of motion theory of Ref. mootz2020 , as
outlined in the supplement (Section S6). 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 inter-band interaction limits. For this, we
calculate the APS and quantum transport nonlinearities mootz2020  driven by
two intense phase–locked THz E-field pulses for an effective 3-pocket BCS
model of FeSCs fernandes:2016 . 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
$\Delta E_{\mathrm{}}/E_{0}$ for two phase–locked THz pulses, which allows for
a direct comparison of our theory with the experiment (supplementary).

The inset of Fig. 4a presents the calculated $\Delta E/E_{0}$ (black line),
shown together with E${}_{\mathrm{pump}}^{2}(t)$ 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 $\omega_{\mathrm{\mathrm{H},1}}$. This
resonance grows up to pump fields $E_{\mathrm{THz}}\approx 320.0$ kV/cm for
the parameters used here, with minimal redshift and without any significant SW
at $\omega_{\mathrm{\mathrm{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
$\omega_{\mathrm{\mathrm{H},2}}$ Higgs mode. The latter mode is strongly
suppressed due to damping induced via electron-hole asymmetry (supplementary).
Specifically, the ellipticity of the e-pockets increases the DOS along the
pump field direction and thus increases the damping of mainly the
$\omega_{\mathrm{\mathrm{H},2}}$ resonance, which leads to a transfer of
oscillator strength to the continuum. This damping has a much smaller
influence on the $\omega_{\mathrm{\mathrm{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 $\omega_{\mathrm{\mathrm{H},1}}$
resonance amplitude, with SW reduction accompanied by a persisting mode
frequency This behavior of the multi–band 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 $\omega_{\mathrm{\mathrm{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 (Figs. 3c-3e).

To scrutinize further the critical role of the strong interband interaction
$U$, we show the fluence dependence of the coherent Higgs SW close to
$\omega_{\mathrm{\mathrm{H},1}}$ in Fig. 4c. SW${}_{0\rightarrow
14~{}\mathrm{meV}}$ differs markedly between the calculation with strong
$U\neq 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\rightarrow 14~{}\mathrm{meV}}$ for strong $U$ grows at low
pump fluences (regime I), followed by a saturation and then decrease at
elevated $E_{\mathrm{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\neq 0$ (blue
circles) vs. $U=0$ (red circles)). Without inter-pocket $U$ ($U=0$ , red
circles in Figs. 4c-4d), the SW of Higgs mode $\omega_{\mathrm{\mathrm{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 $\Delta_{1}$
and a decrease of the SW of Higgs mode $\omega_{\mathrm{H,1}}$ to zero, due to
ransition from a damped oscillating Higgs phase to an exponential decay. Based
on the calculations in Figs. 4c-4d, the decrease of the spectral weight with
interband coupling appears at a $\Delta_{1}$ quench close to 15$\%$, while
without interband coupling the decrease of spectral weight is only observable
close to the complete quench ($\sim$90$\%$) of the SC order parameter
$\Delta_{1}$. 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_{\mathrm{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 Figs. 4e and 4f for $U=0$ (red circles) and $U\neq 0$ (blue
circles). With interband coupling, the SW is strongly suppressed, by about
60$\%$ up to a temperature of $0.6~{}T_{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 $\omega_{\mathrm{H,1}}$ to the higher mode
$\omega_{\mathrm{H,2}}$ with increasing temperature, since the higher SC gap
$\Delta_{2}$ experiences stronger excitation by the applied pump $E^{2}$ 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, Figs. 4e-4f). Moreover, our calculation without light-induced changes
in the collective effects (only charge-density fluctuations) produces a
significantly smaller $\Delta E_{\mathrm{}}/E_{0}$ signal in the non-
perturbative excitation regime (Fig. S11, supplementary). 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 Higgs modes and coherent
excitations in FeSCs, which differ significantly from any previously observed
collective mode in other superconducting materials: 2$\Delta_{\mathrm{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
a compelling mechanism of hybrid Higgs from the two band entanglement in
FeSCs. Our results also warrant further investigation of Higgs collective
modes through broadband THz 2D spectroscopy in the coherent driven regime.

## References

  * (1) I. I. Mazin, J. Schmalian, Pairing symmetry and pairing state in ferropnictides: Theoretical overview. Physica C: Superconductivity 469, 614-627 (2009).
  * (2) Johnston D C, The puzzle of high temperature superconductivity in layered iron pnictides and chalcogenides. Adv. Phys. 59, 803 (2010).
  * (3) Hanaguri T, Niitaka S, Kuroki K and Takagi H, Unconventional s-Wave Superconductivity in Fe(Se,Te). Science 328, 474 (2010).
  * (4) A. Patz et al., Ultrafast observation of critical nematic fluctuations and giant magnetoelastic coupling in iron pnictides. Nat. Commun. 5, 3229 (2014).
  * (5) A. Patz et al., Critical speeding up of nonequilibrium electronic relaxation near nematic phase transition in unstrained Ba(Fe1-xCox)2As2. Phys. Rev. B 95, 165122 (2017).
  * (6) A Charnukha, Optical conductivity of iron-based superconductors. J. Phys.: Condens. Matter. 26, 253203 (2014).
  * (7) R. Matsunaga et al., Light-induced collective pseudospin precession resonating with Higgs mode in a superconductor. Science 345, 1145–1149 (2014).
  * (8) T. Cea, C. Castellani, L. Benfatto, Nonlinear optical effects and third-harmonic generation in superconductors: Cooper pairs versus Higgs mode contribution. Phys. Rev. B 93, 180507(R) (2016).
  * (9) M. Udina, T. Cea, L. Benfatto, Theory of coherent-oscillations generation in terahertz pump-probe spectroscopy: From phonons to electronic collective modes. Phys. Rev. B 100, 165131 (2019).
  * (10) F. Giorgianni et al., Leggett mode controlled by light pulses. Nat. Phys. 15, 341–346 (2019).
  * (11) R. Matsunaga et al., Higgs amplitude mode in the BCS superconductors Nb1-xTixN induced by terahertz pulse excitation. Phys. Rev. Lett. 111, 057002 (2013).
  * (12) S. Rajasekaran et al., Probing optically silent superfluid stripes in cuprates. Science 359, 575-579 (2018).
  * (13) X. Yang et al., Terahertz-light quantum tuning of a metastable emergent phase hidden by superconductivity. Nat. Mater. 17, 586-591 (2018).
  * (14) X. Yang et al., Lightwave-driven gapless superconductivity and forbidden quantum beats by terahertz symmetry breaking. Nat. Photon. 13, 707-713 (2019).
  * (15) X. Yang et al., Nonequilibrium pair breaking in Ba(Fe1-xCox)2As2 superconductors: evidence for formation of a photoinduced excitonic state. Phys. Rev. Lett. 121, 267001 (2018).
  * (16) C. Vaswani et al., Terahertz Second-Harmonic Generation from Lightwave Acceleration of Symmetry-Breaking Nonlinear Supercurrents. Phys. Rev. Lett. 124, 207003 (2020).
  * (17) F. Yang, M. W. Wu, Gauge-invariant microscopic kinetic theory of superconductivity: Application to the optical response of Nambu-Goldstone and Higgs modes. Phys. Rev. B 100, 104513 (2019).
  * (18) A. Kumar, A. F. Kemper, Higgs oscillations in time-resolved optical conductivity. Phys. Rev. B 100, 174515 (2019).
  * (19) L. Schwarz et al., Classification and characterization of nonequilibrium Higgs modes in unconventional superconductors. Nat. Commun. 11, 287 (2020).
  * (20) H. Chu et al., Phase-resolved Higgs response in superconducting cuprates. Nat Commun 11, 1793 (2020).
  * (21) M. Buzzi et al., Higgs-mediated optical amplification in a non-equilibrium superconductor. Preprint at https://arxiv.org/abs/1908.10879 (2020).
  * (22) R. M. Fernandes, J. Schmalian, Competing order and nature of the pairing state in the iron pnictides. Phys. Rev. B 82, 014521 (2010).
  * (23) A. Akbari, A. P. Schnyder, D. Manske, I. Eremin, Theory of nonequilibrium dynamics of multiband superconductors. Europhys. Lett. 101, 17002 (2013).
  * (24) T. Cea, P. Barone, C. Castellani, L. Benfatto, Polarization dependence of the third-harmonic generation in multiband superconductors. Phys. Rev. B 97, 094516 (2018).
  * (25) T. Cea, L. Benfatto, Signature of the Leggett mode in the A1g Raman response: From MgB2 to iron-based superconductors. Phys. Rev. B 94, 064512 (2016).
  * (26) Y. Murotani, N. Tsuji, H. Aoki, Theory of light-induced resonances with collective Higgs and Leggett modes in multiband superconductors.  Phys. Rev. B 95, 104503 (2017).
  * (27) H. Krull, N. Bittner, G. S. Uhrig, D. Manske, A. P. Schnyder, Coupling of Higgs and Leggett modes in non-equilibrium superconductors. Nat. Commun. 7, 11921 (2016).
  * (28) Blumberg G., et al., Observation of Leggett’s collective mode in a multi-band MgB2 superconductor. Phys. Rev. Lett. 99, 227002 (2007).
  * (29) S. Lee et al., Template engineering of Co-doped BaFe2As2 single-crystal thin films, Nat. Mater. 9, 397 (2010).
  * (30) X. Yang et al., Ultrafast nonthermal terahertz electrodynamics and possible quantum energy transfer in the Nb3Sn superconductor. Phys. Rev. B 99, 094504 (2019).
  * (31) J. J. Tu, et al., Optical properties of the iron arsenic superconductor BaFe1.85Co0.15As2. Phys. Rev. B. 82, 174509 (2010).
  * (32) K.W. Kim, et al., Evidence for multiple superconducting gaps in optimally doped BaFe1.87Co0.13As2 from infrared spectroscopy. Phys. Rev. B. 81, 214508 (2010).
  * (33) T. Maag et al., Coherent cyclotron motion beyond Kohn’s theorem. Nat. Phys. 12, 119 (2016).
  * (34) E. A. Yuzbashyan, M. Dzero, Dynamical vanishing of the order parameter in a fermionic condensate. Phys. Rev. Lett. 96, 230404 (2006).
  * (35) T. Papenkort, V. M. Axt, T. Kuhn, Coherent dynamics and pump-probe spectra of BCS superconductors. Phys. Rev. B 76, 224522 (2007).
  * (36) Y.-Z. Chou, Y. Liao, M. S. Foster, Twisting Anderson pseudospins with light: Quench dynamics in terahertz-pumped BCS superconductors. Phys. Rev. B 95, 104507 (2017).
  * (37) M. Mootz, J. Wang, and I. E. Perakis, Lightwave terahertz quantum manipulation of nonequilibrium superconductor phases and their collective modes, Phys. Rev. B 102, 054517 (2020)
  * (38) R. M. Fernandes, A. V. Chubukov, Low-energy microscopic models for iron-based superconductors: a review. Rep. Prog. Phys. 80, 014503 (2016).

## Acknowledgments

This work was supported by National Science Foundation 1905981 (THz
spectroscopy and data analysis). 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. Theory work at the University of
Alabama, Birmingham was supported by the US Department of Energy under
contract # DE-SC0019137 (M.M and I.E.P) and was made possible in part by a
grant for high performance computing resources and technical support from the
Alabama Supercomputer Authority.

## Author Contributions

C.V., L.L. and X. Y. performed the THz spectroscopy measurements with J.W.’s
supervision. J.H.K, C.S. and C.B.E. grew the samples and performed crystalline
quality and transport characterizations. M.M. and I.E.P. developed the theory
for the hybrid Higgs mode and performed calculations. Y. G. C and E. E. H made
Ba122 target for epitaxial thin films. J.W. and C.V. analyzed the THz data
with the help of L.L., D.C., C.H., R.J.H.K. and Z.L. The paper is written by
J.W., M.M., I.E.P. and C.V. with discussions from all authors. J.W. conceived
and supervised the project.

## Correspondence

Correspondence and requests for materials should be addressed to J.W.
<EMAIL_ADDRESS>jgwang@iastate.edu).