Theory of driven Higgs oscillations and third-harmonic generation in unconventional superconductors

Higgs spectroscopy is a new field which aims for the identification of the gap symmetry by analyzing the Higgs modes in superconductors. One experimental setup in which the Higgs mode in s-wave superconductors was observed is periodic driving with THz light, which shows resonances in the third-harmonic generation (THG) signal if twice the driving frequency matches the energy of the Higgs mode. We derive expressions of the driven gap oscillations for arbitrary gap symmetry and calculate the THG response. We demonstrate that the possible Higgs modes for superconductors with non-trivial gap symmetry can lead to additional resonances if twice the driving frequency matches the energy of these Higgs modes and we disentangle the influence of charge density fluctuations (CDF) to the THG signal within our clean-limit analysis. With this, THG experiments on unconventional superconductors allow for a detection of their Higgs modes and may be used to identify the gap symmetry of materials where the gap symmetry is unknown.


I. INTRODUCTION
In recent years with the emergence of THz spectroscopy, studies on matter were possible in regimes which were inaccessible before [1][2][3] . Due to the energy of THz radiation in the meV range, gentle excitations of materials can be performed without destroying the quantum coherence of the whole system. This allows for controlled experiments in nonequilibrium situations from which ground state properties can be uncovered.
One interesting field within the THz studies are excitations of collective modes like the Higgs mode in superconductors [4][5][6][7][8][9] . In equilibrium this mode arises due to the spontaneous U (1) symmetry breaking in the superconducting state which is characterized by a Mexican hat-shaped free energy potential. In principle there are two different collective modes possible: A phase excitation azimuthally around the brim of the Mexican hat -the Goldstone mode and an amplitude excitation radiallythe Higgs mode. However, due to the Anderson-Higgs mechanism 10 , the originally massless Goldstone mode gets shifted to the plasma energy, whereas the Higgs mode remains with an energy of 2∆, i.e. the energy of the superconducting gap 5 .
The Higgs mode does not possess a dipole moment therefore it does not couple linearly to light. This makes its experimental excitation and detection difficult. Except in the special situation, where a charge-density wave order coexist with the superconducting order on a similar energy range providing a coupling between the two orders and making the Higgs mode Raman-active [11][12][13][14][15][16] , more effort has to be expended for an observation of the Higgs mode. Only recently with ultrafast THz laser technology an excitation and observation of the Higgs mode in a pump-probe experiment for the s-wave system NbN and Nb 1−x Ti x N was possible 17,18 and first experiments on cuprate d-wave systems were reported 19,20 . In these experiments, oscillations of the optical conductivity after a short THz pulse resulting from Higgs oscillations in nonequilibrium could be observed. Alternatively, it was shown theoretically that time-resolved ARPES experiments could be be used in a pump-probe setup as well for the observation of Higgs oscillations 21,22 . These oscillations in nonequilibrium were predicted already before the experiment 23-31 and lead to a great variety of further research considering coupling to phonons 32,33 , quasi-1d systems 34 , multiband superconductors with additional Leggett modes 35,36 , strongly-coupled regimes beyond BCS theory 37 , current-carrying states 38,39 and condensates with d-wave symmetry 40,41 .
Especially the consideration of unconventional superconductors may open a new field of spectroscopy. For superconductors with non-trivial gap symmetry, there are multiple Higgs modes possible which can be decomposed into the irreducible representation of the underlying lattice symmetry 7 and can be understood as asymmetric oscillations of the condensate. Therefore one may infer the symmetry of the gap from a careful analysis of the Higgs modes in a chosen system 41 .
Recent studies on light-induced superconductivity in different systems [42][43][44][45][46][47] raise the question on how to define superconductivity in a short-lived nonequilibrium state. One criteria for superconductivity is the Meissner effect, which is induced on a microscopic level by the Anderson-Higgs mechanism. Hence, a measurement of the Higgs mode is an equivalent probe and it is therefore important to understand Higgs oscillations for all kind of unconventional superconductors.
As pump-probe experiments are difficult to conduct due to the requirement of ultrashort single-cycle THz pulses within the energy range of the superconducting gap, it is worth to consider alternative experimental setups for measuring Higgs oscillations. Instead of quenching the superconducting condensate with a THz pump pulse and observing the intrinsic Higgs oscillations directly, a periodic driving scheme with multicycle pulses can be used. The nonlinear coupling of light to the condensate induces oscillations of the energy gap with twice the driving frequency 48 , which leads to third-harmonic generation (THG) in the induced current and transmitted electric field. If the driving frequency is tuned to the energy of the Higgs mode a resonance in the THG intensity occurs which indicates the existence of the Higgs mode.
This setup was already successfully demonstrated in an experiment for the s-wave superconductor NbN 49,50 . More recently, experiments on several different cuprates were performed, where a new mode at an energy below the symmetric 2∆ Higgs mode was observed 51 . This paper tries to give a possible explanation of such a new collective mode below the energy gap, namely an additional Higgs mode resulting from the nontrivial gap symmetry.
There are several experimental difficulties while performing THG experiments. With current technology, the frequency of the driving field in the THz range cannot be tuned arbitrarily. The experiment is performed in such a way that the driving frequency is fixed and the temperature is changed to obtain the resonance condition.
Second and more seriously is the fact that the Higgs mode is not the only process leading to a resonance in the THG spectrum. Due to the driving, charge density fluctuations (CDF) are excited as well which resonate at the pair-breaking energy 2∆ coinciding with the Higgs mode. Depending on the system considered, these contributions may exceed the Higgs contribution by several orders 52 , however an analysis beyond BCS approximation shows that the situation may reverse again 53 . Studying the polarization dependence of the THG signal can help to distinguish between the contributions, which is however not as clear as to be wished as the polarization dependence depends strongly on the dispersion and pairing interaction 54 . However, other papers pointed out that for superconductors in the dirty-limit the paramagnetic coupling to light plays a crucial role which leads to a dominant Higgs contribution in the THG signal [55][56][57][58][59] . This shows that a careful analysis of an experiment has to be performed in order to interpret the data correctly.
In this paper, we generalize the analysis of 48 within the Anderson pseudospin formulation of BCS theory to arbitrary gap symmetry and calculate the induced gap oscillations and the induced current. As already shown in 60 for a two-band superconductor, the THG intensity contains resonances for each band and for the relative phase oscillation -the Leggett mode. Thus, more complex systems can lead to additional features in the spectrum mapping out the underlying structure.
We show, how an arbitrary gap symmetry affects the THG intensity and that composite gap-symmetries can lead to additional resonances. In addition we show how an asymmetric driving scheme may resonate with additional asymmetric Higgs modes leading to multiple resonances in the THG spectrum. Such a driving scheme was inspired by recent experimental findings 51 and tries to model the asymmetry in driving due to small in-plane components of the external field.
This paper is organized as follows. In Sec. II we introduce the model and explain the time-evolution via Bloch equations. In Sec. III we derive analytic expressions for the driven Higgs oscillations and compare the s-and dwave case. In Sec. IV we calculate analytic expressions for the THG response as a function of the driving frequency. In Sec. V we calculate numerically the temperature dependence of the THG response for s-, d-and d+s-wave. The following Sec. VI evaluates the polarization dependence of the THG signal. In Sec. VII we propose an asymmetric driving scheme and show that an additional resonance peak in the THG signal occurs. Finally we summarize and discuss the results in Sec. VIII.

