Classification and characterization of nonequilibrium Higgs modes in unconventional superconductors

Recent findings of new Higgs modes in unconventional superconductors require a classification and characterization of the modes allowed by nontrivial gap symmetry. Here we develop a theory for a tailored nonequilibrium quantum quench to excite all possible oscillation symmetries of a superconducting condensate. We show that both a finite momentum transfer and quench symmetry allow for an identification of the resulting Higgs oscillations. These serve as a fingerprint for the ground state gap symmetry. We provide a classification scheme of these oscillations and the quench symmetry based on group theory for the underlying lattice point group. For characterization, analytic calculations as well as full scale numeric simulations of the transient optical response resulting from an excitation by a realistic laser pulse are performed. Our classification of Higgs oscillations allows us to distinguish between different symmetries of the superconducting condensate.


Introduction
The Higgs mode in superconductors is a collective oscillation of the order parameter ∆ with the characteristic frequency of 2∆.It can be understood as a massive excitation along the radial direction in the Mexican hat potential of the free energy (see Figure 1a) [1][2][3][4][5].The charge neutral Higgs mode does not couple to linear optical probes and therefore was expected to be observable only in materials with competing orders [6,7], for which it was measured in Raman experiments [8][9][10].However, an impulsive excitation of Higgs oscillations in nonequilibrium is possible via a nonlinear process by quenching the Mexican hat potential with an ultrafast THz light pulse.Such a quantum quench was demonstrated for the first time in the s-wave superconductor Nb 1−x Ti x N [11][12][13][14].
There are several works indicating that the spectrum of Higgs modes can be more complex if nontrivial gap symmetries are involved.Studies on multiband superconductors like MgB 2 show that the Higgs oscillation spectrum contains Higgs modes for both gaps as well as the Leggett mode, the relative phase mode [15,16].Additional Higgs modes with lower energies, representing oscillations of the gap in different symmetry channels, were proposed for d-wave superconductors under the assumption of a composite pairing interaction [17].Besides first quench-probe experiments on cuprates [18][19][20] a recent experiment on several types of cuprates shows clear fingerprints of a 2∆ Higgs mode and a so far unknown additional mode below 2∆ [21].These findings require both a deeper understanding and a classification and characterization of Higgs modes in nonequilibrium.
So far all descriptions on how to excite Higgs oscillations with a quench pulse are working within the dipole approximation, i.e. neglecting the small wave momentum q of the external field.There are other studies which show that a linear coupling of the vector potential to the condensate is possible if momentum transfer is involved, either by impurity scattering in dirty superconductors [22][23][24][25][26] or in current-carrying states [27,28].
Independent of the actual coupling to the external field, the following instructive picture can be drawn to understand the excitation process by an ultrashort THz pulse, where Higgs oscillations are excited by taking the superconductor out of equilibrium [15,[29][30][31][32][33][34][35][36][37][38][39].Hereby, Cooper pairs are partially broken and the landscape of the free energy changes suddenly, i.e. the Mexican hat shrinks.Thus the THz laser pulse acts like a quantum quench [40][41][42][43][44][45], reflecting the impulsive character of the light pulse.As long as this process is faster than the time scale of the condensate, given by τ ∆ = h/(2∆), where 2∆ is the energy gap of the superconductor, the condensate is unable to follow the minimum of free energy adiabatically [15,36].Consequently, collective Higgs oscillations of the gap are excited, as sketched in Figure 1a.
In order to excite Higgs oscillations, the laser pulse must fulfill two conditions.On the one hand, the pulse should only excite a small fraction of the Cooper pairs, enough to generate a significant quench of the Mexican hat, but not too many that the superconducting signatures would be screened by hot electrons.More specifically, a short optical pulse far above gap frequencies induces a strong Drude-peak in the optical conductivity, which would overlap with the weak signal of the Higgs oscillations.Instead a suitable pulse corresponds to a peak located in or close to the gap, which only slightly overlaps with the continuum of quasi-particles, as depicted in Figure 1c.On the other hand, the pulse must fulfill the nonadiabaticity condition, which implies a short laser pulse (Figure 1d) and hence requires the broad spectrum in energy domain (rather than a narrow band multicycle pulse tuned close to the gap).For typical gaps in the meV regime, a single cycle THz laser pulse is exactly on the brink of these two regimes [46], allowing for an excitation of the Higgs oscillations without heating or photo-doping the system too much.
In this article we show how in a nonequilibrium setup for superconductors with pairing interaction in a single channel, e.g.pure d-wave, oscillations of the condensate in other symmetry channels can lead to additional Higgs modes as well.We classify these oscillations of the condensate based on the irreducible representations of the point group of the underlying lattice.The resulting Higgs modes depend on the excitation symmetry and the ground state symmetry.Our detailed analysis shows that a full description of the excitation process requires to go beyond the dipole approximation in a Raman-like process and to retain the wave momentum q (see Figure 1b), which plays a crucial role in breaking the symmetry of the condensate and exciting non-A 1g oscillations of the superconducting condensate.We show that nonequilibrium Higgs oscillations offer a unique way to investigate the symmetry and collective excitation spectrum of superconductors, which allows to completely characterize the nature of a superconducting condensate with a single class of experiments.

Results
Quantum quenches.While the nonequilibrium probe of collective excitations in conventional s-wave superconductors has been studied intensively [29-31, 33, 34, 47-50], the response for unconventional superconductors is still in its infancy.These systems often exhibit very complicated correlations [51,52] and a variety of different mechanisms which can lead to superconductivity.If we want to examine the nonequilibrium response of the coherent condensate of such superconductors in general, we have to go back to the fundamental properties of these systems: The symmetry of the lattice.
According to group theory, every configuration of the condensate at a given time can always be decomposed with respect to the different irreducible representations of the point group symmetry of the lattice.As an example, we take a superconductor on a lattice with D 4h space group symmetry, which is the lattice symmetry of cuprate high-T c superconductors.Based on this argument, there are four different irreducible representations: A 1g , A 2g , B 1g and B 2g .The condensate oscillations can always be decomposed into the contributions from these sectors.
We start our theoretical description by considering a quench of the initial state.Every quantum quench deforms the condensate from its equilibrium value, which then can be decomposed into contributions from different irreducible representations.Taking a momentum independent quench for example would only probe the A 1g channel of the condensate [45], but does not couple to other allowed symmetries.The solution to address also the other possible symmetries is to modify the quench and make it momentum dependent, so that we can probe other symmetry sectors as well.For example in case of  In order to illustrate this concept, we perform numerical simulations for s-and d-wave BCS superconductors to study the nonequilibrium response to momentumdependent quantum quenches.The Hamiltonian we are investigating is given by with f k describing the symmetry of the interaction, V the interaction strength, and the normal state Hamiltonian H 0 = kσ k c † kσ c kσ , where c † kσ creates electrons with momentum k and spin σ.Within the BCS solution, the gap is determined by Now we perform a quantum quench by changing the symmetry of the condensate c −k↓ c k↑ slightly away from its equilibrium value with f k = f k + δf q k , where f q k is the quench symmetry and δ the quench strength.The equilibrium value of the condensate has the symmetry of the gap f k which is not changed and always remains in a single symmetry sector.After the quench we calculate the Higgs oscillation of the order parameter as a function of time by evaluating the time-dependent gap equation (2), which sums the oscillations of the condensate in momentum space.For the temporal evolution we use Anderson pseudospin formalism [53], where the time-evolution is governed by Bloch equations.More details are given in the methods.
For a given gap symmetry, there are different oscillations possible for the condensate, which can be excited depending on the symmetry of the quench.We use a new notation to describe this oscillation symmetry, which takes the gap symmetry into account.We add as an additional information the gap symmetry as subscript to the group theoretic notation of the irreducible representation name.We observe that depending on gap and quench symmetries not only the well known 2∆ Higgs mode occurs in the spectrum of the Higgs oscillation, which appears independent on the quench due to coupling between the modes, but also a second mode at lower energy.
Two examples for this observation are shown in Figure 3.In Figure 3a, we see the Higgs oscillations of the d x 2 −y 2 gap after a f q k ∼ x 2 − y 2 and f q k ∼ 1 quench which excites A 1g x 2 −y 2 or B 1g x 2 −y 2 oscillations of the condensate.The f q k ∼ 1 quench for a s-wave superconductor is shown in red for comparison.We highlight that the d-wave oscillations decay much faster than the s-wave oscillations.This can be traced back to the stronger dephasing of the mode due to coupling to the gapless quasi-particles in the d-wave case.The final value of the gap, i.e. ∆ ∞ , depends on the strength of the quench, i.e. how strongly the initial states deviate from the equilibrium state [42,43].Figure 3b shows the Fourier transform of the Higgs oscillations.The large peak at 2∆ ∞ corresponds to the symmetric A 1g x 2 −y 2 oscillation of the condensate.
Most importantly a second low energy mode is visible for the d-wave superconductor after the f q k ∼ 1 quench.This mode does not exist for pure s-wave superconductors and it is also not excitable by the f q k ∼ x 2 − y 2 quench in the d-wave case, as the quench has the same symmetry as the ground state gap.Similarly for other combinations of gap and quench symmetries additional modes can be identified.Thus there exists a direct connection between the symmetry of the gap and the existence of low energy Higgs modes.
To understand the nature of the second mode in more detail, we perform a linear analysis of the gap dynamics after a quantum quench.Specifically, we analytically compute the dynamics of the gap according to the expansion in the first order of δ∆(t) for different initial states.Here ∆(0) is the gap at time t = 0 directly after the quantum quench.Transforming into Laplace space with complex frequency s allows us to identify the leading contributions ) quench as a function of time.The red solid line shows an f q k ∼ 1 quench for an s-wave superconductor for comparison.b) Fourier spectrum |∆(ω)| = | FT |∆(t)|| of the Higgs oscillations.The oscillation for the d-wave gap excited by the f q k ∼ x 2 − y 2 quench shows a single peak, similar to the s-wave case.The peak position corresponds to 2∆∞ which is two times the maximum of the gap for t → ∞ after the quench.For the f q k ∼ 1 quench a second peak at low energy appears resulting from B 1g x 2 −y 2 oscillations of the condensate.
to the gap dynamics and leads to the expression with where the dots imply additional weighting factors.For F 1 (s), we find the same denominator in the integrand.
In these expressions we assume that the symmetry functions primarily depend on the angle ϕ between k and the k x -axis.For further details on the calculation, please see Supplementary Note 1.We can see that the spectrum of Higgs oscillations is controlled in a nontrivial way by f (ϕ) 2 − f (ϕ) 2 integrated over ϕ weighted by additional factors.Thus we can trace back the second mode to the difference in the symmetry between the quench and the condensate.A second mode in the Higgs oscillation spectrum is only visible if this difference leads to a second minimum in the denominator.Particularly, this happens if there are changes of the nodal directions in the condensate symmetry compared to the equilibrium value.Due to the nonexistence of nodes in the s-wave case, there will be no second mode visible in the Higgs oscillations for all possible condensate oscillations.
Realistic pulse.So far we concentrated our analysis on quantum quenches, which we classified according to the deformation symmetry from the equilibrium value.To show that these results also carry over to more realistic scenarios, we calculate the response of a d x 2 −y 2 -wave BCS superconductor coupled to a laser field.So called pump-probe experiments have been used to study the excitation and relaxation processes of superconductors [51,54].In a pump-probe experiment the pump pulse excites the system and after a delay time the probe pulse measures various properties of the transient dynamics of the system.Varying the delay time, the temporal evolution of the system after a perturbation can be studied.
As the purpose of the pump pulse in our setup is to quench the system suddenly, we call this pulse a quench pulse.The Hamiltonian describing the quench and probe pulses is given in the methods.We use the density matrix formalism [29] to calculate the dynamics, which is exact for the Hamiltonian in (1).We assume a short and intense THz laser pulse which excites the condensate in an anisotropic fashion and the superconducting gap starts to oscillate.For all pulses we fixed the pulse duration to τ p = 0.4 ps.With this choice we are in the nonadiabatic regime, where a generation of collective modes is possible.
Further, we varied the direction of the light wave vector q to study the dependence of the optical conductivity on the quench pulse.Thus we define the angle φ between q and the k x -axis of the superconductor.The light wave vector is small compared to the Fermi wave vector |q| |k F | such that there is no excitation of Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) oscillations.However, it is large enough to break the condensate symmetry as it couples offdiagonal elements in the quasi-particle distribution (see Equation (19)).This is possible due to the Raman-like excitation, where the photon momentum can be transferred to the condensate.By choosing the angle of the quench pulse, different oscillation symmetries can be addressed selectively (see Table I).
We compare both methods, quantum quench and quench pulse, in the Supplementary Figure 1.Besides the time evolution of the gap, the quench-probe optical conductivity provides an experimental fingerprint to observe Higgs oscillations as well.This was explicitly shown in case of s-wave symmetry, i.e. as the oscillation of the conductivity depending on the delay time [33].Thus we use a quench pulse to induce the Higgs oscillations and a much weaker probe pulse in the same direction to measure the optical conductivity.
In Figure 4 the real part of the optical conductivity Re σ(∆t, ω) is shown for the angles φ = 0, along the antinodal direction, and φ = π/4, along the nodal direction, as a function of frequency ω and time delay ∆t.These angles correspond to the pulse direction with maximum response in Table I.
For φ = π/4 the symmetry breaking happens along the diagonal axis which excites in addition to the symmetric A 1g x 2 −y 2 also the B 2g x 2 −y 2 oscillation of the condensate and most of the weight is located at the energy of the 2∆ Higgs mode.There is no low-lying peak visible in the spectrum as the B 2g x 2 −y 2 oscillation does not lead to a second mode.For φ = 0 the A 1g x 2 −y 2 and B 1g x 2 −y 2 oscillations are excited, resulting in the low energy Higgs mode and the 2∆ Higgs mode and the spectrum is dominated by an in-gap response.Most importantly the signal oscillates with respect to the time delay between quench and probe pulse, reflecting the excitation of Higgs oscillations.This is demonstrated in Supplementary Figure 3 in more detail.Thus the angle-resolved quench-probe experiments should be able to see Higgs oscillations in the optical conductiv-ity, which can address different oscillation symmetries of the condensate depending on the incident angle.
Classification of Higgs oscillations.Using the results from the quench and quench-probe calculations, an extension to other symmetries and evaluation of the Higgs oscillations enables us to write down a classification scheme for nonequilibrium Higgs modes.As an example, for the D 4h lattice point group, the condensate oscillations and resulting Higgs modes as well as their excitation symmetries are shown in Table I for all fundamental gap symmetries allowed by the point group.
For each gap symmetry f k in column one, quench symmetries f q k in all irreducible representations of the point group are listed in column two.Quenching via pump pulses at an arbitrary incident angle φ results in an excitation of all modes.Choosing a direction along the high symmetry directions φ = 0 or φ = π/4, a selective excitation of B 1g or B 2g oscillations is possible.These directions, which correspond to the respective quench symmetry with maximum intensity are listed in the third column.As the light pulse always breaks the symmetry it is in principle not possible to excite the symmetric A 1g mode alone, even in the s-wave case.Induced by the quench, the resulting oscillations of the condensate c −k↓ c k↑ are shown in column four.Independent of the quench symmetry, the symmetric A 1g oscillation is always excited as any disturbance leads to a global change in the quasi-particle distribution.The time-dependent amplitude of the energy gap is calculated from the gap equation (2), where the condensate for each momentum Table I.Classification of Higgs oscillations Possible Higgs oscillations for a lattice with D 4h point group symmetry shown for s, d x 2 −y 2 , dxy and g xy(x 2 −y 2 ) gap functions (column one).A quench can be applied to the condensate with a certain symmetry f q k (column two), which disturbs the ground state condensate symmetry.These quenches can be controlled by an incident THz pulse with angle φ.Pumping at an arbitrary angle corresponds to a quench in all symmetry channels.Choosing high symmetry direction (column three) allows for a selective excitation.Such a quench excites oscillations of the condensate (column four), classified by the notation of the irreducible representations of the lattice symmetry.Oscillations of the condensate lead to amplitude oscillations of the gap and the qualitative Fourier spectrum of these Higgs oscillations is illustrated in the last column, showing the possible Higgs modes.An animation on how each quench deforms a given symmetry can be found in the Supplementary Movie 1.
point is summed.This results in Higgs oscillations of the gap and a schematic picture of the spectrum is shown in the last column.For all gap symmetries the 2∆ Higgs mode is visible, which corresponds to the symmetric A 1g oscillation of the condensate.Depending on how the non-A 1g oscillations change the condensate symmetry from its equilibrium symmetry, i.e. the gap symmetry, a second Higgs mode is visible in the spectrum.This is not always the case, e.g. the B 2g x 2 −y 2 oscillation is not visible in the spectrum as a second Higgs mode despite its asymmetric deviation from the ground state symmetry.As this oscillation only shifts weight inside the positive and negative lobes of the d x 2 −y 2 symmetry but doesn't move the nodal directions, it will not lead to a second mode in the summation process for the calculation of the gap oscillations.Hence, for a full analysis of the gap symmetry information from multiple quench symmetries are required.Yet, if we can obtain this information, nonequilibrium Higgs oscillations can be used as an efficient tool to completely classify the ground state symmetry of a superconductor.

