Phase signatures in third-harmonic response of Higgs and coexisting modes in superconductors

Third-harmonic generation (THG) experiments on superconductors can be used to investigate collective excitations like the amplitude mode of the order parameter known as Higgs mode. These modes are visible due to resonances in the THG signal if the driving frequency matches the energy of the mode. In real materials multiple modes can exist giving rise to additional THG contributions, such that it is difficult to unambiguously interpret the results. In this paper, we additionally analyze the phase of the THG signal, which contains microscopic details beyond classical resonances as well as signatures of couplings between modes which are difficult to observe in the amplitude alone. We investigate how the Higgs mode, impurities or Coulomb interaction affects the phase response and consider exemplary two systems with additional modes. We argue that extracting this phase information could be valuable in future experiments.


I. INTRODUCTION
Recent progress in ultrafast THz laser technology lead to an increasing interest in studying collective excitations of superconducting systems. Especially the investigation of the amplitude (Higgs) mode of superconductors lead to a new emerging field termed Higgs spectroscopy, where the Higgs mode is used as a probe for intrinsic properties of the system [1][2][3][4]. Intrinsically, a superconductor possesses two collective modes due to the spontaneous U (1) symmetry breaking: An amplitude oscillation of the order parameters, known as the Higgs mode and a phase oscillation, known as the Goldstone mode [5,6]. The Goldstone mode is a massless mode in the long-wave limit for uncharged systems. However, the Anderson-Higgs mechanism in a charged superconductor shifts its energy to the plasma energy. The Higgs mode is a massive mode with an energy at the quasiparticle energy 2∆ and thus is a low-energy excitation in the range of THz frequency.
Experiments to excite the Higgs mode are usually performed in either of two ways. One option is to quench the system with an ultrafast, single-cycle THz pump pulse to abruptly change the system's parameter and bring it out of equilibrium [7,8]. The order parameter starts to automatically oscillate around its new equilibrium state with the Higgs mode frequency. This oscillation is experimentally measured in a pump-probe geometry, where the probe pulse scans the dynamics of the system with a variable time-delay after the pump pulse [1].
The second option is to drive the system periodically with a multi-cycle THz pulse at frequency Ω. This enforces the order parameter to oscillate with twice the driving frequency 2Ω due to the quadratic excitation process [9][10][11]. Furthermore, this leads to a third-harmonic generation (THG) process, which can be measured in the transmitted electric field. Tuning the driving frequency into resonance with the Higgs mode energy, i.e. 2Ω = 2∆, a resonance peak is visible in the signal. This can be achieved either by varying the driving frequency or, as it is currently done experimentally, by changing the value of the order parameter ∆(T ) by sweeping the temperature T . The resonance can be used as a signature for the collective Higgs mode as it was demonstrated for the s-wave superconductor NbN [1,12].
Our work is motivated by a recent THG experiment on several cuprates [3]. The experiment revealed an interesting phase signature containing antiresonance behavior which cannot be explained by the excitation of a single collective mode. In our work, we therefore take into account the existence of another mode and we investigate the THG signal for such systems, where we concentrate on the phase of the 3Ω oscillations, which was not discussed theoretically so far. As it is well known, a driven oscillator shows an abrupt phase change if the driving frequency is varied across the eigenfrequency of the system. Furthermore, multiple coupled oscillators show an antiresonance, resulting from the interplay of driving force and coupling, where a minimum in the driving amplitude occurs accompanied by a negative phase change. This behavior occurs in many physical systems and we show that it is also visible in the THG response following from a microscopic calculation with coupled modes. Yet, the response is more complex compared to a classical model as modified mode propagators and susceptibility terms occur.
After a general analysis, we provide two detailed microscopic calculations of a two mode scenario. First, a coupling of the Higgs mode to a charge density wave and second a Higgs mode with a coexisting Bardasis-Schrieffer mode. With this paper, we propose that analyzing the phase of the THG signal in addition to the amplitude yields additional information valuable for understanding the interplay of superconductivity and other modes.
The paper is organized as follows. In Sec. II, we generally discuss the phase response of a single driven oscillator. We start from a classical model, proceed to a phenomenological Ginzburg-Landau theory for superconductivity and finally show a microscopic calculation in an effective action formalism. In Sec. III, we extend the single mode analysis to a second mode. We start again with a classical model and then discuss the general features of a microscopic theory. In Sec. IV, we explicitly calculate the THG response of a superconducting system with coexisting CDW. In Sec. V, we explicitly calculate the THG response of a superconducting system with Bardasis-Schrieffer modes. Finally, we summarize and conclude in Sec. VI

II. PHASE SIGNATURE OF A SINGLE MODE
Before studying the full microscopic quantum mechanical model for superconductors and its collective modes, let us first consider a simple classical system. This will allow us to define and observe the crucial features which are important for the later discussion. Hereby, we investigate classical driven oscillators which represent the collective modes of the system.