II. HAMILTONIAN AND TIME EVOLUTION
The starting point of our study is the minimal model for Higgs oscillations, i.e. the mean-field BCS-like Hamiltonian, where we allow a k-dependent energy gap Hereby, k = −k is the energy dispersion measured from the Fermi energy F , c † kσ and c kσ are the electron creation or annihilation operators and the momentum-dependent energy gap, where we assume a separable pairing interaction V kk = V f k f k with strength V and the symmetry function f k such that We implicitly restrict all summations over k in the range −k c < k < k c with the momentum cutoff k c and number of k-points N . For the sake of simplicity we will restrict all of our calculations to two dimensions having also the quasi-2d layered cuprates in mind. Coupling to an electric field is incorporated via minimal substitution k → k−eA(t) with the electron charge e and the vector potential where A 0 is the driving amplitude, Ω the driving frequency and θ the polarization angle. The choice of our vector potential is such that its wave vector is perpendicular to the superconducting plane without an in-plane component. This choice matches the experiments done so far for the conventional s-wave superconductor NbN 49,50 . However, small in-plane components could be induced due to non-perfect alignment, tilted laser pulses or higher-order nonlinear couplings leading to an asymmetric driving scheme. This idea and its effects will be addressed in Sec. VII. For our further analysis we will make use of Anderson pseudospin formalism 61 . We introduce the Nambu-Gorkov spinor Ψ † k = (c † k↑ c −k↓ ) and define Anderson Hereby, we assumed a real equilibrium gap, i.e. ∆ = 0 and ∆ = ∆ . We use the following ansatz to describe the deviations x k (t), y k (t), z k (t) and δ∆(t) from the equilibrium pseudospin expectation values and the gap The time evolution of the expectation values is governed by Heisenbergs' equation of motion which take the form of Bloch equations Note that we have set = 1 here and in the following.
For the further calculation, we expand the z-component of the pseudomagnetic field up to second order in A 0 with the definition where we used the short notation . Due to the assumption k = −k , the term linear in A 0 vanishes and only the quadratic coupling remains. Written out explicitly, we find with The term D k (θ) contains the polarization dependence and the second derivatives of the dispersion. We can see that a non-parabolicity is crucial for the coupling as a parabolic dispersion would lead to a k-independent coupling term describing only a time-dependent variation of the Fermi energy. This changes in the dirty-limit, where also the paramagnetic coupling term remains leading to a much stronger coupling 58,59 . The explicit form of the Bloch equations used in the following reads