Discussion
To summarize, we introduce a classification scheme for nonequilibrium Higgs oscillations which allows to characterize the ground state of superconducting condensates.Our analytical calculations show that depending on the symmetry of the quench and of the gap function, lowlying modes exist, which can be directly identified with the different oscillations of the condensate.We introduce a new notation to combine the information of the ground state with the quench symmetry in order to distinguish the different Higgs oscillations.Simulations of quenchprobe experiments using realistic pulses in a microscopic model show that the usually ignored wave momentum in the dipole approximation plays an important role in the excitation of non-A 1g oscillations of the condensate.Despite its small value compared to the Fermi wave vector, it is large enough to break the ground state symmetry and can lead to additional Higgs modes implementing the proposed analytic quench setup.
It is important to note that the proposed experimental excitation of the Higgs mode is a Raman-like excitation and should not be confused with an infrared-active excitation [28].In the latter case, a driven ac current would occur and thus the strength and polarization of the electric field is more important than the small momentum of the photon.In the former case, the polarization of the electric field plays a minor role and the photon momentum becomes much more important.
We find that the Higgs modes are visible in the optical conductivity of the proposed quench-probe experiment, paving the way for investigations and classifications of the dynamics of known and unknown superconductors directly within this framework.This analysis is applicable for all superconductors and requires only the knowledge of the symmetry of the crystal.It is a natural extension of the group theoretical notation to the case of nonequilibrium excitation of the system.To demonstrate the approach, we fully characterize all possible Higgs oscillations for the important D 4h point group, relevant for example for high-temperature cuprate superconductors.This main result is summarized in Table I, which goes beyond a simple product table of gap symmetry and light pulse direction.The number of excited fundamental condensate oscillations does not directly correlate with the number of observed Higgs modes, which depend in a nontrivial way on the phase and nodal structure of the order parameter.
The experimental realization of the proposed momentum transfer to break the condensate symmetry will be challenging.If light couples to the Higgs mode only indirectly via electrons, momentum scattering on timescales faster than the oscillation period might wash out or destroy the preferred direction of the pulse.Depending on the strength of this momentum distribution effect, the second Higgs mode could be damped, might no longer follow its predicted angular dependency or may be even completely suppressed.On the other hand, recent studies have shown that impurity scattering in dirty superconductors even enhances the coupling of light to the condensate and the excitation of the Higgs mode [22][23][24][25][26].As no experiments exist so far which allow to measure the transferred photon momentum in detail, the field is open for further experimental and theoretical investigations.However, the current efforts to measure Higgs oscillations on different cuprates [18][19][20][21] show already fingerprints of collective oscillations.
It is important to note that we do not introduce additional energy scales nor other degrees of freedom, such as subdominant channels: the observed oscillations and the corresponding frequencies are intrinsic to the pure d-wave superconductor and do not require composite pairing symmetries.Note that we assume that no competing order, such as a charge density wave, exists.Otherwise, the spectroscopic signatures of the Higgs oscillations could be modified due to the interplay between the two phases [55].Furthermore, effects which can modify the Higgs spectrum as well are superconductors in the strongly coupled regime [37], coupling to Leggett modes in multiband systems [15] or collective excitations of pair states in subleading channels, i.e. an excitation of the Bardasis-Schrieffer mode [56][57][58][59].Other details of the normal state, such as Fermi arcs, play an unimportant role after a quantum quench of the superconducting condensate, as they only modify the scattering processes of the broken Cooper pairs.This could potentially change the damping of the Higgs oscillations, but has no effect on our proposed classification scheme.
There are different possibilities how a symmetry breaking momentum transfer to the condensate is realized in an experiment.Tilting the quench pulse direction towards the superconducting plane induces a finite in-plane momentum of the photons.More controlled momentum dependent excitations and probes are possible in a THzfour-wave mixing [60] or transient grating [61,62] setup.Other possibilities include momentum-dependent scattering processes as well as coupling to other finite-momentum modes.This is discussed for example for superconductors under external current [27,28] or as a possibility for phonon-coupled amplitudon dynamics in excitonic insulators [63,64].
The classification and characterization of Higgs oscillations open the possibility to perform spectroscopic studies on superconductors to determine the symmetry of the order parameter.Compared to other types of measurements like ARPES or interferometry experiments using Josephson junctions, which can retrieve either amplitude or phase information, spectroscopy of Higgs oscillations with phase-stable THz lasers allows to determine amplitude and phase within a single type of quench-probe experiment.In principle our theory for Higgs spectroscopy is not limited to characterizing equilibrium condensates but may also be used to investigate possible light induced superconducting states in transient states of matter [65][66][67][68].Beyond that, Higgs spectroscopy could also be extended to investigate collective excitations of non-superconducting, symmetry broken phases, such as order parameter oscillations in excitonic insulators [63,64] or Higgs modes in antiferromagnets [69].

