Momentum-resolved analysis of condensate dynamic and Higgs oscillations in quenched superconductors with tr-ARPES

Higgs oscillations in nonequilibrium superconductors provide an unique tool to obtain information about the underlying order parameter. Several properties like the absolute value, existence of multiple gaps and the symmetry of the order parameter can be encoded in the Higgs oscillation spectrum. Most studies so far concentrate on experiments in which momentum averaged quantities like the optical conductivity or third-harmonic effects in the transmitted light field are investigated, which allows to access the collective Higgs dynamic but does not give access to additional information contained in the oscillation spectrum of the underlying condensate. Here, we study the time-resolved spectral function measured in angle-resolved photo emission spectroscopy for different quench protocols. We analyze the induced amplitude oscillations in the whole Brillouin zone to understand the creation process of collective Higgs oscillations. These momentum-dependent amplitude oscillations of the spectral function provide more information than the usual considered oscillations in energy and an evaluation of the oscillation phase allows to trace the collective condensate movement.

Higgs oscillations in nonequilibrium superconductors provide an unique tool to obtain information about the underlying order parameter. Several properties like the absolute value, existence of multiple gaps and the symmetry of the order parameter can be encoded in the Higgs oscillation spectrum. Most studies so far concentrate on experiments in which momentum averaged quantities like the optical conductivity or third-harmonic effects in the transmitted light field are investigated, which allows to access the collective Higgs dynamic but does not give access to additional information contained in the oscillation spectrum of the underlying condensate. Here, we study the time-resolved spectral function measured in angle-resolved photo emission spectroscopy for different quench protocols. We analyze the induced amplitude oscillations in the whole Brillouin zone to understand the creation process of collective Higgs oscillations. These momentum-dependent amplitude oscillations of the spectral function provide more information than the usual considered oscillations in energy and an evaluation of the oscillation phase allows to trace the collective condensate movement.