A. Harmonic oscillator
It is well known that a driven harmonic oscillator has a characteristic amplitude and phase response which depends on the driving frequency. With the eigenfrequency ω 0 , damping factor γ, driving amplitude F 0 and driving frequency Ω, the equation of motion for the displacement x(t) readsẍ (t) + ω 2 0 x(t) + γẋ(t) = F 0 cos(Ωt) .
The steady-state solution can be written as x(t) = A cos(Ωt − φ), where the frequency-dependent amplitude A and phase φ are given by One observes that the amplitude has a resonance peak at Ω = ω 0 which is accompanied by an abrupt phase change from 0 to π. Thus, the oscillation is in-phase with the driving frequency below the resonance and lags behind with opposite phase above the resonance. The amplitude and phase is plotted in Fig. 1 for different damping values γ. While for small damping a pronounced resonance peak is visible in the amplitude, for large damping, the resonance peak is heavily suppressed and broadened. In contrast, the phase still shows a phase change from 0 to π, even though it is broadened as well. This means that both amplitude and phase have a signature of the resonance, yet the phase change signature is more robust against the influence of damping. Hence, in a strongly damped system with suppressed resonance peak, the eigenmode would still be identifiable via the phase signature.

B. Ginzburg-Landau model
Let us investigate now whether we can observe such a behavior for THz-driven collective modes in superconductors as well. Hereby, the oscillator corresponds to a collective mode which is driven by a THz light field. In the experiment, the driven collective mode is not measured directly. Instead the induced current proportional to the transmitted electric field is recorded.
As a first step, we investigate the phenomenological Ginzburg-Landau model, where we will consider amplitude and phase fluctuations. The time-dependent Lagrangian of a superconductor coupled to a gauge field is given by where ψ is the superconducting order parameter, D µ = ∂ µ + ieA µ the covariant derivative with the four-vectors ∂ µ = (∂ t , −∇) and A µ = (Φ, −A) and electromagnetic field tensor F µν = ∂ µ A ν − ∂ ν A µ in units where c = 1. In principle, the Lagrangian could also contain additional linear derivative terms. Yet, we assume perfect particle-hole symmetry, such that the time dynamics of a superconductor is described only by a second-order derivative term [5,6]. The potential V (ψ) = α|ψ| 2 + β 2 |ψ| 4 is the free energy of a superconductor with β > 0 and α = α 0 (T − T c ) such that for T < T c the potential takes the form of a Mexican hat with the ground state ψ 0 = −α/β. We introduce amplitude (Higgs) fluctuations H(r, t) and phase (Goldstone) fluctuations θ(r, t) via ψ(r, t) = (ψ 0 + H(r, t))e iθ(r,t) , and choose a gauge A µ → A µ + 1 e ∂ µ θ and ψ → ψe −iθ . Then, the Lagrangian up to second order in the fluctuations reads Hereby, the phase fluctuations are removed from the Lagrangian by the chosen gauge and are implicitly included in the longitudinal component of the transformed gauge field A µ which obtains an additional mass term ∝ A µ A µ . This effect is known as the Anderson-Higgs mechanism [19]. Calculating the equations of motion for the Higgs mode H, neglecting spatial fluctuations for q → 0 and choosing a gauge with Φ = 0, yields The dynamics of the Higgs oscillations is governed by a harmonic oscillator with frequency ω 0 = √ −2α. The driving term is quadratic in the vector potential A(t). With a periodic light field A(t) = A 0 cos(Ωt), the system is effectively driven by 2Ω such that the resonance in the system occurs at 2Ω = ω 0 . Thus, on a phenomenological level, the collective Higgs oscillations of a superconductor and its amplitude and phase signature is exactly described by the classical model discussed before. The measured transmitted field is described by the induced current given by [9] A nonlinear third-harmonic component in the current is induced as A(t) · H(t) ∝ cos(3Ωt − φ) + . . . The resonance behavior of the amplitude and phase in the current j(t) is directly given by the Higgs response H(t).

C. Microscopic BCS model
While in the phenomenological model the coupling of light to the system contains no further details, in a microscopic model additional effects with frequency-dependent susceptibilities occur. Furthermore, there are quasiparticles in the microscopic model which render the Higgs mode less stable due to the additional decay channel.
To address these effects, we proceed to the full microscopic theory using an effective action approach [10,20].