Methods
Extended BCS model.The Hamiltonian we are investigating is given by The normal state Hamiltonian H 0 is taken to be a free electron gas with an effective mass m.The pairing interaction V kk = V f k f k is assumed to be separable with the interaction strength V .The energy dispersion k = 2 k 2 /(2m) − F is measured relative to the Fermi level F .We apply the BCS solution in order to describe the superconducting phase.The superconducting gap equation reads The sums in Equations ( 7) and ( 9) are taken over the set W of all k vectors with | k | ≤ ω c , ω c being the frequency cutoff.For a phononic glue, this corresponds to the Debye frequency.The function f k is the gap symmetry function, where in the case of an s-wave superconductor f k = 1 is a constant.In general, the symmetry function can be decomposed into the basis functions of the irreducible representations of the point group of the underlying lattice.In case of the D 4h group, the basis functions are shown in Supplementary Table 1, where the d-wave symmetry belongs to the B 1g representation.For all of our calculations we assume that there is only a polar angle dependency ϕ on the momentum k in the vicinity of the Fermi energy.Therefore we use the functions shown in the third column.
In the ground state the expectation values for the electron and quasi-particle distribution read where Anderson pseudospin description.We define the Nambu-Gorkov spinor and Anderson pseudospin [53] where τ are the Pauli matrices.The BCS Hamiltonian takes the form with where we assume a fixed phase of the gap such that ∆ ∈ R.
In equilibrium, the y-component of the pseudospin is zero σ y k = 0, whereas the x-and z-component read At t = 0, we apply a state quench where we change the symmetry of the condensate by changing the pseudospin expectation values where k with the quench symmetry f q k and strength δ.This changes the initial ground state symmetry of the condensate, which is the same as the gap symmetry, to the quenched symmetry f k .Note that the gap ∆(0) at t = 0 is different from the equilibrium gap ∆ due to the sudden change of the system state.For arbitrary times, the gap equation reads This renders the Hamiltonian time-dependent as the pseudomagnetic field depends on the gap.The time-evolution of the pseudospin in the quenched system is described by Bloch equations [47] The Bloch equations (18) can then be solved together with the time-dependent gap equation ( 17) self-consistently.
Coupling to vector potential.The Hamiltonian describing the coupling between superconductor and quench pulse, which brings the system out of equilibrium, is modeled by where A q (t) is the transverse vector potential [29,33].
Working within the Coulomb gauge, the quench pulse is expressed in terms of the transverse vector potential The quench pulse is of Gaussian shape with photon frequency ω p , photon wave vector q p , full width at half maximum (FWHM) τ p and amplitude A p .For our simulations we consider various directions of the photon wave vector q p and with this concomitantly various directions of the quench induced by the pulse.
Optical conductivity.To calculate the optical conductivity, we calculate the temporal evolution of the current density as function of the time delay between the quench and probe pulse where V is the normalization volume and q pr is the wavevector of the probe pulse [29].We neglect the second term, because it only leads to an offset of the imaginary part of the optical conductivity.Then, the optical conductivity can be calculated [33] by computing σ(∆t, ω) = j q pr (∆t, ω) iωA q pr (ω) .
Density matrix formalism.In order to simulate the evolution of the system, we use methods based on an expansion of Heisenberg's equation of motion.For the temporal evolution of the order parameter we use the density matrix formalism [70].The main task of this technique is to derive equations of motion for quasi-particle densities.Within this formalism it is advantageous to perform a Bogoliubov transformation of the electron operators which diagonalizes the Hamiltonian H BCS in the initial state.We introduce new fermionic operators α k and β k , with where We emphasize that the coefficients u k and v k do not depend explicitly on time, i.e., the temporal evolution of the quasi-particle densities is computed with respect to a fixed time-independent Bogoliubov-de Gennes basis in which the initial state is diagonal.All physical observables, such as the order parameter amplitude |∆ k (t)| can now be expressed in terms of the new Bogoliubov quasi-particle densities α and α k β k .Applying the density matrix formalism for these quasi-particle densities, we get a closed set of differential equations.The ensuing differential equations are then solved on a finite size grid in momentum space.More details about the implementation can be found in [29,33].
The pulse solution for a pumping angle of φ = 0 is shown in Supplementary Figure 1.It shows a stronger broadening than the analytical calculations, but is in qualitative agreement.Note that we do not expect quantitative agreement, since a quantum quench is different from a laser pulse.
Numerical implementation.In our simulations we use the parameters ∆ = 1.35 meV, E F = 9470 meV and m = 1.9 m e , which are motivated by the parameters for lead [29].However all our computations can be rescaled to any energy scale for the gap.The numerical equations are computed on a finite size grid in momentum space in 2 dimensions similar to [29,33].To obtain the required accuracy to resolve the small wave momentum q p , we restrict our grid in a small region around the Fermi energy with a cutoff of E c = 8.3 meV.We have ensured by varying the cutoff energy that our results do not depend on the discretization range.The x-direction is discretized with a step size of the wave momentum q p , which results in 1000 -2000 points.This is advantageous as we can resolve directly the coupling between the offdiagonal elements like α k β k+q .As the coupling in y-direction is only indirect via the energy gap and therefore much smaller, we choose between 100 and 500 points for this direction.In total we have of the order of 10 6 grid points.To reduce computational effort, we consider offdiagonal elements like α k β k+nq only up to n = 4 as larger offdiagonal elements only contribute in order O A 5 p .We focus first on the nominator and evaluate the integral × tanh −1 s 2 − 4(∆ 2 f (ϕ) 2 − ∆(0) 2 f (ϕ) 2 ) The expression for F 1 (s) reads Without explicitly solving the integral over ϕ, we already see that possible modes, i.e. maxima in the function F 2 (s), are controlled by the difference f (ϕ) 2 − f (ϕ) 2 .This expression in the denominator determines the main oscillation frequencies in combination with the other weighting factors in the integrand.Depending on the gap and the quench symmetry, a low energy peak may appear as observed in the main text.