I. INTRODUCTION
Angle-resolved photo emission spectroscopy (ARPES) is a powerful method as it allows to measure the electronic structure of materials directly [1]. Combined with additional time resolution (tr-ARPES), dynamic processes of systems out of equilibrium can be studied in great detail [2]. Hereby, a pump pulse excites the system in a nonequilibrium state prior to the photo emission probe pulse, which is applied after a variable time delay to scan the dynamic of the system.
In recent years, there was an evolving interest in studying collective excitations of materials in nonequilibrium as these can provide fundamental insights into the internal symmetries and properties of a system. This has been demonstrated on charge-density wave systems [3][4][5] and high-temperature cuprate superconductors [6][7][8]. In superconductors, an intrinsic collective amplitude (Higgs) mode of the order parameter exists due to spontaneous U (1) symmetry breaking at the transition to the superconducting state [9,10]. An excitation of the Higgs mode is challenging with the result that only in the last years, an exploration of this area started. This progress was heavily supported by the upcoming of new technologies like ultrafast THz spectroscopy [11], which possess the required excitation energy to not fully deplete the superconducting condensate in the pump process.
The main reason for the difficulties to excite the Higgs mode is due to the fact that the Higgs mode is a scalar mode with neither net charge nor electric dipole, which requires to go beyond the linear excitation regime. While the first predictions of the Higgs mode in superconductors are relatively old [12,13], only indirect measurements in systems with competing charge-density wave orders could be realized in the beginning [14]. The first direct observation was performed in a THz pump-probe experiment on the s-wave superconductor Nb 1−x Ti x N [15], where Higgs oscillations of the order parameter, reflected in oscillations of the electromagnetic response, could be observed. Since then, only a few more experiments were reported [16][17][18][19][20], but theoretical works on the subject became more popular as it became clear that Higgs oscillations in nonequilibrium systems provide rich information about the superconducting ground state.
It was found that much information about the ground state order parameter is encoded in the Higgs oscillation frequency. First of all, the main Higgs oscillation frequency corresponds to the value of the order parameter itself [21][22][23][24][25]. In addition, multiband systems show frequencies for each band and also for the frequency of the relative phase mode (Leggett mode) [26,27]. Excited in an asymmetry way, nontrivial gap symmetries can show additional oscillation frequencies as a result of an oscillation of the condensate in a different symmetry channel [28]. Finally, composite order parameters [29], excitations of subleading pairing channels [30] as well as coupling to other coexisting modes [20] might show up as additional frequencies.
As of now, two main approaches to study Higgs oscillations exist, namely measurements of the optical conductivity in pump-probe experiments [31][32][33] and resonances in third-harmonic generation in periodically driven systems [34,35], where both methods allow to deduce the intrinsic Higgs modes. The only different approach, also in a pump-probe setup, was the prediction of Higgs oscillations in tr-ARPES [36][37][38], where the spectral function shows oscillations in energy with the same frequency as the Higgs mode.
In this article, we study the spectral function measured in tr-ARPES in a more general approach. While previous papers concentrated on the oscillation of the spectral function in energy, we also evaluate the amplitude oscillation, which contains more information about the condensate dynamic. We allow arbitrary gap symmetry and study the effect of quenches in symmetry channels different from the ground state. The idea behind such quenches is to model the net effect of pump pulses in a controlled way, which act on the condensate momentum-dependently. These might be realized experimentally by tuning polarization and pulse direction [28] or could be implemented in more complex approaches like transient grating [39,40] or four-wave mixing [41] setups.
Quenching superconductors with such an approach, where the symmetry of the condensate is altered with respect to the ground state, can result in dynamically created additional Higgs modes [28]. It is important to note that we ensure that the order parameter always keeps its ground state symmetry and only the underlying condensate is modified. This is plausible as a modification of the gap symmetry and thus, the pairing interaction itself, in a conventional pump-probe experiment is unlikely. However, a modification of the condensate, i.e. the Cooper pair or quasiparticle distribution might be controllable with light excitation.
By evaluating the induced oscillation of the spectral function, both in amplitude and energy, we explore the creation process of additional Higgs modes. Moreover, we compare superconductors with the same nodal gap structure but different signs, like d-wave and nodal s-wave, quenched in the same symmetry channel. A calculation of the oscillation phase in different lobes shows opposite phase oscillations, reflecting the different collectively induced condensate movements.

A. Hamiltonian
We consider the mean-field BCS Hamiltonian where k is the electron dispersion relative to the Fermi level and c † kσ , c kσ the creation and annihilation operators for electrons with momentum k and spin σ. The momentum-dependent energy gap ∆ k is defined by where the pairing interaction V kk = V f k f k is assumed to be separable such that we can write ∆ k = ∆f k with Hereby, V is the pairing strength and f k the gap symmetry function, which can be chosen according to the symmetry groups of the underlying lattice. We want to study the system in an out of equilibrium state excited by a laser pulse. As the net effect of an ultrashort THz pump pulse is that it changes the system abruptly, we model this effect in a general manner by considering an effective quantum quench, i.e. the pump pulse acts like a quench pulse. In literature [23], the usual way to go is an interaction quench, i.e. changing the interaction strength abruptly from its equilibrium value V to a new value V = gV , where g < 1. With this, the old equilibrium state is no longer an eigenstate of the quenched system. The disadvantage of this method is that an interaction quench cannot easily be realized experimentally in condensed matter. Only in cold atom systems, where the interaction strength can be tuned with Feshbach resonances [42], an implementation is feasible. Besides that, an interaction quench is not directly equivalent to the effect of a quench pulse and in addition, it always acts on the whole system isotropically, whereas a momentum-dependent excitation is more interesting to study, especially for unconventional superconductors with nontrivial pairing symmetry.
Considering all of this, we implement a momentumdependent state quench. Namely, the effect of a quench pulse is to deplete a small portion of the condensate and create a nonequilibrium quasiparticle distribution. What is happening in particular, is that the equilibrium distribution of the condensate at T = 0 with the quasiparticle energy E k = 2 k + (∆f k ) 2 , is modified as a result of the quench pulse, where the symmetry of the quenched state c −k↓ c k↑ does not necessarily have to be the same as in equilibrium. To implement such a quench in a controlled way, we introduce the quenched distribution which is an equilibrium distribution for a different symmetry f k with ∆ k = ∆f k and E k = 2 k + (∆f k ) 2 . This distribution is of course no longer the equilibrium distribution for the original system and the deviation can be controlled by the modified symmetry function f k = f k + gf q k , which is a small deviation of strength g from the original symmetry f k with the quench symmetry f q k . To solve the induced time-evolution after a quench, it is advantageous to perform a Bogoliubov transformation to new quasiparticle operators α k and β k with In the expression above, we choose the phase of the equilibrium order parameter ∆ to be zero. The Bogoliubov transformation is performed for this fixed value of the gap. The BCS Hamiltonian at arbitrary times for timedependent values of the order parameter ∆(t) in Eq. (1) then reads where we neglect constant terms and define In equilibrium, ∆ k (t) = ∆ k and the Hamiltonian becomes diagonal with R k = E k and C k = 0. We can express the anomalous expectation value c −k↓ c k↑ from the gap equation (3) with the new Bogoliubov quasiparticle operators and find