The BCS Hamiltonian reads
Hereby, k = ξ k − F is the electron dispersion ξ k measured relative to the Fermi level F and c † k,σ or c k,σ the electron creation or annihilation operators. The separable BCS pairing interaction is given by V k,k = V f k f k with pairing strength V and symmetry f k . A coupling to light represented by the vector potential A(t) is realized by minimal coupling k → k−A(t) . An expansion in powers of A(t) yields the lowest non-vanishing diamagnetic coupling term shown above, while the linear paramagnetic coupling ∝ ∂ i A i (t) vanishes due to parity symmetry. In the expression, we have introduced the short-hand notation ∂ 2 ij = ∂ 2 kikj . Here, we initially neglect long-ranged Coulomb interaction and the coupling to phase fluctuations which is important in real materials. We will show later in Sec. II D that including Coulomb interaction does not affect the phase signature. The action of the system in imaginary time τ is given by We perform a Hubbard-Stratonovich transformation introducing the bosonic field ∆, with amplitude fluctuations ∆(t) = ∆ + δ∆(t). After integration of fermions, we split the action in a mean-field and fluctuating part, which we expand up to fourth order in A. For more details about the calculation see Appendix A. The effective action with Matsubara frequencies iω m in fourth-order of the vector potential reads Hereby, H −1 (iω m ) is the inverse Higgs propagator defined as the renormalized pairing interaction V The susceptibilities are given by with and the BCS Green's function G −1 0 = iω m τ 0 − k τ 3 +∆ k τ 1 where τ i are Pauli matrices. The indices ∆ and A 2 in the susceptibilities represent the vertices, i.e. the coupling to the Higgs propagator via f k τ 1 or the coupling to light via ∂ 2 ij k τ 3 , respectively. Integrating out the amplitude fluctuations and after analytic continuation iω m → ω+i0 + one obtains There are two contributions in the action, one containing the Higgs oscillations and one the quasiparticle response. These contributions are shown diagrammatically in Fig. 2(a) and (b). For simplicity, we will only consider linear-polarized light in x-direction, such that we can neglect the polarization indices in the following. With this, the thirdharmonic response is given by (15) with the Higgs (H) and quasiparticle (Q) contribution Comparing the response j (3) with the phenomenological Ginzburg-Landau model in Eq. (7), we can observe several differences which modify the response. First, the Higgs propagator H(ω) is a more complex object and does not have a simple resonance pole as we will see. Second, light does not directly couple to the Higgs mode but through the susceptibility χ ∆A 2 (ω). Third, there is an additional quasiparticle response given by χ A 2 A 2 (ω). In the following, let us disentangle these effects. Evaluating the Matsubara sum and rewriting the momentum sum as integral assuming s-wave symmetry, the Higgs propagator can be analytically evaluated at T = 0. Concentrating on the pole structure one obtains It does not have a simple pole but a square root term in the denominator. Transformed into time-domain, this leads to a power-law decay of the Higgs mode. It can be understood as a decay into quasiparticles as the Higgs mode energy overlaps with the quasiparticle continuum at 2∆. In addition to the obvious consequence of stronger damping, it also affects the phase response. The square root reduces the π phase change at the resonance frequency to π/2. Thus, the driven amplitude oscillation only lags behind a quarter cycle at high frequencies instead of being completely anti-phase as found in the phenomenological model. Next, we check how the electron bubbles χ ∆A 2 (ω) generating the light-Higgs coupling, affects the phase response. In the expression in Eq. (16), the term occurs twice evaluated at opposite frequency. It can be written as its absolute squared value and thus is not affecting the phase. Finally, let us examine the quasiparticle response which is actually known to be much larger than the Higgs response [10]. Evaluating the Matsubara sum of the respective susceptibility and solving the momentum sum (see Appendix A) one obtains for the pole structure The quasiparticle response has the same square root pole structure as the Higgs mode, leading to the same π/2 phase change at the resonance frequency. In Fig. 3(a) and (b) the amplitude and phase of the diamagnetic Higgs and quasiparticle response is shown using k = −2t(cos k x + cos k y ) − µ, t = 10 meV, µ = −10 meV, ∆ = 1 meV and a residual broadening ω → ω + i0.05 meV. Hereby, the momentum sums are evaluated numerically on ) Ω (∆) Confirming the analytic study, we can see that phase shows a π/2 phase change at the resonance frequency 2Ω = 2∆. Above the resonance, a drift is observable to higher values for the quasiparticles and lower values for the Higgs mode.
As it has been emphasized in literature [10], the Higgs mode is much smaller in the clean-limit BCS theory.

D. Influence of Coulomb interaction
As a next step, we discuss the influence of Coulomb interaction given by an additional term in the Hamiltonian where V (q) is the Coulomb potential. We follow [10] and decouple the Coulomb interaction by means of an additional Hubbard-Stratonovich transformation introducing the density field ρ(q, τ ) = ρ 0 + δρ(q, τ ) and allow amplitude and phase fluctuations in the superconducting order parameter ∆(τ ) = (∆ + H(τ ))e iθ(τ ) . With this, one obtains for the fourth order action The susceptibilities are given in Appendix B. Integrating the fluctuations and using 1/V (q) → 0 for q → 0 one obtains for the third-order current The Coulomb interaction renormalizes the Higgs and quasiparticle response. Yet, due to obtained structure, the phase signature is not changed as the expressions in the nominator do not contribute due to the absolute square and the χ ρρ term has the same phase behavior as the unrenormalized propagators. This can be seen in Fig. 3(c),(d), where the respective expressions are numerically evaluated with the same parameters of the previous section. Except global scaling factors and small deviations resulting from the 1/χ ρρ contribution the result is basically unchanged with a phase change of π/2 at the resonance.