III. HIGGS OSCILLATIONS
To gain a first insight into the driven dynamics and see the differences to the known s-wave case, we start by deriving analytic expressions for the gap oscillations. To this end, we assume a small driving amplitude eA 0 1, such that we can neglect terms in second order of the deviation from equilibrium. For simpler notation, we will show only the results for the zero temperature case T = 0. A generalization for finite temperature is straight forward by including the factor tanh(E k /(2k B T )) in the respective expressions and numerical results for this case are shown in the later sections. The linearized equations of motion reaḋ We perform a Laplace transform from time t to complex frequencies s according to to obtain algebraic equations Solving for the deviation terms, we find .
To proceed with a further analytic analysis we use several assumptions and approximations. The x-and y-directions are equivalent for f 2 k and k , i.e. f (k x , k y ) 2 = f (k y , k x ) 2 and (k x , k y ) = (k y , k x ). It follows for any function a k whose k-dependence stems only from f 2 k and k that This assumption is true for functions like k = (|k|) or k ∝ cos k x + cos k y and f k = 1 or f k ∝ cos k x − cos k y . As we are mostly interested in values close to the Fermi energy F , i.e. values of k ≈ 0, we expand the Laplacian of the dispersion to first order in the dispersion As an example, for the single-band tight-binding dispersion on the square lattice with nearest neighbor hopping t it follows and the expansion to first order with α 0 = − F /2 and α 1 = −1/2 becomes exact. Further we assume that the mixed derivatives of the dispersion vanishes, i.e. ∂ 2 xy k = 0. With these assumptions, we can simplify the interaction term A k (t) under a summation which becomes polarization independent To evaluate the momentum sum in the gap equation Eq. (8), we use the approximation that the dispersion depends only on the absolute value of k near F , i.e. k = (|k|), whereas the gap symmetry function only depends on the polar angle ϕ, i.e f k = f (ϕ). Then, we can replace the momentum sum with an integral over the energy and an integral over the polar angle ϕ where the density of states D( F ) is assumed to be constant near the Fermi energy. In the following we will replace the integral borders of the integral with ±∞ which is valid for c ∆. For further notation simplification, we will drop the energy and angular dependencies of the functions in the integrands.
Next, we insert the expressions for x k (t) and y k (t) into the real and imaginary part of the gap equation We start by solving the equation for the real part of the gap. The terms ∝ δ∆ (s) and ∝ α 0 in Eq. (42) vanish as due to the asymmetry of the integrand. The resulting expression reads To obtain this result we identified the equilibrium gap equation by adding a zero 4∆ 2 f 2 +s 2 −4∆ 2 f 2 −s 2 in the nominator of the respective expression in the integrand. The function F (s, ϕ) is defined by Next, we use and define (51) to write δ∆ (s) = α 1 ∆e 2 A 2 0 (I 1 (s) + I 2 (s)). The solution in the time domain is obtained by performing the inverse Laplace transform. For the first term this is trivial and we get For the second term I 2 (t), the Bromwich integral has to be explicitly evaluated The integrand has three poles at s = 0, ±2iΩ and depending on f ∈ [−1, 1] a continuous line of branchpoints between s = −2i∆ and s = 2i∆. We evaluate the integral by extending the path into the complex plane considering the contour in Fig. 1, which is chosen such that the poles contribute with their residue and a branch cut between s = −2i∆ and s = 2i∆ on the imaginary axis is excluded.
The closed loop integral vanishes as no poles are included, the outer integral c 1 vanish if the radius goes to infinity. The integral denoted by c 2 , i.e. the paths from infinity to the poles and between the poles and the end of the branch cut along the imaginary axis cancel each other.
We are left with the contributions around the poles c 3 , c 4 , c 5 and the paths left and right of of the branch cut c 6 and c 7 Each of the small circles around the poles contribute with their residue in the limit of the radius going to zero n=3,4,5 Res r (e st I 2 (s)) .
For the residues it follows The integrals along c 6 and c 7 can be parametrized with Combining the results, we finally get an expression for δ∆ (t) consisting of three contributions Phase ζ(Ω) Driving frequency 2Ω/2∆ Figure 2. Amplitude A and phase ζ of the driven gap oscillation (63). In the case of s-wave, a sharp resonance at 2Ω = 2∆ occurs, whereas for d-wave a broad peak slightly below 2∆ is apparent. In addition, the phase of s-wave shows a sharp jump from 0 for 2Ω < 2∆ to π/2 at 2Ω = 2∆ and a drifting for 2Ω > 2∆. For d-wave, there is a broad phase change from 0 for Ω = 0 up to a value slightly below π/2 for 2Ω = 2∆. At this point there is a sharp kink and a drifting for 2Ω > 2∆ similar to the s-wave case. where δ∆ The term δ∆ D (t) is a forced oscillation of the gap with twice the driving frequency Ω due to the nonlinear A(t) 2 driving. The term δ∆ H (t) is the intrinsic Higgs oscillation with a frequency of ∼ 2∆ induced by an effective interaction quench resulting from the periodic driving 48 . The amplitude of the Higgs oscillation depends on the driving frequency Ω. Without solving the integral explicitely, we can see that there will be a resonance when Ω = ∆ due to the prefactor Ω 2 /(4Ω 2 − r 2 ). We can understand this as a coincidence of the two oscillations 2Ω and 2∆, when Ω is tuned to ∆. This resonance can be found also in the third term, which is again a forced oscillation of 2Ω where the amplitude is also frequency dependent. We rewrite it where A(Ω) represents the oscillation amplitude and ζ(Ω) the phase relative to the driving. Let us first recapitulate the known result for the s-wave case 48 , where f (ϕ) = 1 and λ → λ/(2π). In this case, there are two isolated branch points at s = ±2i∆ which lead to an exact resonance in the amplitude if 2Ω = 2∆, i.e. the driving frequency resonates with the intrinsic Higgs oscillation at 2∆ The amplitude A(Ω) and phase ζ(Ω) of the resonance term is plotted in Fig. 2. One can observe a sharp resonance and also a sharp π 2 phase jump at the resonance condition. Next, let us consider the d-wave case where f (ϕ) = cos(2ϕ). Here, we do no longer have an exact resonance condition but there is still a maximum of the amplitude given by the minimum of the integral in the denominator of (62). This leads to a broad peak in the amplitude accompanied by a smooth phase change of ζ(0) − ζ(2∆) < π/2. The peak in the amplitude is at an energy slightly below 2Ω = 2∆, whereas are sharp kink in the phase at 2Ω = 2∆ is found.
The peak position corresponds to the intrinsic Higgs oscillation with its energy determined by the interplay between the terms f 2 and ∆ 2 f 2 − Ω 2 under the integral over ϕ in (62). For values Ω ≈ ∆, the terms with the highest weights, i.e. f ≈ 1, vanish which decreases the value of the integral and creates its minimum. Therefore, the maximum of the amplitude can be found at 2Ω ≈ 2∆. For values 2Ω > 2∆, the square root is always imaginary for each value of ϕ, which leads to the sharp edge in the phase at 2Ω = 2∆.
After studying the real part, we will finally also calculate the imaginary part of the gap within Eq. (43). Similar to the real part, the terms ∝ δ∆ (s) and ∝ α 1 vanish and we are left with Solving for δ∆ (s) yields and the solution in time-domain is given by We can see that the imaginary part is independent of the gap symmetry. It follows the oscillation of the driving frequency with a frequency independent amplitude. In addition it contains a drift linear in time. Like in the s-wave case, for small driving amplitudes, the imaginary part only contributes to the phase ϕ ∆ of the gap ∆(t) = The absolute value of the gap is independent of the imaginary part and is solely determined by the real part