B. Iterated equation of motion method
We calculate the time-evolution of the Bogoliubov operators with the help of Heisenberg's equation of motion where the HamiltonianH(t), α † k (t) and β † k (t) are the operators in the Heisenberg picture. To evaluate these equations, we make use of the iterated equation of motion approach [43] to write the Bogoliubov operators in the Heisenberg picture as a time-dependent superposition of the equilibrium operators Inserted into Heisenberg's equation of motion, such an ansatz leads in general to an infinite hierarchy of equations for the prefactors a ik (t) and b ik (t), which have to be truncated at a certain order. However, for a bilinear Hamiltonian like the BCS Hamiltonian, the approach from Eq. (12) using only a second order ansatz for each operator is exact. This allows to derive a closed set of differential equations for the time-dependent prefactors. The derived equations are given in appendix A.
The coupled differential equations (A1) are solved numerically in a self-consistent manner by evaluating the time-dependent gap equation in each time step. All the used parameters are listed in Tab. I in the appendix. In order to increase the precision while keeping the numerical effort in a reasonable order, the 2d momentum space is discretized only in a small region around the Fermi level given by an energy cutoff c . This is justified by the fact that superconductivity only occurs near the Fermi momentum k F , while contributions far away are negligible. The momentum grid was chosen in polar coordinates, with N k points in radial and N ϕ points in azimuthal direction.

C. ARPES
To calculate the ARPES intensity, we use the often used three-step model and the sudden approximation [1], where the intensity is calculated via Hereby, I 0 (k, ω p ) is the transition probability of the initial state with momentum k to the final state after absorbing a photon with energy ω p . The Fermi distribution f (ω) ensures that only occupied states are considered. A(k, ω, t ) is the spectral function which determines the electron momentum and energy distribution and is the only quantity we calculate here. The spectral function can be computed from the lesser Green's function which, in the time-domain, is defined as [2] The required electron expectation value at different times can be computed from the time-evolution of the Bogoliubov quasiparticle operators The Fourier transform is calculated with respect to the time difference τ = t − t , i.e. the time between quench and probe pulse Following literature [2,36,38], in the expression for the Fourier transform above, we incorporated the finite width of the probe pulse by multiplying the Green's function with a Gaussian function  with full width at half maximum τ p prior to the Fourier transform.