SUPPLEMENTARY NOTE 2: COMPARISON BETWEEN ANALYTICAL CALCULATIONS AND NUMERICAL RESULTS
Since we cannot solve the integral (10) analytically, we calculate the gap dynamics of the linearized equations numerically and compare it to the solution of the full equation in Supplementary Figure 1.We assume that the order parameter has d-wave symmetry, i.e. f (ϕ) = cos(2ϕ) according to Supplementary Table I, and we quench in the A 1g channel, i.e. f q (ϕ) = 1.Both solutions show a low energy peak, however the full nonlinear equations further amplify the intensity of the low energy mode.The energy of the lower peak depends on the quench strength, as a stronger quench increases the difference between the ground state symmetry f (ϕ) and the quenched symmetry f (ϕ).To classify all possible condensate oscillations, we quench all fundamental gap symmetries with all fundamental quenches allowed by the D 4h point group from Supplementary Table I.In Supplementary Figure 5, we show the quenched symmetry functions for all combinations of gap and quench symmetry.After the quench, we perform a time-evolution, where we extract for each step of the time-evolution the symmetry of the condensate c −k↓ c k↑ for each step of the time-evolution and classify its oscillation.According to Equation (10) of the main text, the expectation value in equilibrium reads Using our assumption that k = (|k|) and f k = f (ϕ), for k-values far away from k F it follows We therefore use c −k↓ c k↑ at a k-value far from k F as an approximation to trace the dynamics of the symmetry of the condensate.An animation of the condensate dynamics can be found in the Supplementary Movie 1.We name the oscillations according to its symmetry.Finally, we calculate the gap oscillation and its Fourier spectrum to obtain the collective Higgs modes resulting from the condensate oscillations.