IV. THIRD HARMONIC GENERATION
The resonance of the forced gap oscillation with the Higgs mode can be found experimentally in the transmitted light field. The nonlinear coupling of the vector potential to the condensate leads to higher-harmonic generation where the lowest non-vanishing order is third-harmonic generation. Its intensity I THG , which is proportional to the amplitude squared of the induced current j (3) (3Ω) shows the same resonance as the gap oscillation and will be calculated in the following. However, not only the driven gap oscillations contribute to THG but also charge density fluctuations, i.e. single particle excitations resonating at the pair-breaking energy 2∆. In the clean-limit, depending on material parameters, these can be significantly larger in an experiment and overlay the contribution from the gap oscillation. Therefore, this contribution is considered and calculated in this section as well and compared to the contribution from the Higgs oscillation.
Driving the superconductor periodically will induce an electric current with the group velocity v k = ∇ k and the charge density To calculate the lowest order response, we expand the velocity The charge density can be expressed with the z-component of the pseudospin and we obtain for the current where The first term j (0) (t) vanishes due to parity, the second term j (1) (t) represents the induced current oscillating with the driving frequency Ω. The third term j (3) (t) is the lowest order of higher-order generation, which oscillates with 3Ω due to the proportionality ∝ A j (t)z k (t) as we have seen in the previous section that z k (t) oscillates with 2Ω. The induced current for an arbitrary angle relative to the polarization of the vector potential can be decomposed into a parallel and perpendicular component which we will calculate separately. We insert the expression for the vector potential and expand the summation of the components To obtain the spectrum, we perform a Fourier transform The expression for the vector potential (4) is inserted and the convolution at ω = 3Ω is evaluated We make use of the solution z k (s = 2iΩ) from the linearized Bloch equations (35) including also the temperature dependence. There are three different contributions originating from the terms ∝ δ∆ (s), ∝ δ∆ (s) and ∝ A k (s) such that we can write j We define and obtain where the real and imaginary part of the gap can be written as and The first term (92) is the Higgs contribution following from the amplitude oscillation of the order parameter. The second term (93) follows from the oscillation of the imaginary part of the gap, i.e. oscillations of the phase. The third term (94) depends on the coupling of the vector potential to the band dispersion and describes charge density fluctuations as the expression has the form of a density-density correlation function 52 .
To understand these expression, we make use of the same approximations as in the previous sections. It follows x ⊥ 4 (s) = 0 , x ⊥ 5 (s) = 0 , (99) x ⊥ 6 (s) = where we identified sin 4 θ + cos 4 θ = 1 − 1 2 sin 2 2θ, 2 sin 2 θ cos 2 θ = 1 2 sin 2 2θ and sin θ cos θ(cos 2 θ − sin 2 θ) = For these expressions we defined the integrals Each of the three terms for the current contribute to the THG intensity Both expressions G(2iΩ) occurring in the phase and CDF terms as well as 1/H(2iΩ) in the Higgs term diverge for 2Ω = 2∆ in the s-wave case or show a maximum in the d-wave case. Therefore, the resonance of the amplitude oscillation with the Higgs mode is not the only cause of the peak found in the THG. In particular, the individual strengths are determined by the parameters λ and α 0 . For α 0 = 0, e.g. the half-filling case in the tight-binding model, the resonance due to the G term vanishes, which removes the phase contribution completely and strongly suppresses the CDF term. For a polarization value of θ = π/4, the CDF contribution perpendicular vanishes and the diverging part due to G in the parallel contribution exactly cancels with the phase contribution. The polarization dependence will be discussed in Sec. VI in more detail. For large λ, i.e. large interaction strength, the CDF contribution is enhanced over the Higgs contribution. It can be understood in this way that the resonance term for the Higgs contribution scales with 1 λ due to the 1 H term, whereas the resonance term for the CDF contribution scales with λ due to the G term. Therefore one can roughly estimate that j In Fig 3, a comparison of the individual contributions for the s-and d-wave case is shown by evaluating the expressions (106), (107) and (109) for typical parameters. We can see that the CDF term exceeds the Higgs term by more than two orders of magnitude. Within the chosen approximations and parameters, the phase term is roughly 1/2 of the CDF term. Both terms show a sharp peak at 2Ω = 2∆ for the s-and d-wave case. This results from the resonance in the G(2iΩ) at the pair-breaking energy 2∆ even for the d-wave case due to the much higher weight at the antinodes relative to the nodes. The shape of the Higgs term, originating from the resonance in the amplitude oscillation 1/H(2iΩ), follows the shape shown in Fig. 2, i.e. a sharp peak for s-wave and a broad peak for d-wave. Again for d-wave, there is also a lot of weight in the range 2Ω < 2∆ as for any Ω there is always a ∆ k where 2Ω = 2∆ k leading to an enhancement of the amplitude. The phase change in the s-wave case is very sharp for all contributions, while for the d-wave case there are small differences. The phase change in the Higgs term is similar to the phase change of the amplitude oscillation, while the phase change of the CDF and phase term is smooth in the beginning but contains a steep step around 2Ω = 2∆.
Despite the fact that the CDF contribution may exceed the Higgs contribution in our simple analysis, it is still useful for the understanding of the physical mechanisms. The actual weighting of the terms in an experiment depends strongly on the material by further effects not considered in our analysis, like retardation effects in materials with phonon-mediated interaction 53 or the paramagnetic coupling for superconductors in the dirty-limit 55-59 .  After gaining a first understanding of the terms contributing to the THG intensity under the chosen assumptions for k and f k , we will drop these in the following sections and solve the summations numerically without approximation and include the temperature dependence. This allows for arbitrary dispersion and gap symmetries to be considered which can introduce new features like additional resonances and polarization dependencies.