E. Influence of impurities
Recently it was pointed out by several papers [14,15,[21][22][23][24] that nonmagnetic impurities allow an additional paramagnetic coupling of light to the condensate. This is shown diagrammatically in Fig. 2(c) and (d) where the light-coupling vertices are in the τ 0 -channel. While these diagrams vanish in the clean limit, they have been shown to dominate the optical nonlinear response even for small disorder. Here, we adopt the Mattis-Bardeen approximation first applied to the nonlinear response in [14] and subsequently formulated in the effective action framework in [15].
Following [15], we implement a 3D continuum model, where we can express the THG current within the Mattis-Bardeen approximation as with the Higgs (H) and quasiparticle (Q) susceptibilities The triangle and square bubbles are defined as and are shown in Fig. 2(c) and (d). The transition matrix element J kk = k| ep m |k is approximated by a Lorentzian distribution with impurity scattering rate γ, Fermi velocity v F , and density of states at Fermi surface N (0). We choose the parameters ∆ = 2 meV, mass m = 0.78m e of the parabolic band dispersion, F = 1 eV, and impurity scattering rate γ/∆ = 10. We evaluate Matsubara sums analytically and numerically compute the momentum integrals. For further details about the calculation see [15]. While the Mattis-Bardeen approximation may not quield qualitatively accurate results in the nonlinear response, it serves well to discuss qualitative differences of the phase response compared to the clean limit. The resulting amplitude and phase of the dirty superconductor are shown in Fig. 3(e) and (f). We find a pronounced resonance peak at 2Ω = 2∆. Here, the Higgs contribution is no longer subdominant but instead gives the main contribution to the THG signal. The resonance peak is accompanied by a positive phase jump of roughly π across the resonance. The detailed structure of this phase response as well as the value of the phase jump show some weak dependence on material parameters.
The more complex phase structure in the dirty limit can be understand as follows: While the clean phase response is only given by the Higgs propagator due to the cancellation of the phase of electronic susceptibilities χ ∆,A 2 in Eq. (18), the phase of the susceptibilities in the dirty case no longer cancels but gives an additional additive contribution to the phase. It is represented by the fermionic triangles χ AA∆ (2Ω, −Ω)χ AA∆ (−2Ω, −Ω) and shown in Fig. 2(c). Thus, the phase response in the dirty limit is not only given by the Higgs propagator but has an additional contribution from the electron-mediated microscopic coupling of light to the Higgs mode.

III. PHASE RESPONSE OF TWO MODES
Now, we will consider systems which contain two modes and study the interaction between these. Again, we start by an analysis of the classic analogon of two coupled oscillators to understand the fundamental properties before proceeding to a microscopic model.

A. Coupled oscillators
If there are two modes in a system, interference effects occur in the driven system which can lead to the so-called antiresonance phenomenon. The usual way to understand this effect is based on the assumption that there are two modes in the system which are coupled and only one of these modes is externally driven. For a particular driving frequency, the external force on the driven mode cancels exactly with the force induced by the other coupled mode such that the amplitude of the oscillation of this mode vanishes -thus the name antiresonance. Furthermore, the antiresonance is accompanied by a negative phase jump of π, therefore it goes in the opposite direction compared to a resonance.
The same phase signature can also be obtained when both oscillators are driven and the observed signal is comprised of the sum of both oscillation amplitudes. Here, this effect is a trivial consequence of a destructive interference and does not necessarily rely on a coupling between the modes. An additional coupling between the modes allows for a tuning of the antiresonance frequency. We refer to this scenario as antiresonance behavior as well.
To make this effect more clear, let us first investigate again the classic model where we consider now two coupled and driven oscillators described by the following equations of motion There are two oscillators x 1 (t) and x 2 (t) with individual eigenfrequencies ω i , dampings γ i and driving amplitudes F i but same driving frequencies Ω. The coupling between the modes is controlled by the constant g. Using the complex variable method ansatz wherex i (t) =Â 1 (Ω)e iΩt withÂ i (Ω) = A i (Ω)e −iφi(Ω) , we write the equations in matrix form where we define the "propagator" of the oscillators as We also consider the total response One obtains for the complex amplitudeŝ In Fig. 4 we show a numerical evaluation of the individual and total amplitudes A i and phases φ i for two distinct cases (see Appendix C for the exact expressions). In the first column, the two oscillators are coupled, i.e. g = 0, but only the first oscillator is driven F 2 = 0. In the second column, the two oscillators are uncoupled, i.e. g = 0, but both oscillators are driven F i = 0. The first scenario (left column) corresponds to the usual definition of the antiresonance, namely a destructive interference between the driving force and the force from the second oscillator due to the coupling. The dip between the two resonance peaks and the negative π phase change, is clearly visible for the first oscillator (red curve). The energy of the antiresonance ω A is determined byP 1 = 0, which leads to ω A = ω 2 , i.e. the antiresonance occurs at the energy of the other undriven mode. The total response A T and φ T (blue curve) also shows this behavior resulting from the antiresonance of the first oscillator. Yet, the energy of the antiresonance is shifted as a result of the second superposition scenario.
We can further see that the finite coupling shifts the resonances frequency with respect to the uncoupled eigenfrequencies ω i . The resonance frequency for the lower modes is decreases, while the resonance frequency of the higher mode is increased. The energies are given by the poles of the renormalized propagators For the shown parameters, this results in a resonance peak below ω 1 = 1 and a resonance peak above ω 2 = 2. The total response (blue curve) of the second scenario (right column) shows a very similar behavior, namely a dip in between the two resonance peaks and a negative π phase change. However, in this case the negative phase jump does not result from an individual oscillator, both individual oscillators (red and green curve) do not show this behavior. It rather results from the superposition of the two oscillations where the sum of both cancel out in an intermediate position between the resonances. This energy is determined by P 1 + P 2 = 0 leading to As the sum of both undergoes a sign change from negative to positive a negative π phase change occurs naturally. The resonance frequencies in this scenario are not changed and still occur at ω i .
To summarize, the antiresonance behavior of the total response of two oscillators can have different origins. It can be either controlled by the coupling between the oscillators or the interference if both modes are driven.

B. Microscopic theory
Let us now investigate whether this behavior is observable in a microscopic model as well. For now, we will make some general arguments assuming that there are two modes 1 and 2 in the system, for example the Higgs mode and a second collective mode. In Sec. IV and Sec. V we will consider specific examples.
Taking into account the general form of the effective action and the analysis of the classical oscillator system, we are anticipating the results of the next sections and postulate the general structure of the response. The fourth order effective action for two modes reads where Hereby, P i stands for the propagator of mode i, χ i,j for different coupling susceptibilities and A the vector potential, where the polarization indices are not included. The tilde denotes a renormalization due to the other mode. This can be understood as an RPA series renormalization of the propagators shown in Fig. 5(d) and expressed as which leads toP The fourth-order kernel K (4) explicitly reads with These terms are diagrammatically shown in Fig. 5 After these general remarks, let us now consider specific examples of two microscopically coupled modes in the next sections.

IV. HIGGS AND CHARGE DENSITY WAVE
As a first example of two coupled collective modes, we will consider a coexisting superconducting and charge density wave (CDW) system. The amplitude modes are schematically shown in Fig. 6(a) in the picture of the free energy. An example for such a scenario is NbSe 2 , where the coupling of the Higgs mode to a CDW phonon was observed in Raman response [25,26] and theoretically investigated by several authors [18,27,28]. Another relevant system are cuprates, where superconductivity and fluctuating charge order has been reported in the underdoped regime [29,30]. This scenario might be a possible explanation of the antiresonance behavior observed in recent THG experiments [3].
To model the system, we follow [18] and start from the BCS Hamiltonian in Eq. (8)    of momentum Q responsible for creating the charge order and a coupling to electrons with strength g. The Hamiltonian is given by with Hereby, b † q or b q are the phonon creation or annihilation operators and ω Q the energy of the CDW phonon. The electron phonon coupling is controlled by g · g k with strength g and momentum dependence g k .
To simplify the calculation, we will make the following assumptions. We consider a 2d square lattice with a tight-binding dispersion and nearest-neighbor hopping t at half-filling, namely k = −2t(cos k x + cos k y ). As it was shown in [18], a finite chemical potential leads to a qualitative similar result. Choosing Q = π, π we have perfect nesting and a commensurate CDW with k+2Q=k where k+2Q = k and k+Q = − k . We assume an swave superconductor with f k = 1 and an anisotropic s-wave CDW with g k = | cos k x − cos k y |.
We start from the action of the system, where we introduce a CDW order parameter D k = Dg k with D = g b Q + b † −Q and the superconducting order parameter ∆ using a Hubbard-Stratonovich transformation. Details of the calculation can be found in Appendix D. Please note that we neglect here the Coulomb interaction and phase fluctuations as we have shown in Sec. II D that they do not affect the phase signature. Furthermore, in the half-filled case, as considered here, its influence vanishes completely as the system has perfect particle-hole symmetry [18].
After integration of the fermions and expansion of the action at Gaussian level for amplitude fluctuations ∆(t) = ∆ + δ∆(t) and D(t) = D + δD(t), the effective fourth order action reads where H is the Higgs propagator and P the renormalized Phonon propagator. It is defined as with the bare phonon propagator P 0 = −2ω Q /(ω 2 Q − (iω m ) 2 ) and χ DD the susceptibility describing the influence of the CDW. The susceptibility on the off-diagonal in M leads to a coupling between Higgs and CDW and the expression is equivalent to the coupled oscillator system Eq. (30) in the previous section. An integration of the amplitude fluctuations finally leads to with where we identify the renormalized Higgs and phonon propagator as shown diagrammatically in Fig. 5(d) The expression of the action follows exactly the general structure as shown in the previous section and the resulting diagrams are the ones shown in Fig. 5(a)-(c).
Let us analyze the result in more detail. First, we evaluate the expression for the bare Higgs propagator. One finds For D = 0, we would retain expression Eq. (17). However, for finite CDW gap D, the Higgs mode energy 2∆ no longer coincidences with the quasiparticle excitation gap ∆+D and, as already pointed out by [18], the Higgs mode becomes stable as its energy is below the gap. This has consequences for the pole structure, as the Higgs mode now has a simple pole without square root such that we can expect a π phase shift when varying the driving frequency ω from below to above the Higgs mode energy. The CDW phonon propagator can be evaluated and reads with The phonon propagator has a simple pole leading to a phase change of π.
Assuming linear polarized light in x-direction as in Sec. (II C), we can write the action as where the kernel is given by K (4) = χ H + χ P + χ M with the Higgs (H), phonon (P) and mixed (M) contributions It has the same structure as Eq. (40c) introduced as the general response for coupled modes. With the insight of the previous sections, we can expect an antiresonance behavior with a negative phase change of π in between the two resonances where a phase change of positive π should occur. To confirm, we calculate numerically the total THG response as function of frequency and temperature. The temperature dependence is necessary to compare with experimental results where only the temperature can be varied for fixed driving frequency. The calculation for a set of parameters ∆ = 2.5 meV, ω 0 = 15 meV, D = 3 meV, t = 10 meV, ω → ω + i0.1 meV is evaluated on a 2d grid with 2000 × 2000 points and shown in Fig. 7. Hereby, the CDW phonon energy is above the Higgs mode energy. The top row shows the THG intensity, and the bottom row the THG phase. The first column is a 2d plot of the THG signal as function of frequency and temperature. Thus, vertical cuts in this plot, shown in the second column, correspond to varying the frequency for fixed temperature and horizontal cuts, shown in the third column, correspond to varying the temperature for fixed frequency.
The result fulfills our expectation of the previous general analysis. Looking at the 2d plot in the first column, we can see the Higgs mode as a sharp resonance peak which follows the temperature dependence of the gap. However, due to the coupling to the CDW, the energy of the Higgs mode is renormalized and shifted to lower frequencies. This resonance peak is accompanied by a positive phase jump of π as the Higgs mode is a stable mode below the total gap as discussed before. At a slightly higher energy, we observe an antiresonance behavior with a dip in the amplitude and a negative phase jump of π. At higher energy, we observe a second resonance peak at the renormalized CDW phonon frequency with an associated positive phase jump of π. This resonance peak follows the temperature dependence of the CDW gap. Please note, due to the residual broadening, the positive phase change at the Higgs mode and the negative change at the antiresonance is slightly lower than π.
The temperature dependence of both modes are similar and follow roughly a quarter circle as can be seen in the left column in Fig. 7. Thus, vertical cuts along the frequency and horizontal cuts along the temperature reveal, in principle, the same information. Resonance peaks and phase changes are visible in both signals. Yet, obtaining single cuts at unfavorable positions, e.g. in between the modes, or limited variation range of parameters might miss several features. Experiments with a large variation of temperature and frequency are therefore necessary to reveal the full information.

V. HIGGS AND BARDASIS-SCHRIEFFER MODE
As a second example, we consider a superconducting system where the ground state is dominated by one symmetry, yet fluctuations in a subleading pairing channel are allowed. These fluctuations are known as Bardasis-Schriefer modes [31,32] and might exist for example in iron-based superconductors [33][34][35], where multiple pairing instabilities occur on different bands. While most studies investigated the signature of a Bardasis-Schrieffer mode in Raman spectra, recent work has shown that such modes can also be excited with THz light [17]. This mode is schematically shown in Fig. 6(b) and will be discussed in the following.
Here, we will consider an s-wave ground state and allow fluctuations in the d-wave channel. We use the BCS Hamiltonian in Eq. (8) with a sum of separable pairing interactions with the (anisotropic) s-wave symmetry f s k = (cos k x + cos k y )/2 and d-wave symmetry f d k = (cos k x − cos k y )/2. Choosing k = −2t(cos k x + cos k y ) − µ with t = 6 meV, we solve the gap equations with ∆ k = i ∆ i f i k for i = s, d for varying Fermi level µ and symmetry ratio V d /V s . This phase diagram is shown in the Appendix in Fig. 10. In the following, we choose parameters µ = −12 meV, V d = V s and ∆ = 2 meV, where the ground state is s-wave but still close to the d-wave transition. The residual broadening is ω → ω + i0.05 meV. Please note, for the chosen parameters, i.e. the anisotropic s-wave and the energy dispersion, the maximum of the gap at the Fermi level is ∆ max ≈ 1 meV, such that the Higgs mode energy is at ω H ≈ 2∆ max and not at ω H = 2∆.
If the system is excited, we allow fluctuations in both symmetry channels (56b) Hereby, δ∆ s (t) corresponds to the usual Higgs mode of the dominant symmetry channel, while δ∆ d (t) are amplitude fluctuations of the subleading channel orthogonal, i.e. in the imaginary axis direction. This is the Bardasis-Schrieffer mode. As shown in [32], subleading fluctuations in the parallel or real channel do not lead to a finite energy mode. After integrating out the fermions (see Appendix E), we obtain the same structure of the effective action as Eq. (43) with where H(iω m ) is the usual Higgs propagator and B(iω m ) the Bardasis-Schrieffer propagator defined as The susceptibilities are defined as with X αβ defined in Eq. (13). In analogy to the previous sections, the fluctuations are integrated out which leads to with and the renormalized propagators For monochromatic, linear polarized light with polarization angle θ, the THG current parallel to the light polarization is given by (see Appendix E) j ) Ω (meV) Ω (meV)