Figure 1 .
Figure 1.Illustration of Higgs oscillations in a superconductor.a) Free energy landscape F of a superconductor as a function of the real and imaginary part of the superconducting gap ∆.After a quench at t1 > t0, the free energy is suddenly changed, exciting the superconducting condensate and leading to collective Higgs oscillations, indicated by a black arrow.The red arrow indicates a quench by a THz light pulse.b) Feynman diagram describing the excitation of a Higgs mode H by a light field A using the Raman vertex.An infrared excitation of the Higgs mode (not considered here) is only possible if an external current is present.c) Higgs excitation mechanism using a THz quench pulse.The quench pulse only slightly overlaps with the quasi-particle continuum indicated in blue.The Mexican hat shrinks due to the breaking of Cooper pairs.d) To excite the Higgs oscillation the pulse must fulfill the nonadiabaticity condition in time domain.

Figure 2 .
Figure 2. Illustration of d-wave condensate oscillation symmetries.Possible condensate oscillation symmetries for a d x 2 −y 2 -wave superconductor with point group symmetry D 4h of the underlying lattice.The arrows indicate the motion of the lobes as a function of time.The notation of the gap symmetry in the subscript stresses the initial state, from which the oscillations of the condensate are excited.

d x 2
−y 2 -wave superconductivity, the possible oscillations of the condensate are shown in Figure2.