V. TEMPERATURE DEPENDENCE
In all current experiments so far, the driving frequency cannot be tuned continuously like done in the theoretical analysis from the previous section. In order to find the resonance, the driving frequency is fixed and the temperature is varied until 2Ω = 2∆(T ) is fulfilled 49 . For the following, we solve the sums (85)-(91) numerically without further approximations. The interaction strength V is calculated for a chosen initial energy gap at T = 0. Then, the temperature dependence of the energy gap ∆(T ) is determined self-consistently for each temperature. To handle the divergences in the summations and keep the required momentum grid resolution in a reasonable range, we introduced a residual broadening of 2iΩ → 2iΩ+0.01∆. This slightly broadens the resonance peaks and washes out the sharp phase jumps but does not change the overall result. We discretize the momentum space around the Fermi energy with an energy cutoff of E c = 2∆ with N k = 2000 points in k and N ϕ = 2000 points in the angular direction. For the following calculation we use a d-wave gap function ∆(ϕ) = ∆ cos(2ϕ) with ∆ = 20 meV and the tight-binding dispersion (38) with t = 200 meV and F = −400 meV. In Fig. 4, the temperature dependence of the THG intensity and phase for d-wave in comparison with s-wave is shown for different driving frequencies. For driving frequencies Ω > ∆(T = 0) no resonance occurs and the THG intensity follows roughly the cubed temperature dependence of the energy gap ∼ ∆(T ) 4 . As soon as Ω < ∆(T = 0), there is a temperature, where Ω = ∆(T ) and a resonance occurs. There are two main differences in the intensities of s-and d-wave. First, for the same driving frequency the resonance peaks are broader for d-wave as there is no single resonance point like in the s-wave case. Second, due to the continuous variation of the gap from 0 to ∆ for d-wave, there is always some resonant gap for temperatures T < T R , where T R is the temperature for which Ω = ∆(T R ). This leads to a larger THG intensity background in that temperature range. One can also see that for the chosen parameters the CDF contribution dominates as there is still a sharp resonance  peak in the d-wave intensity. In comparison, the Higgs contribution alone shows a much broader peak.
The phase analysis confirms the analytic result in Fig. 3. For s-wave a sharp phase jump of π/2 occurs at the resonance, whereas for d-wave the phase jump is much broader.
As an example for what new features a more exotic gap symmetry may introduce, we will study a case with d+s-symmetry, i.e. d-wave with an admixture of s-wave as it was observed for overdoped YBCO [62][63][64][65] . Following 64 , we choose the gap ∆f (ϕ) = ∆ d cos(2ϕ) + ∆ s with ∆ d = 0.9∆ and ∆ s = 0.1∆. The temperature dependence of the THG intensity is shown in Fig. 5. As the absolute value of the d+s gap contains two local gap maxima, i.e. ∆ max 1 = ∆ d + ∆ s = ∆ and ∆ max 2 = ∆ d − ∆ s = 0.8∆, we show the temperature dependence of these two curves. We can see that if Ω < ∆ max 2 , two resonances occur when the driving frequency matches these two maxima. This can be seen both in the THG intensity as two peaks as well as in the phase, where a broad phase transition occurs over the range of the two resonance temperatures with sharp kinks at the resonance points.
It is interesting to note that the two-peak structure is an effect originating alone from the Higgs contribution despite its smaller value. While the CDF contribution only shows a single peak, it is the Higgs contribution which shows the two-peak structure. We can understand this as two Higgs modes at energies ∆ max 1 and ∆ max 2 for each local gap maxima which resonate with the driving frequency. This shows that even if the CDF contribution dominates, the Higgs contribution may still contribute to specific features visible in the spectrum. If the polarization is tuned to θ = π/4, the CDF contribution contains a two peak structure as well as the G term in Eq. (109) is suppressed due to the equivalent term with opposite sign in Eq.(107) and the smaller H term with the two-peak structure becomes visible. We can conclude this section by stating that composite gap symmetries can show additional resonances if there are multiple local gap maxima with different amplitudes.