III. RESULTS
In the following, we will investigate the dynamics of two quantities. On the one hand, we will consider the position of the maximum of the energy distribution curve (EDC) at k = k F , i.e.
This quantity, i.e. the position of the maximum with respect to the Fermi level reflects the energy gap. On the other hand, we will follow the dynamic of the amplitude of the spectral function at ω = 2∆ ∞ , where the strongest dynamic can be expected, i.e.
While the first quantity is typically used to trace the time-evolution of the system [36,37], it will appear that an investigation of the second quantity yields even more information.
A. Interaction quench for s-wave superconductor Before studying superconductors with nontrivial gap symmetry, we first consider the s-wave case with f k = 1, where we perform a simple interaction quench to trigger the time dynamic. In Fig. 1, the time-evolution of the spectral function can be found. One observes that after the quench ( Fig. 1(b)), spectral weight is shifted towards and above the Fermi level compared to the equilibrium case ( Fig. 1(a)). This can especially be seen by determining E as defined in (19). This position of the maximum, indicated by a circle in Fig. 1(a) and Fig. 1(b) is shifted from its initial value |E(t = −∞)| = ∆ to a smaller value |E(t = 0)| < ∆, which indicates a suppression of the energy gap after the quench. To trace the induced dynamic, we look at the time-evolution of the EDC at k = k F shown in Fig. 1(c) and Fig. 2, where we extract the two quantities E(t ) and A(t , k = k F ).
For a sufficient short probe pulse (small τ p ), the timeresolution is high enough to directly observe the induced Higgs oscillations of the order parameter as oscillations of E(t ), i.e. an oscillation of the spectral function in energy. The extracted curve E(t ) is shown in Fig. 1(d) in comparison to the calculated dynamic of the order parameter ∆(t ), where a close accordance can be observed.
The pulse width plays a crucial role for the resolution of the gap dynamic in E. If the probe pulse is too wide, it cannot resolve the oscillations of the system and E will be constant in time at the mean value ∆ ∞ , the value of the order parameter after long time. The influence of the probe pulse width is discussed in more detail in  appendix B.

Symmetry function
In addition to the oscillation of E(t ), we also look at the dynamic of the amplitude of the spectral function A(t , k = k F ) as defined in Eq. (20), shown as well in Fig. 1(d) for comparison. We find oscillations with the same frequency as the order parameter. This is confirmed in Fig. 1(e), where the Fourier transform of the order parameter ∆(t ), the maximum of the EDC E(t ) and the spectral function A(t , k = k F ) are compared. All three quantities show the same oscillation with a frequency of ω = 2∆ ∞ , the energy of the Higgs mode. Thus, the spectral function provides a direct measure of Higgs oscillations for superconductors in nonequilibrium, both in energy and amplitude.