Figure 3 .
Figure 3. Higgs oscillations of a d-wave superconductor.a) Numerical simulation of the Higgs oscillations induced by various quench symmetries for a d x 2 −y 2 -wave superconductor.The solid (dotted) blue line shows the gap oscillations after af q k ∼ 1 (f q k ∼ x 2 − y 2) quench as a function of time.The red solid line shows an f q k ∼ 1 quench for an s-wave superconductor for comparison.b) Fourier spectrum |∆(ω)| = | FT |∆(t)|| of the Higgs oscillations.The oscillation for the d-wave gap excited by the f q k ∼ x 2 − y 2 quench shows a single peak, similar to the s-wave case.The peak position corresponds to 2∆∞ which is two times the maximum of the gap for t → ∞ after the quench.For the f q k ∼ 1 quench a second peak at low energy appears resulting from B 1gx 2 −y 2 oscillations of the condensate.

Figure 4 .
Figure 4. Optical conductivity of a d-wave superconductor after excitation with a quench pulse.Real part of the optical conductivity Re σ(∆t, ω) after a realistic quench pulse for incident angles a) φ = 0 and b) φ = π/4 for a d x 2 −y 2 -wave superconductor.The pulse parameters are τp = 0.4 ps, |Ap| = 7 • 10 −8 Js C −1 m −1 and ω = 3 meV for the quench pulse and τp = 0.25 ps, |Ap| = 1 • 10 −8 Js C −1 m −1and ω = 2.5 meV for the probe pulse.The gap value in the simulation is 2∆ = 2.7 meV.The vertical axis denotes the time delay between the excitation of the system with the quench pulse and the probe pulse.The horizontal axis denotes frequency ω.The oscillation frequencies in ∆t correspond to the frequencies of the Higgs modes as shown in Supplementary Figure3.

Supplementary Figure 1 :Supplementary Figure 2 :
Comparison between different solutions for d-wave Higgs modes Fourier spectrum |∆(ω)| = | FT |∆(t)|| of the gap oscillations induced by a quench f q (ϕ) = 1 and a quench pulse.The dynamics in the quenched case is calculated once in the linear approximation and once with the full solution.All three results show the 2∆ Higgs mode and a second mode below.Comparison between experimental pulse and theoretical parametrization.The electrical field of a single cycle THz quench pulse as used in experiment.The parameters of a fit with a Gaussian shape, as given in Equation (20) of the main text, are given by ωp = 1.33 meV and τp = 1.63 ps.
Supplementary TableI: Basis functions for D 4h point group Even parity basis functions for the D 4h point group in Cartesian coordinates and as function of the azimuthal angle ϕ in polar coordinates[2].