VI. POLARIZATION DEPENDENCE
As stated above, one possibility in an experiment to gain more insight about the relative weight of Higgs and CDF contributions is the polarization dependence. As one can see in Eq. (109), the CDF contribution has a very characteristic polarization dependence, whereas the Higgs and phase terms do not dependent on the polarization independent on the gap symmetry. If there is no polarization dependence in an experiment it can be a hint that the Higgs contribution is stronger than the CDF part. This was observed in 50 where it was concluded that in NbN the Higgs contribution dominates the THG intensity. If we look at the expressions (109) and (107) for the CDF and phase contribution of the current, we can see that for θ = π/4 the G(2iΩ) term in the CDF expression cancels exactly with the G(2iΩ) term in the phase expression, which means that for this particular polarization angle only the Higgs contribution remains.
However, the derived formulas for the current are only valid under the chosen special assumptions f 2 k and k are symmetric under the exchange of k x ↔ k y , as it is the case for s-or d-wave symmetry and a tight-binding dispersion. For the d+s case this is no longer valid. In Fig. 6 the polarization dependence for s-, d-and d+s-symmetry calculated numerically is shown in comparison with the analytic result. To show the influence of the dispersion, we use a slightly distorted square lattice dispersion in the d+s-wave case with δ 0 = −0.03 and t = 200 meV, t = −80 meV and F = −240 meV 64 . One can see that for s-and d-wave the result closely follows the approximation, i.e. the Higgs contribution is polarization independent and the CDF contribution has a ∝ 1 − 1 2 sin 2 2θ dependency with a small offset resulting from the chosen approximations. This changes for d+s, where the Higgs and phase terms also becomes polarization dependent and the polarization dependence of the total THG intensity deviates from the analytic formula. Choosing a different dispersion can  lead to another polarization dependence. We can see that the polarization dependence alone is not enough to unambiguously distinguish between Higgs and CDF contributions if the band dispersion is not exactly known. This was also concluded in 54 , where it was shown that the Higgs contribution can get polarization dependent as well while on the other hand the polarization dependence of the CDF contribution can be suppressed. Additional deviations from the derived formulas for the polarization dependence are possible, e.g. due to interplay with the pseudogap phase or multiband models. This is however beyond the scope of this paper.