B. State quench for d-wave superconductor
After having considered the s-wave case with an isotropic interaction quench, we will proceed to a d-wave superconductor, which we will quench in different symmetry channels with the help of a state quench as defined in Sec. II A. The symmetry function of a d-wave superconductor close to the Fermi level can be written as f k = cos(2ϕ), where ϕ is the polar angle. In order to excite the system in a systematic way, we choose quench symmetries which are fundamental with respect to the lattice point group. To this end, the superconductor is quenched in all different symmetry channels for the D 4h point group, which is the underlying lattice point group for the important d-wave example of cuprate superconductors. We define  the following four quenches: The symmetry function and the quenched symmetry function, as well as the spectral function after the quench, can be found in Fig. 3. First of all, in the equilibrium case, the spectral function closely resembles the absolute value of the symmetry function. In the antinodal directions at ϕ = [0, π/2, π, 3π/2], a full open gap can be observed with |E| = ∆. Moving ϕ towards the nodal directions ϕ = [π/4, 3π/4, 5π/4, 7π/4], the maximum E approaches zero.
The influence of the different quenches is also reflected in the change in the spectral function. The A 1g -and B 2gquenches shift the position of the nodes. In the case of the A 1g -quench, the shift is in opposite directions relative to the lobe maxima, which creates lobes with different sizes. The B 2g -quench shifts all nodes in the same direction, which results in a rotation of all lobes. This can be seen by a transfer of spectral weight above the Fermi level at exactly these points in the spectral function. The B 1g -quench does not change the symmetry at all and the A 2g -quench only shifts weight inside the lobes, such that for these two quenches, no obvious change in the symmetry of the spectral function is observable.
This fact is crucial for the induced Higgs oscillations of the gap, as for the A 1g -and B 2g -quenches, the shift of the nodal lines creates a second low-lying Higgs mode Amplitude (a.u.) dynamically, whereas for the other two quenches, it will not. This can be seen in Fig. 4, where the oscillations of the order parameter and their Fourier transform is shown for the different quench symmetries. While the B 1g -and A 2g -quenches only show a single broad peak around ω = 2∆ ∞ , the A 1g -and B 2g -quenches show a two peak or kink structure in the spectrum with one peak around ω = 2∆ ∞ and a second peak at lower frequencies.
Now, we will analyze the amplitude oscillations of the spectral function A(t ) in more detail to gain a deeper understanding of the dynamic creation of the low-lying Higgs mode. According to the gap equation Eq. (3), the value of ∆(t) is obtained by a summation over the condensate c −k↓ c k↑ (t). Both, the anomalous Green's function c −k↓ c k↑ (t) and the normal Green's function c † k↑ c k↑ (t) share a similar dynamic as their equations of motion are coupled. This can be seen in Eq. (D7) in appendix D. Therefore, a momentum-resolved analysis of the spectral function, i.e. A(t , k, ϕ), which is proportional to the normal Green's function, can reveal important information about the dynamic processes. In comparison, an angle-resolved evaluation of the maximum of the EDC curve, i.e. E(t , ϕ), as considered in [37], does not give further insight as it traces only the momentum-averaged quantity ∆(t ).
To this end, we extract A(t , k, ϕ) and its frequency ω A k not only at k = k F but for all points in momentum space and plot its distribution in Fig. 5(a) for a d-wave superconductor quenched in the A 1g -and B 1g -channels. The resulting frequencies are momentum-dependent yet do not depend on the quench. Namely, the frequency at each momentum point k is given by This result is derived in appendix D. The same result for s-wave symmetry and an interaction quench was also found in [23]. A comparison of the extracted frequencies and the analytic formula shows perfect agreement.
In the summation process of the gap equation, the frequency with the largest weight will dominate the oscillations. In the case of s-wave with f k = 1, this turns out to be at k = k F , where k = 0, which results in a single main frequency ω = 2∆ ∞ . However, as the other frequencies still contribute, the superposition of all different frequencies drives the oscillations out of phase in the long-time limit, which results in the well-known 1/ √ t decay [21]. In the case of d-wave, there is a much larger variety of frequencies due to the angular dependence of the gap ∆ k . Therefore, the summation process creates a much stronger decay [24].
The second low-lying frequency for certain quenches can be explained in the same picture. In Fig. 5(b), we show the extracted amplitudes for ω A k at each momentum point, which correspond to the actual weight for each frequency. While there is no difference in the frequency distribution between different quenches in Fig. 5(a), we can clearly see a difference in the amplitudes between the A 1g -quench, which shows a second peak in the Higgs oscillations, and the B 1g -quench, which shows no second peak. The change in symmetry due to the quench creates a strong asymmetry in the weights regarding the positive and negative lobe direction in the case of the A 1g -quench, while the amplitudes for the B 1g -quench are symmetric. This can result in the summation process in an additional enhancement of frequencies other than ω = 2∆ ∞ . In Fig. 5(c), we show the histogram of the frequencies from all momentum points. The raw histogram (red curve) looks exactly the same for both quenches, with a peak at ω = 2∆ ∞ already without additional weighting. Now, if we weight the histogram with the amplitude of the oscillation, the picture changes (green curve). In the case of the A 1g -quench, a two peak or two kink structure becomes visible, very similar to the Fourier transform of the gap oscillation (blue curve), whereas the symmetric weighting in the case of the B 1g -quench only shows a one peak structure. Hence, a careful analysis of the momentumresolved frequencies in the ARPES spectrum reveals the dynamic of the condensate and explains the underlying processes leading to the collective Higgs oscillations.