mixed (M) and quasiparticle (Q) contribution
The response has again the same structure of coupled modes as discussed before.
To get a first insight into the Bardasis-Schrieffer mode, we evaluate the expression for the propagator analytically for T = 0 and assuming a constant density of states at the Fermi level. One obtains for continuum implementation (see Appendix E) Using these simplification, one can see that for V d → V s the pole of the propagator shifts to zero, while for V d → 0 the pole approaches 2∆. Furthermore, due to the √ 4∆ 2 − ω 2 term in the numerator, the expression is always zero at ω = 2∆, which leads to a negative phase change of π/2. Thus, we expect a positive phase change of π at the Bardasis-Schrieffer mode energy below 2∆ and a negative phase change of π/2 at 2∆.
As the coupling term χ ∆B contains the product of the two symmetry functions f s k f d k , which are orthogonal, the term vanishes and the Higgs and Bardasis-Schrieffer modes are not coupled but contribute individually to the response. However, for strong pulses beyond the Gaussian level, as used in pump-probe experiments, a finite coupling between Higgs and Bardasis-Schrieffer mode can exist [17]. In addition, due to the symmetry function f d k in χ ij BA 2 , the coupling strength of light to the Bardasis-Schrieffer mode will depend on the polarization. For our band structure, θ = π/4 will correspond to the A 1g symmetry, such that we expend a vanishing of the d-wave (B 1g ) Bardasis-Schrieffer mode. This polarization sensitivity has also been discussed in [17].
Next, we evaluate the THG response at T = 0 numerically for two different polarization angles θ = 0, π/4 and show the individual contributions in Fig. 8. As we have anticipated, there is no coupling between Higgs and Bardasis-Schrieffer mode and the mixed contribution is zero (not shown). Thus, the quasiparticle (blue curve) and Higgs (red curve) response is the same as in the pure system discussed in Sec. II C. The Higgs contribution is polarization independent, while the quasiparticle contribution gets reduced for θ = π/4. The Bardasis-Schrieffer (green curve) mode has a strong polarization dependence. For θ = 0, the expected phase behavior originating from the Bardasis-Schrieffer propagator B(ω) is visible. At the resonance energy a positive π phase change occurs and at the energy of the Higgs mode a negative π/2 phase change occurs. The intensity shows a strong peak at the Bardasis-Schrieffer resonance. The small peak at the Higgs mode energy is a result of the susceptibility χ BA 2 . For θ = π/4, the Bardasis-Schrieffer mode is not excited.
Having no coupling between the modes, the resulting phase signature is influenced only by the interference of the pure, individual contributions. The total THG signal (yellow curve) consists of two resonance peaks at the original, unrenormalized frequencies each accompanied by a positive phase change. In between we see an antiresonance behavior with a dip in the intensity and a negative phase change.
In analogy to the previous section, we also calculate the frequency and temperature dependence of the total THG signal. The result is shown in Fig. 9. A before, the temperature dependence follows roughly a quarter circle (left column). Thus, both vertical and horizontal cuts along the frequency or temperature axis show a similar result and the resonance and antiresonance behavior is visible in both cases.