VII. ASYMMETRIC DRIVING
In an experiment, the coupling of light to the superconducting condensate may contain more subtleties and consists not only of a symmetric diamagnetic A 2 term. For examples, small in-plane components of the wavevector induced by non perfect perpendicular alignment of laser and crystal or on purpose tilted lasers as well as higher order couplings to other finite-momentum modes may induce an asymmetry while driving. Such a momentum dependent driving can lead to asymmetric oscillations of the condensate with respect to the origin. One finds 41 that an asymmetric oscillation of the condensate can show up as a second frequency in the gap oscillation below the well-known 2∆ Higgs mode. This asymmetric Higgs mode depends on the gap symmetry and the respective asymmetric deviation.
We propose a phenomenological asymmetric driving scheme to describe such effects in an experiment. Due to a momentum dependent driving, the gap symmetry gets altered by an additional symmetry component f A k . We add such a term to the pseudomagnetic field altering the gap symmetry with the same time-dependence as the usual driving term where δ A determines the strength of the asymmetric driving. There is some similarity to 7 , where, however static and ab-initio, a composite pairing interaction leads to multiple Higgs modes in the different pairing channels. Our approach is purely phenomenological and neglects additional polarization dependencies which are likely to occur due to the asymmetric driving. As the polarization dependence is hard to predict without an in detail understanding of the actual microscopic coupling, we neglect it for our approach.
Using the modified pseudomagnetic field and performing a linearization in the same way as in Sec. III one would neglect important contributions from the products between the asymmetric driving term and the deviations, e.g. δ sin 2 (Ωt)f A k x k (t), etc. Therefore, we solve the Bloch equations for this section numerically without any approximation by integrating the differential equations in time. As we haven't added an additional polarization dependence, we can suppress the CDF contribution by choosing θ = π/4 to obtain the contribution from the Higgs channel.
In Fig. 7, we show the temperature dependence of the THG intensity for a d-wave gap with an asymmetric driving f A k = 1, i.e. a distortion in the s-wave channel. This corresponds to driving the B 1g mode of the d-wave gap. For low driving frequencies one can observe a second resonance peak in the THG signal below the 2∆ peak, which is also accompanied by a phase change of π/2. Higher driving frequencies do not show such a resonance. This can be understood better, by calculating the frequency dependence of THG for different temperatures from which one can extract the temperature dependence of the modes. For the chosen set of parameters, the resulting curve of the second mode follows approximately ω B1g (T ) = 0.44∆(T ) which is also shown in the upper row of Fig. 7. As one can see, the chosen higher driving frequency of 2Ω = ∆ 0 stays always above the second mode which explains that there is no second resonance peak. For the lower driving frequency of 2Ω = 0.4∆ 0 , the second resonance peak appears at the point where 2Ω = ω B1g (T ) Indeed, in recent THG experiments on cuprates a collective mode below 2∆ was observed 51 . An asymmetric driving and resonantly excitation of a B 1g mode (or other non-A 1g modes) could be an explanation of this experiment.