C. Analysis of oscillation phase
In the previous section we have seen how the oscillation of the spectral function can provide a deep insight into the creation process of the collective Higgs oscillation of the order parameter. Now, we will analyze the phase of these oscillations at different points in momentum space, to resolve the collective condensate dynamic. We have seen that an A 1g -quench for d-wave changes the sizes of the positive and negative lobes, i.e. the lobes of one sign increase in size while the other decrease. The resulting movement of the condensate will therefore be an oscillation, where the smaller lobe increases while the larger lobe decreases. Thus, the movement of the lobes is out of phase (see Fig. 6(a)).
The same quench applied to a nodal s-wave superconductor with f k = cos 2 (2ϕ), which has the same nodal structure, will change its lobes all equally due to the same sign of the gap. The resulting oscillations are therefore in phase (see Fig. 6(b)). For comparison, the quenched symmetry function and the resulting spectral functions analogous to the d-wave case in Fig. 3 can be found in the appendix in Fig. 8. Even though the oscillation of the lobes is in phase for this quench, the gap oscillation shows a two peak structure (Fig. 4(d)) due to a shift of nodal lines. Thus, a comparison of the Higgs oscillation spectrum does not allow to distinguish between d-wave and nodal s-wave in this case.
To see whether we can observe the different condensate dynamic for these two gap symmetries in the spectral function, we extract the amplitude oscillations A(t , k = k F , ϕ) at ϕ = 0 and ϕ = π/2, i.e. at points in momentum space rotated by ϕ = π/2, which should oscillate with opposite phase in the case of d-wave and with the same phase for the nodal s-wave. The result can be found in Fig. 6(c) and (d) and confirms the expected oscillatory behavior. Oscillations of the spectral function's amplitude provide therefore not only information about the amplitude of the gap but relative phase information from different points in momentum space can be extracted as well. This shows that the information obtainable from the time-dependent amplitude of the spectral function allows to clearly trace the collectively induced condensate movement.

IV. SUMMARY AND DISCUSSION
In this article we demonstrate how the collective Higgs oscillations of the superconducting order parameter can be traced back and be observed as oscillations of the spectral function measured in THz tr-ARPES experiments. There are two kinds of dynamics, namely momentumindependent oscillations of the EDC maximum E(t ), which directly correspond to oscillations of the energy gap ∆(t ), and momentum-dependent oscillations of the amplitude of the spectral function A(t , k, ϕ). While the EDC dynamic only reveals information about the momentumaveraged quantity ∆(t ), the amplitude oscillation of the  spectral function allows to perform a momentum-resolved analysis of the condensate dynamic. Experimental conditions might dictate which quantity can be observed more easily. The EDC maximum oscillations require short enough probe pulses to resolve the Higgs oscillations. Furthermore, they are restricted in their amplitude by the oscillation amplitude of the gap, which can be in the sub-percentage range depending on quench strength and gap symmetry. In contrast, the amplitude oscillations of the spectral function can be up to several percentages in size and thus, be much larger for certain conditions. However, they might be more difficult to resolve in general as oscillations in energy.
A tracking of the amplitude oscillation frequencies ω A k of the spectral function A(t , k, ϕ) at different momenta leads to a deeper understanding of the microscopic creation process of collective Higgs oscillations ∆(t ). While the frequencies can be described with a single formula independent of gap-and quench symmetry, the relative weight of each frequency contributing in the collective amplitude oscillation of the order parameter is determined by the gap-and quench-symmetry in a nontrivial way. If there are asymmetries in the weights, additional Higgs modes can appear, which are dynamically created due to quenching the condensate in a symmetry channel different from its equilibrium value. Furthermore, quenching the condensate in defined symmetry channels and observing the oscillation phase of the spectral function allows to resolve the exact condensate movement. This condensate dynamic depends on the gap and quench symmetry. In this paper, the quench symmetry for a certain gap was given such that the resulting condensate movement, like an in-phase or out-of-phase oscillation of lobes of a d-wave or nodal s-wave superconductor are clearly defined. Importantly, this explicit construction of a quench symmetry allows to distinguish between the two different gaps with same nodal structure but different phase.
Whether and how all considered quenches can be realized experimentally is an open question and beyond the scope of this paper; however, there are several promising possibilities. In-plane excitations can lead to symmetry breaking due to the small but non-vanishing photon momentum [28], which allows to quench the system asymmetrically with respect to its ground state symmetry. Such momentum transfer can be enhanced and controlled in more detail in transient grating setups [39,40] or fourwave mixing experiments [41]. Once tailored quenches are realized experimentally, a new spectroscopic tool will be available to gain phase-sensitive information about the gap symmetry of unconventional superconductors.