VI. CONCLUSION
In this paper, we investigated the phase signatures of the third-harmonic (TH) signal generated by driving superconductors with THz light. Hereby, the phase of the oscillatory TH signal is measured with respect to the phase of the first-harmonic. While it is well known  that resonances in the intensity of the signal occur if the driving frequency matches with the energy of intrinsic modes, in this paper we showed that the phase change is a robust feature as well and additionally encodes intrinsic properties of the system.
From a classical point of view, the resonance peak of a driven oscillator is accompanied by a positive phase change of π. This phase signature is more robust against damping than the peak itself, as it is still observable even if the peak is heavily suppressed. Coupled modes lead to antiresonance behavior where the interference of driving force and coupling leads to a dip in the oscillation amplitude and a negative phase change of π. Usually, the antiresonance is understood as a feature of a single oscillator which is externally driven and coupled to an undriven second oscillator. Yet, for oscillators which are both driven but uncoupled a similar feature also occurs in the superposition of the oscillation amplitudes as a destructive interference of two oscillations.
In a microscopic BCS model for superconductors, the amplitude (Higgs) mode can be driven by a nonlinear, quadratic process such that the effective driving frequency is doubled which leads to an induced third-harmonic current. This driving of the Higgs mode is similar to a classical oscillator, yet, the microscopic details lead to a modification of the phase signature.
First of all, as the energy of the Higgs mode coincidences with the quasiparticle excitation energy at the gap 2∆ it is intrinsically damped. This leads to a propagator with a square root pole structure such that the phase change across the pole is reduced to π/2. Secondly, in the microscopic theory, frequency-dependent susceptibility terms occur which govern the coupling of light to the condensate. These terms are missing in the classical theory and can, in principle, modify the phase signature of the observed signal. In the clean limit, the influence of the diamagnetic coupling terms cancel out such that the total phase response is dominated solely by the propagator. Yet, in the dirty limit, the paramagnetic coupling terms modify the phase signature and approximately restore the π phase change. Finally, in addition to the Higgs mode response, quasiparticles contribute to the TH signal as well which also show a π/2 phase signature.
Considering long-ranged Coulomb interaction, the TH response of Higgs and quasiparticles is renormalized. Yet, one finds that these modifications do not change the phase signature. Thus, in a BCS superconductor one can expect a π/2 phase change in the clean limit and a π phase change in the dirty limit at the resonance frequency.
If in addition to the Higgs mode other modes exist, the coupled mode scenario with the antiresonance behavior is applicable. Yet, again, microscopic details may modify the classical analysis. To get some insight, we investigate two specific scenarios in this paper: A coupling of the Higgs mode to a charge density wave (CDW) and a superconducting system which additionally hosts a Bardasis-Schrieffer mode, an amplitude fluctuation of a subleading pairing channel.
In the superconductor-CDW scenario, the propagator of the Higgs mode itself is modified as the energy of the Higgs mode no longer overlaps with the total energy gap of the system which consists of the superconducting and CDW gap. Therefore, the propagator obtains a proper pole and the Higgs mode becomes stable. This restores the usual π phase change at the Higgs mode energy. In this system, both modes are driven by light and additionally couple to each other such that an antiresonance effect is expected both from the usual scenario but also due to the interference of the individual contributions. The structure of the analytically evaluated response exactly corresponds to the classic coupled oscillator system. An evaluation of the response as function of frequency and temperature shows the antiresonance behavior and thus, acts as fingerprint of the existence of two modes.
The second scenario with the Bardasis-Schrieffer mode serves as an example of two uncoupled modes. As the Bardasis-Schrieffer mode is a fluctuation of a subleading pairing channel orthogonal to the dominant pairing channel, the coupling element vanishes. Nevertheless, the superposition of the individual contributions leads to an antiresonance behavior where a dip in the intensity and a negative phase change occurs.
The considered scenarios are of course not exhaustive but serve as a proof-of-principle how a careful examination of the THG signal phase gives insight about the existence and nature of collective modes. Furthermore, these scenarios may be applicable in real systems. The superconductor-CDW scenario is relevant for example for NbSe 2 , where superconductivity and CDW coexists [18,26] but might also be relevant for cuprates, where an antiresonance behavior was observed in recent THG experiments [3]. On the other hand, Bardasis-Schrieffer modes might exists in iron-based superconductors. As our study shows, this coupled mode scenario in a microscopic theory is as generally applicable as in the classic theory. Thus we expect that it can describe any other collective mode including Leggett modes in multiband systems [13], Josephson-Plasma modes in layered systems [16], or generally phonon and magnon excitations in the THz regime.
So far, THG experiments on superconductors have only been performed in a setup where the temperature was varied to reach the resonance condition [2,3,36]. As our calculation shows, the resonance and antiresonance signature of the coupled modes is also visible in this case. Yet, to obtain a full picture, it would be necessary to obtain full temperature and frequency data to map out the temperature and frequency dependence of the modes. A further experimental difficulty in this case is the extraction of the phase from the measured signal. As the screening of the THz light is temperature dependent, the first-harmonic signal might be shifted such that a comparison of the relative TH phase could be unreliable. While the phase change is generally more robust than resonance peaks as a signature, strong damping may decrease or wash out the antiresonance behavior, especially if a resonance and antiresonance are positioned close to each other.
To conclude, we have shown that the phase of the THG signal is an interesting quantity to study as it serves as a signature of microscopic details and coupled modes in superconductors. It is a further step in the new field of Higgs spectroscopy and extracting phase information in future experiments will help to reveal more details of the investigated systems.
with electron dispersion k = ξ k − F measured relative to the Fermi level and c † k,σ or c k,σ the electron creation or annihilation operators. We assume that the system is parity symmetric, i.e. k = −k . The separable BCS pairing interaction is given by V k,k = V f k f k with pairing strength V and symmetry function f k . The coupling to light is obtained by the expansion of the minimal coupling term in powers of the vector potential A up to second order Hereby, the paramagnetic term linear in A vanishes due to parity symmetry, i.e. ∂ i −k = −∂ i k , while only the diamagnetic term quadratic in A remains due to ∂ 2 ij −k = ∂ 2 ij k . The partition function of the system is given by with the action in imaginary time τ We decouple the pairing interaction with the help of a Hubbard-Stratonovich transformation. Furthermore, we allow amplitude fluctuations via ∆(τ ) = ∆ + δ∆(τ ). Introducing the Nambu spinor ψ † k = c † k,↑ c −k,↓ , the action can be written in the compact form with the inverse Green's function After integration of the fermions, one obtains in frequency representation where the trace include summation over momentum and frequency and Expanding the logarithm for small Σ, one obtains Relevant for THG is the fourth order action. Thus, we consider the second order term in the sum 1 2 tr G 0 ΣG 0 Σ which leads to with the inverse Higgs propagator and the susceptibilities Finally, integrating the fluctuations using and after analytic continuation iω m → ω + i0 + , one obtains For the following, we will only consider linear polarized light in x-direction A(t) = A 0êx cos(Ωt), such that we can neglect all polarization indices and the action reads where the kernel is given by The third-order current is computed via and is proportional to the fourth-order kernel evaluated at 2Ω. To analytically evaluate the momentum sums, we assume T = 0, s-wave symmetry, i.e. f k = 1, and a constant density of states at the Fermi level such that we can write k → λ d . We use and expand the derivative term k ∂ 2 ij k = k α 0 + α 1 k , which is valid for our band structure. We obtain