VIII. SUMMARY AND DISCUSSION
In this work we analyze the induced gap oscillations in unconventional superconductors due to periodic driving with THz light within BCS theory in the Anderson pseudospin formalism. We show exemplary for d-wave how a non s-wave symmetry broadens both the resonance peak and the phase change in the oscillation amplitude and derive analytic expressions for the oscillations. The gap oscillations lead to third-harmonic generation whose intensity shows a resonance peak if twice the driving frequency coincidences with the energy of the Higgs mode. We compare contributions from CDF and Higgs and illustrated that the Higgs contribution to THG follows closely the amplitude of the gap oscillation, whereas the CDF contribution shows still a sharp peak even in the non s-wave case. The temperature dependence of the THG signal as measured in the experiment shows the same characteristics as the THG signal as a function of the driving frequency.
In addition one finds that more complex gap symmetries, shown exemplary for d+s-wave, which contain multiple local gap maxima show resonances for each maxima. Hereby it is the contribution of the Higgs oscillation which shows this feature as the composite gap contains two Higgs modes with frequencies at these local gap maxima.
For d-wave, or in general for gap symmetries where the squared symmetry function is isotropic in x-and y-direction, the polarization dependence of the THG signal is the same as for s-wave. For a simple squarelattice tight-binding dispersion, the Higgs contribution is polarization-independent, whereas the CDF contribution has a characteristic dependency. However, deviations both from the isotropy of the gap symmetry as well as the simple tight-binding dispersion lead also to deviations in the polarization dependence and the Higgs contribution can become polarization dependent. Yet, if the banddispersion is known, a measurement of the polarization dependency allows to distinguish the individual contributions from CDF and Higgs.
In our analysis for the clean-limit, where the coupling to the vector potential is exclusively through the nonparabolicity of the band dispersion, the CDF contribution exceeds in general the Higgs contribution. Nevertheless, this work gives an interesting insight into the driven dynamics of unconventional superconductors. On the one hand, the Higgs contribution can still show its unique features as seen in the d+s-wave case with the two-peak structure, which is a result from the Higgs contribution alone. On the other hand as it was shown elsewhere 55-59 , effects beyond BCS theory and impurity scattering in the dirty-limit can strongly enhance the Higgs contribution over the CDF contribution.
Finally, we propose an asymmetric driving scheme to model experiments where the coupling of the driving field acts non-symmetrically with respect to the groundstate symmetry on the condensate. Such asymmetry induces asymmetric oscillations of the condensate which can show up as an additional oscillation frequency of the gap, dependent on the symmetry of the gap and the deviation. With such a setup, THG experiments should be in principle able to measure asymmetric Higgs modes and provide therefore the same information as pump-probe experiments. However, as the proposed phenomenological asymmetric driving scheme may be difficult to realize experimentally, the information obtained by THG experiments are limited in this respect.
In a more general sense, THG experiments may also serve as a new measure for defining superconductivity in nonequilibrium. Recent experiments and theoretical studies on light-induced superconductivity 42-47 raise the question on how one defines superconductivity in a shortlived nonequilibrium state. So far, the criteria only include the vanishing resistivity property of superconductors measured by a divergent imaginary part of the optical conductivity for ω → 0, however the expelling of a magnetic field, i.e. the Meissner effect, has not yet been considered. As the Meissner effect is induced by the Anderson-Higgs mechanism, a measurement of the Higgs mode should be an equivalent fingerprint of superconductivity. While the repulsion of a magnetic field on an ultrashort timescale is difficult to measure or even impossible, a resonant behavior of the THG signal in the light-induced superconducting state could potentially be realized.
In addition to pump-probe experiments, where the gap is quenched by a short pulse and the following intrinsic Higgs oscillations can be observed, THG experiments can serve as an alternative tool for identifying Higgs modes of a superconductor. These driven experiments have some advantages over the pump-probe experiments. No ultra-short single-cycle pulses are required and the strongdamping of Higgs modes in gaps with nodes are partly overcome due to the forced periodic driving and oscillation of the gap. Thus, in the context of Higgs spectroscopy, i.e. the classification of gap symmetries by identifying Higgs modes in a superconductor, THG experiments extend the range of possible experimental setups.