ACKNOWLEDGMENTS
We are grateful to S. Kaiser for fruitful discussions. We also thank N. Cheng and H. Krull for initial numerical studies. We further thank the Max-Planck-UBC-UTokyo Center for Quantum Materials for collaborations and financial support.

Appendix A: Equations of motion
Heisenberg's equation of motion (11) for the timedependent quasiparticle operators α † k (t) and β † k (t) can be evaluated by using the ansatz (12). From the resulting equations, a comparison of coefficients leads to a set of coupled differential equations for the prefactors The initial values for the prefactors can be derived from the initial values of the quasiparticle expectation values. The quasiparticle expectation values expressed within the iterated equation of motion ansatz read where we used the fact that for T = 0, all equilibrium quasiparticle expectation values are zero. If the system is in equilibrium for t = t = 0, we easily find i.e. α † k (0) = α † k and β † k (0) = β † k . The solution for a nonequilibrium distribution is straight forward as well. However, we found it advantageous for more numerical stability to first use the relation following from the anticommutation rules for the timedependent Heisenberg operators, to rewrite Eq. (A2) as such that the solution reads This prevents dividing by zero in the calculation of a 1k (0) and b 1k (0).

Appendix B: Influence of probe pulse width
The probe pulse width determines the energy resolution in time [44]. If the probe pulse is very long compared  to the intrinsic oscillations of the system, it can only detect the average value of the oscillations, which results in a constant in energy position of the spectral function.

Symmetry function
Only amplitude oscillations of the spectral function can be observed in this case. However, if the probe pulse is short in time, i.e. shorter than the intrinsic oscillation period of the system, it can scan the oscillations which leads to additional oscillations of the spectral function in energy.
In Fig. 7, the spectral function of an s-wave superconductor after an interaction quench is shown for different probe pulse widths. An infinite long probe pulse with τ p = ∞ leads to a very sharp spectral function in energy. However, only amplitude oscillations (red curve) can be observed in this case, while there are no oscillations in energy (black curve). For decreasing probe pulse width, the energy resolution of the spectral function decreases and the spectral function peak becomes broader. Simultaneously, the time-resolution in energy increases and oscillations of the EDC maximum can be observed, which follow the oscillations of the energy gap (blue curve) closely.
Appendix C: State quench for nodal s-wave superconductor In analogy to Sec. III B and Fig. 3, we quench a nodal s-wave superconductor in all symmetry channels of the D 4h lattice point group. The result can be found in Fig. 8. In the areas, where the quench shifts nodal lines, we can observe a corresponding shift in spectral weight. If one calculates the resulting gap oscillations shown in Fig. 4, one can observe a two peak or kink structure for all quenches. We can understand this by the fact that all quenches shift the nodal lines such that there is an asymmetry in the weights of the k-dependent frequencies.
In the case of the B 1g -quench, this shift is smaller compared to the other quenches, which results in a lower frequency of the second Higgs mode. In general, the energy of the dynamically created Higgs mode depends on the quench strength, i.e. the deviation from the equilibrium state. , (D10a) . An inverse Laplace transform yields where ω k = 2 2 k + ∆(0) 2 f 2 k .
In the linear regime, we can approximate ∆(0) ≈ ∆ ∞ , such that we arrive at the formula Eq. (22) from the main text. We can see that the main oscillation frequency of the summand in the gap equation and both the anomalous and normal Green's function, i.e. the x-and z-component of the pseudospin in Eq. (D8) is given by The frequency is independent of the quench, however the individual weighting in momentum space, i.e. the amplitudes of each frequencies depend heavily on the quench. Namely, the prefactors, controlled by the gap symmetry f k and quenched symmetry f k , determine the dominating frequency in the summation process in a nontrivial way.