Appendix B: Effective action with Coulomb interaction
The susceptibilities for the effective action in Eq. (21) are given in Eq. (A16) and Appendix C: Coupled oscillator The explicit expressions for the complex amplitudesÂ i defined in Eq. (33) are given bŷ We define V 5 = (ω 2 1 − Ω 2 )(ω 2 2 − Ω 2 ) − γ 1 γ 2 Ω 2 − g 2 , V 6 = γ 1 Ω(ω 2 2 − Ω 2 ) + γ 2 Ω(ω 2 1 − Ω 2 ) (C5) to split the nominator and denominator in real and imaginary part Extracting absolute value and phase and using the definitionÂ i = A i e −iφi , finally yields for the real amplitudes A i and phases φ i The total complex amplitudeÂ T is defined aŝ Thus, it follows for the real amplitude A T and phase φ T G −1 0 (k, iω m , iω n ) = (iω m τ 0 ⊗ σ 0 − k τ 3 ⊗ σ 3 + ∆τ 1 ⊗ σ 0 + Dg k τ 3 ⊗ σ 1 )βδ ωm,ωn , where τ i are Pauli matrices in Nambu space and σ i Pauli matrices in the CDW channel. The saddle point equations are A diagonalization of the Hamiltonian yields the quasiparticle energy E k = 2 k + ∆ 2 + |D k | 2 . After integration of the fermions and expansion of the logarithm as in Appendix A, the action is split into a mean-field part and a fluctuation part After evaluating the sum to second order, one obtains the fourth order action Eq (43). The Higgs propagator is defined as Analogously, the renormalized phonon propagator is defined as with tanh(βE k /2) .
Integration of the amplitude fluctuations finally leads to Eq (46).  Figure 10. Phase diagram showing ground state symmetry for system with two possible pairing channels described in Sec. V as function of chemical potential and ratio V d /Vs. In the blue region, the s-wave channel is dominant and in the red region the d-wave channel.

Appendix E: Higgs-Bardasis-Schrieffer model
Using the ansatz in Eq. (54) for V kk including the two pairing channels, the action in imaginary time τ after decoupling of the quartic interaction is given by The usual Higgs mode lives in the τ 1 channel, while the Bardasis-Schrieffer mode lives in the τ 2 channel. In analogy to the previous sections, the fermions can be integrated out which leads to After expansion of the logarithm for small Σ it follows S(δ∆ l ) = S mf + S fl (δ∆ l ) , yy . We evaluate the Bardasis-Schrieffer propagator analytically for T = 0, and in the limit of constant density of state at the Fermi level. We assume f s k = 1 and f d k = cos(2ϕ) and rewrite the sum over momentum as integral . (E24) Using Eq. (A22) one finds