Calculation of an enhanced A1g symmetry mode induced by Higgs oscillations in the Raman spectrum of high-temperature cuprate superconductors

In superconductors the Anderson-Higgs mechanism allows for the existence of a collective amplitude (Higgs) mode which can couple to eV-light mainly in a non-linear Raman-like process. The experimental non-equilibrium results on isotropic superconductors have been explained going beyond the BCS theory including the Higgs mode. Furthermore, in anisotropic d-wave superconductors strong interaction effects with other modes are expected. Here we calculate the Raman contribution of the Higgs mode from a new perspective, including many-body Higgs oscillations effects and their consequences in conventional, spontaneous Raman spectroscopy. Our results suggest a significant contribution to the intensity of the A1g symmetry Raman spectrum in d-wave superconductors. In order to test our theory, we predict the presence of measurable characteristic oscillations in THz quench-optical probe time-dependent reflectivity experiments.


Introduction.
Collective excitations of superconductors in non-equilibrium are a new emerging field. The most studied collective modes in superconductors, among others, are the amplitude (AM, so-called Higgs mode) and phase (Goldstone) modes [1,2]. Yet, while the Higgs is a collective oscillation having a frequency of twice the superconducting order parameter ∆, i.e. ω Higgs = 2∆, the phase mode is shifted to the plasma frequency at higher energies due to the Anderson-Higgs mechanism [1]. Since the Higgs mode does not carry a dipole moment, it is well known that a linear activation by light is almost impossible [3,4], such that the relevant process is the non-linear Raman effect. In ultra-fast non-equilibrium experiments the Higgs mode has been uncovered both in quench-probes [5] and in periodically driven setups [6][7][8], respectively, where the coupling of light to the charges of the superconductor can be described by a quadratic Raman-like process. In this context, the third harmonic generation has been a widely studied example of this non-linear optical effect [9][10][11]. In the '90s conventional, steady-state Raman spectroscopy on superconductors turned out to be an effective experimental technique to provide evidence via polarization-dependence for d-wave Cooper pairing [12][13][14][15][16][17]. In high-T c cuprates with D 4h crystal symmetry the three most important Raman symmetries (A 1g , B 1g , B 2g ) could be explained by a d-wave order parameter [17][18][19][20][21]. However, in all spectra of cuprates there is still an open question concerning the A 1g peak intensity [21,22]. In conventional Raman scattering theory, due to Coulomb screening of quasiparticles (QP) the response in A 1g symmetry is much smaller than the experimental observations in cuprates [16,21], which is of the same order or larger than the B 1g peak [17][18][19]. This is the so-called long-standing A 1g problem [21]. Over the past three decades, various attempts have been made to reconcile theory with experiments without arriving at conclusive results: these include, but are not limited to, the use of higher harmonics [18,19], magnon resonance [23,24], screening effects [25], neutron magnetic resonance [26] and the interaction of the AM with the η−mode, a spin-singlet excitation which can appear below 2∆ in d-wave superconductors [27]. However, none of these approaches could solve the A 1g problem quantitatively, namely providing the correct intensity and peak position. Recently, Puviani et al. [28] have shown that both equilibrium and non-equilibrium activation of the Higgs oscillations in superconductors correspond to the same physical non-linear Raman process. The conventional equilibrium Raman response of s-wave superconductors should be small for symmetry arguments, with the exception of the presence of competing orders, i.e. superconductivity and charge density wave [29,30]. Eventually, recent terahertz pump-optical probe (TPOP) experiments have provided hints of the presence of the Higgs mode through transient reflectivity change in the d -wave superconductor Bi2212 [31,32]. Despite this, no effective time-dependent oscillation detection, nor a consistent theoretical calculation have been provided for d -wave superconductors to support this. In this Letter, we go beyond the present theoretical background on the amplitude mode in superconductors, including many-body effects of Higgs oscillations on quasiparticles: when light interacts with the charges of the condensate creating an electron-hole pair, these particles undergo the so-called Andreev scattering which leads to creation of pairs of electrons/holes which add to the existing condensate and eventually produce the Higgs (Fig.  1a). However, in addition to this production mechanism, which is generally subdominant, the quasiparticle excitation generated by light can lead to dominant many-body Higgs oscillations due to a mechanism of interaction between Cooper pairs (CPs) generated by Andreev scattering (Fig. 1b), which has been neglected so far. Theoretically, this is done through a systematic diagrammatic approach, adding vertex corrections in the amplitude channel alongside with the usual random phase approximation (RPA) summation for the AM propagator. This calculation gives rise to an effective hybridization of the amplitude and charge degrees of freedom, namely many-body effects of Higgs and quasiparticles, resulting now in a new strong contribution of the Higgs mode. This approach is able to show that there is a significant signature of the Higgs mode in the calculated Raman spectra of d -wave superconductors, as suggested in Ref.s [33,34], generalizing the phenomenological coupling used by Zeyher and Greco in Ref. [35]. In our theory a qualitative and quantitative agreement with the Raman experimental data is possible. Simply speaking, similar to phonons which can be renormalized to polarons due to strong electron-phonon interactions, the Raman-active Higgs mode can be strongly renormalized due to quasiparticles to form a Cooper pair-polaron, resulting in a dressed Higgs mode with strong intensity in A 1g polarization. Note that the conclusions of our work can be generalized to other clean systems (e.g.: SC+CDW, multi-band superconductors) where the Higgs peak is shifted from the quasiparticle continuum and the symmetry allows for a many-body response. In order to further test our theory, we suggest that using an ultra-short THz pump in TPOP experiments our theory would explain the relative intensity of the transient reflectivity change measured for the A 1g and B 1g symmetries, and it predicts the possibility to detect the characteristic time-dependent oscillations due to the amplitude mode contribution in both of these symmetries.
Calculation of the amplitude mode. For describing a d -wave superconducting system, we consider a generalized version of the BCS theory, namely the Hamiltonian: H = H el. +H p +H C . It is the sum of the electronic tightbinding Hamiltonian H el. , a momentum-dependent attractive pairing interaction H p and a long-range Coulomb repulsive energy H C , respectively. The pairing term responsible for the AM is given by with the pair creation operatorP † k,q =ĉ † k+q,↑ĉ † −k,↓ , and the pairing interaction function V (k, k ) [9,10]. For dwave superconductors we can factorize the pairing inter- is the d x 2 −y 2 -wave form factor, V being the pairing strength derived from the self-consistent gap equation.
The Hamiltonian H el. + H P can be treated at mean field level defining the order parameter ∆ k = ∆ max f k , giving rise to a mean-field Hamiltonian for unconventional superconductors. The bare susceptibility response of the condensate coupling quadratically to light is given by where τ i=1,2,3 are the Pauli matrices, τ 0 the identity matrix, γ i,s k the interaction vertex for the incident/scattered light for a given symmetry, G k ,m ≡ G(k , iν m ) and G + k ,m ≡ G(k + q, iν m + iω) the matrices of the Matsubara Green's functions in Nambu-Gor'kov space. The coupling between the light and the amplitude mode occurs indirectly via the charge-amplitude susceptibility, χ γf (q, ω). Then, the AM propagator can be calculated with the RPA summation as in Ref. [10] Here we have introduced the susceptibility χ f f (q, ω) with vertices f k τ 1 −f k τ 1 in the amplitude channel. This allows us to calculate the total Raman susceptibility including both quasiparticles (Eq. (2)) and amplitude mode as where we used the susceptibilitiesχ which include charge fluctuations due to the long-range Coulomb interaction between quasiparticles [36]. We want to stress here that for a given value of the order parameter ∆ max we calculate self-consistently within the gap equation the value of pairing strength V. Note that here we have considered only the diamagnetic susceptibility, since the paramagnetic contributions are vanishing in clean superconductors and in visible light experiments, as well as coupling to the spin channel (resonance peak) is negligible [36].
Many-body Higgs oscillations. As mentioned earlier, we now want to go beyond the RPA summation for the amplitude mode propagator, including many-body dynamic Higgs interactions between Cooper pairs generated by electrons/holes via Andreev scattering ( Fig. 1  b), mixing the amplitude and charge channels. To do that, we defined a dressed light-condensate interaction vertex through the ladder diagram summation Analogously, for the f k τ 1 vertex in the amplitude channel we derived the self-consistent sum Substituting to the bare vertices the dressed ones including the many-body effects of the Higgs mode, the Raman response in Eq. (4) can be replaced with the full manybody expressioñ where the susceptibilitiesχ include charge fluctuations [10,16].
Steady-state Raman spectroscopy. The intensity measured in Raman experiments is given by the photon scattering differential cross section, which is proportional to the imaginary part of the susceptibility, namely ∂ 2 σ/∂ω∂Ω ∝ −χ (q, ω). In order to solve the set of equations for the dressed vertices and the susceptibilities, we considered the lowest terms in a Fermi-surface harmonic expansion of the Raman vertices, restricting ourselves to q = 0 for the light and to the clean limit of superconductor, in agreement with previous studies [18,19]. Using the same optimized parameters fixed for all the Raman symmetries investigated, we additionally checked systematics in order to confirm the stability of the results [36]. In Fig. 2a we show the effect of the bare AM and the many-body dressing contribution in the Raman spectra for different symmetries of Bi2212 (a d -wave superconductor with 2∆ max ∼ 550 cm −1 ): on the one hand, for the A 1g response the dressed AM calculated within the RPA approximation using Eq. (4) is negligible, while the full many-body Higgs oscillations result into a contribution which is more intense and peaked. On the other hand, for the B 1g and B 2g symmetries the charge fluctuations and Higgs oscillations screen the Raman response, reducing its intensity. In particular, the B 1g spectrum becomes broader in frequency and the peak gets shifted to higher frequencies, while the B 2g acquires a shoulder peaked at around ω ∼ 2∆ max . These many-body Higgs effects describe the experimental Raman on Bi2212 in a much improved fashion, thus contributing to the long-standing problem regarding the A 1g peak intensity [18,19]. As we have concluded before, the A 1g response provided by QPs is much lower than the experimental one, and the bare AM contribution is negligible. However, considering many-body effects of the Higgs oscillations, the B 1g (Fig.  2b) and B 2g (Fig.  2c) spectra get screened, while the A 1g response becomes enhanced and peakes now at a frequency ω ∼ 400 cm −1 , in quite good agreement with the experimental result. Note that a simultaneous and quantitative understanding of all relevant polarizations is indeed possible employing the same set of many-body parameters [36] as shown in Fig. 2d. This contributes to explain the long-standing A 1g problem concerning the peak intensity in optimally doped Bi2212 as formulated in the literature [17,18,21].
Ultra-fast Raman-like optics. As discussed earlier, the same Higgs oscillations should be seen in non-equilibrium in a time-dependent experiment. While, in principle, it is possible to measure Fig.  2d time-dependently in an ultra-fast experiment, it is much easier to detect the Higgs oscillations in the transient reflectivity change ∆R(t)/R. It has been widely demonstrated that the coupling of the light to the collective AM can happen only in a nonlinear regime with a Raman-like interaction, thus new evidences should be looked for in nonlinear optical effects, such as the Kerr signal in pump-probe experiments. To this extent, it is useful to look at the timedependent reflectivity change in non-equilibrium superconductors, which is related to the equilibrium Raman susceptibility aforementioned and easily accessible in experiments. The transient reflectivity change due to a pump electric pulse E pump (t) is given by [37] ∆R where is the real part of the third order susceptibility in real time, t being the time at which the reflectivity change is detected, t < t the previous time at which the pump pulse interacts with the condensate. Recently, Katsumi et al. [31,32] demonstrated that the relative intensity between the transient reflectivity change projected onto the symmetries A 1g and B 1g in a TPOP experiment on Bi2212 cannot be explained only with the quasiparticle density fluctuations contribution, deducing that the Higgs should play an important role. However, since they pumped the system in the fully adiabatic regime (ω pump ∆ max , ∆t pump ∼ 4 ps), they were not able to probe any Higgs oscillation in the time-dependent response. We performed the calculation of the theoretically expected time-dependent transient reflectivity change on a d -wave superconductor, using Eq. (8) together with the third-order susceptibility containing the full many-body Higgs contribution according to Eq. (7). Considering the same pump regime of the experiments, we have been able to get a good agreement of the relative intensities for the A 1g and B 1g symmetries [36]. In addition to this, we repeated the numerical calculation with an ultra-short THz pump (duration ∆t pump ∼ 0.1 ps, frequency ω pump = 5.57 THz < ∼ ∆ max ) in order to simulate a quench experiment: the result is shown in Fig. 3. Including the full Higgs contribution to the third-order susceptibility, oscillations with a frequency ω = ω Higgs = 12.5 THz < ∼ 2∆ max appear in the decaying region of ∆R(t)/R, while a lower intensity and no oscillations are present considering only the QPs and The dark red and dark blue curves represent the normalized transient reflectivity of the A1g and B1g symmetries, respectively, calculated including the full many-body Higgs contribution. The light red curve, instead, shows only the QP and CF response in the A1g symmetry. The top right inset shows the differential A1g (red line) and B1g (blue line) signal due to the Higgs contribution. charge fluctuations contributions. Performing such a time-resolved experiment would provide a clear evidence of the presence of the Higgs mode and the effects of many-body contributions in d -wave superconductors: indeed, they are responsible for enhancing the A 1g symmetry response and, at the same time, for the screening effect on the B 1g .
Conclusion and outlook. For the first time, we pro-vide calculations of Higgs oscillations for d -wave superconductors including many-body contributions originating from the interaction with quasiparticles (e.g. broken Cooper pairs), resulting in a renormalized, dressed Higgs mode. We treat the full vertex corrections within a diagrammatic approach, which is able to describe both nonequilibrium and conventional steady-state Raman experiments, respectively. In particular, we are able to provide a contribution to the long-standing A 1g problem in Raman spectra of d-wave superconductors [18,19] in Fig.  2d. Our conclusions are not based on any specific scenario for cuprates. Furthermore, our results are also unchanged if we take into account that Higgs oscillations can be dressed by phonons as well. However, these contributions to the Raman spectra are minor [36]. Furthermore, we explain the non-equilibrium response and the correct magnitude of the transient reflectivity change in recent THz pump-optical probe experiments [31]. Our theory can simply be tested by a smoking gun experiment: we propose an experimental set-up using ultrashort THz pump pulse in a non-adiabatic regime to measure the reflectivity in d -wave superconductors, in order to quench the condensate and thus directly observe the time-dependent oscillations characteristics of the Higgs mode. It is expected that the resonance peak obtained with inelastic neutron scattering can play a role in the A 1g problem because of the energy matching with the A 1g Raman response. However, it is interesting to note that there is no linear coupling of the spin channel to the Higgs mode or to the Raman vertex: as shown by Venturini et al. in [23], only including spin fluctuations at higher order it is possible to match the energy position of the A 1g peak, but without substantially affecting its intensity. We also point out that so far it has been possible to clearly detect the Higgs contribution only in the presence of competing orders such as in 2H − NbSe 2 , namely superconductivity and charge density wave. The latter provides the excited phonons which couple to the superconducting AM [38]. This mechanism allows to shift the Higgs A 1g peak to energies lower than 2∆, becoming thus unambiguously detectable [38][39][40][41]. Similarly, a Higgs response has been reported in E 2g Raman response of pressure-induced 2H − TaS 2 [30], a transition metal dichalcogenide with cohexisting SC and CDW. Our findings provide a new route for calculations of many-body effects of the Higgs mode while interacting with various other modes in superconductors [39,42]. We believe that an extension of our model to these systems can provide further insights and guide future experiments on unconventional superconductors and systems with cohexisting phases.  We provide here some definitions and short notation that is used throughout the text. First, we define averages g(k) in terms of the Tsuneto function: where iω n is the Matsubara imaginary frequency, ε k the electronic band structure, E k = ε k 2 + ∆ k 2 the quasiparticles energy, while the gap function is ∆ k = ∆ max f k . Then, we can introduce the Nambu-Gor'kov spinor Ψ † k = (ĉ † k,↑ĉ −k,↓ ), and distinguish the interaction terms the phase channel and eventually the charge channel is given by Ψ † k τ 3 Ψ k =ĉ † k,↑ĉ k,↑ −ĉ −k,↓ĉ † −k,↓ . Thus, we write the bare BCS Green's functions in Nambu-Gor'kov space as We can usefully rewrite it in its spectral form [1] where G (k, ω + iδ) is the imaginary part of the Green's function. In order to shorten the notation , we define G k,m = G(k, iν m ) and G + k,m = G(k + q, iν m + iω n ). The function in real frequency space is obtained with the analytic continuation iω n → ω + iδ. Finally, we define the susceptibilities according to the interaction vertices, so that χ τiτj are the polarization bubbles with τ i −τ j interaction vertices, while χ γτi is the γ k τ 3 −τ i bubble and χ γf is the susceptibility with vertices γ k τ 3 −f k τ 1 . Analogously we can write all the other susceptibility amplitudes.

B. Bare Raman response
For describing the d -wave superconducting phase, we consider a model Hamiltonian which can be split into three components, namely H = H el. + H p + H C , which is the sum of the electronic tight-binding term, attractive pairing arXiv:2012.01922v2 [cond-mat.supr-con] 4 Oct 2021 interaction and long-range Coulomb repulsion, respectively. The tight-binding term is simply where ε k is the band dispersion; while the Coulomb interaction is given by The pairing term responsible for the amplitude mode is given by with the pair creation operatorP † k,q =ĉ † k+q,↑ĉ † −k,↓ , and V (k, k ) = −V f k f k the factorized pairing interaction function [2,3] calculated as A simple choice for the tight-binding band structure to be used in (4) which contains the basic physics of 2D-like tetragonal systems is given by the following: where t is the nearest-neighbor hopping integral, 2B = t /t the ratio between the nearest-neighbor and the nextnearest-neighbor hopping, µt the chemical potential. The d -wave order parameter is ∆ k = ∆ max f k and the corresponding form factor results to be f k = (cos k x − cos k y )/2. Within the effective mass approximation, the Raman vertex can be written in terms of the curvature of the energy band according to the expression γ k = k e i ∂ 2 ij ε k e j , with ∂ 2 ij = ∂ 2 /∂k i ∂k j , e i,j being the polarization vectors of the incident and scattered light, respectively [4,5]. The derivatives of the band dispersion (8) with respect to k x and k y are ∂ x ε k = 2t(sin k x − 2B sin k x cos k y ) , ∂ y ε k = 2t(sin k y − 2B cos k x sin k y ) , and the second derivatives ∂ 2 xx ε k = 2t(cos k x − 2B cos k x cos k y ) , ∂ 2 yy ε k = 2t(cos k y − 2B cos k x cos k y ) , ∂ 2 xy ε k = 4tB sin k x sin k y .
Therefore the bare Raman vertex functions for different symmetries are given by The bare Raman diamagnetic bubble with the charge-charge vertices which probes the quasiparticles (QPs) response can be written as a sum over all the Brillouin zone and the Matsubara frequencies as where γ i,s k is the Raman interaction vertex for the incident/scattered light of the appropriate symmetry. Making use of the spectral representation in (3) we can write the Matsubara summation in (10) as In the limit q → 0 we obtain the Raman susceptibility which can be written using the compact notation in (1) in the short form Following [4], we can include the Coulomb interaction in the long-range limit giving rise to charge fluctuations (CFs) by means of the expressionχ with C. RPA amplitude mode propagator The pairing interaction energy in (6) can be included by means of the RPA summation of the amplitude mode (AM) propagator where χ τ1τ1 (q, k , ω) is the τ 1 − τ 1 bare polarization bubble. Using the usual factorization for the pairing interaction function by means of the form factor f k , we can write D AM (q, k, k , ω) = D AM (q, ω)f k f k with Here we used the pairing interaction coupling in the amplitude channel with strength V /2 [3], and the f k τ 1 − f k τ 1 susceptibility amplitude The Raman response contribution which includes the RPA AM propagator in (18) can be calculated as with the susceptibility χ γf (q, ω) = k f k γ i k χ τ3τ1 (q, k , ω), and the γ k τ 3 − f k τ 1 polarization bubble In the limit q → 0 we simply get: Adding the CFs [3] in Eq. (20) we obtaiñ Definingχ we can rewrite the expression in (23) asχ

D. Many-body Higgs oscillations
We now consider the many-body Higgs oscillations (HOs) by means of the vertices dressing with the pairing interaction energy terms: this allows for a mixing of the interaction channels described in section I.-A. The dressed γ k τ 3 can be calculated with the Bethe-Salpeter equation in the rainbow-ladder approximation: The simplest form that the dressed vertex function can have in order to satisfy this self-consistent equation is given by  (29), respectively. With the replacements γ k τ3 → Γ and τ1 → Λ in Fig. 1(b) the full many-body Higgs oscillations contribution to the Raman response is obtained.
where we have neglected the spin and phase channels and putting Γ .

(28b)
Analogously, for the dressed τ 1 vertex we adopt the form: and solving the corresponding self-consistent equation we obtain the values .

(30b)
Thus, the full Raman response obtained using the dressed vertices is written as where the susceptibilities screened by the long-range Coulomb interaction are given bỹ We used the definitions of the unscreened susceptibilities as follows: Their full expressions in term of averages are given by:

A. Phonon-Higgs Raman response
Let's consider a phonon with q = 0 and frequency ω 0 : its bare propagator is given by Including the screening due to the RPA electron-phonon interaction, with g D (k) = g 0 g k being the electronphonon coupling function specific of a given phonon symmetry, we get the dressed propagator D ph (Ω) = D 0 with This can be calculated for the superconducting state, obtaining Below T c , we can also include interaction with the AM propagator, obtained via RPA summation,D(Ω) = D ph (Ω) + with The phononic Raman susceptibility is given by with with γ k the Raman vertex. We can also add the explicit interaction of the phonon with the AM in the Raman response (Fig. 3 c-e), obtaining the total phononic contribution We can now include the CFs due to Coulomb interaction in the long wavelength limit (q → 0 in the Coulomb potential 4πe 2 /q) using the following functions rather than the bare ones: (50)

B. Phonons with many-body Higgs interactions
We now include the interaction of the phonons with many-body Higgs oscillations by means of the vertex corrections in the very same way as used for the Raman response calculated in Section I.D. Thus, the electron-phonon interaction vertex g k for a given symmetry gets dressed Then we rewrote Eq. (45) using the dressed vertices with many-body Higgs interactions, namely obtaining where we also took into account of the CFs due to Coulomb long-range interactions into the expressions of the susceptibilities and the phonon and AM propagators.

III. INTEGRATION AROUND THE FERMI SURFACE
The tight-binding band dispersion used so far is derived from a reduction of a three band model which gives a large contribution to the density of states at the Fermi level for cuprates. However, in order to take into account of admixtures of more Fermi surface harmonics to the Raman contributions for different symmetries, we adopted a calculation independent of the details of the band dispersion, namely integrating the response functions around the Fermi surface and including higher harmonics.

A. Integrals around the Fermi surface
Let's assume that the energy band dispersion ε k depends only on the absolute value of k, while the form factor given by the gap symmetry depends only on the polar angle f k = f (ϕ) (and thus ∆ k → ∆(ϕ) = ∆ max f (ϕ)). We can replace the summations over the points k to integrals over the energy ε (around the Fermi level in the interval [−ε c , +ε c ], ε c being the appropriate cut-off energy) and the azimuthal angle ϕ: S where D(ε) is the density of states, such that +εc −εc D(ε) dε = 1. We assume that the density of states is constant around the Fermi energy, D(0) ≈ const. = 1/(2ε c ). Thus we can write We can define λ ≡ V D(0) to shorten the notation. In general we can substitute to the discrete sum over the k points the integral k → D(0) So that we obtain (57) We observe that this value depends upon the energy cut-off ε c , so that λ = λ(ε c ). In the limit ε c ∆ max we can substitute the integral +εc −εc dε with +∞ −∞ dε. Instead, if we consider only the Fermi surface, i.e. the density of states being D(ε) = δ(0), then we would get which does not depend upon any cut-off. We introduce now for the Fermi-surface integration approach the definition of the average x analogous to the one in (1) as follows: If x depends only on the polar angle, then we can solve the integral over the energy and get

B. Raman response
We consider a d-wave superconductor with D 4h lattice point group symmetry, and we express the gap function depending only on the azimuthal angle ϕ as ∆(ϕ) = ∆ max f (ϕ), with the form factor f (ϕ) = cos (2ϕ). We first consider a cylindrical Fermi surface, and then we allow for deviations from cylindricity retaining the lowest terms in a Fermi-surface harmonic expansion depending on azimuthal angle ϕ [5]: Where α 0 , α 1 , α 2 and β 1 = γ 1 are parameters related to the material. If α 0 = 0 we recover the result of the integral over the Fermi circle: in this case no AM is present, because all the integrals over ε are odd and go to zero. If we set α 1 = 0 and α 0 = −1/2 (which is the case for a square lattice, or B = 0 referring to the tight-binding model in (8)), then we are in the circular Fermi surface limit. In the general case with α 0 = −1/2 and α 1 = 0 we obtain the full result [2,6]. We can also add higher harmonics to the A 1g vertex [5], obtaining the final expression: The bare A 1g Raman response without AM and with CFs due to long-range Coulomb interaction is given bỹ Substituting the expressions for the vertex and the gap function, retaining only the nonzero integrals we get The amplitude mode contribution to the A 1g response calculated within the RPA is given by The B 1g and B 2g responses are calculated as

C. Many-body Higgs contribution
Considering the many-body effect of Higgs oscillations through the vertex corrections, for the A 1g symmetry we get the susceptibilities: The vertices are given by while Λ 0 = 0. On the other hand, for the B 1g response we have: and Γ (1) Here we used the value of λ as given in Eq. (57), with the pairing interaction strength λ/2. For the B 2g Raman susceptibility it is enough to replace β 1 f (ϕ) → γ 1 sin(2ϕ) in Eq. (70). Then, we can add the contributions of CFs, analogously to the previous case.

IV. IMPURITIES AND eV -LIGHT RAMAN
As shown by N. Tsuji and Y. Nomura [7], the resonant paramagnetic superconducting response which is activated by impurities is given by the following diagrams: Without writing explicitly the k-dependence of the propagators and the momentum dependence due to the impurity scattering, the quasiparticles' paramagnetic susceptibility is given by the expression where Ω eV and ω R are the eV-light frequency and the Raman frequency shift, respectively. This susceptibility vanishes exactly in the superconducting clean limit [3], while it is activated (i.e.: becomes non-zero) in the presence of impurities. The Green's function expressed in a matricial form in Nambu-Gor'kov space reads: where τ i are the Pauli matrices. From this expression we can realize that the paramagnetic process is in an off-resonant condition, since the pole of the fermionic propagator G is shifted up to the eV range by the light energy Ω eV , far away from the resonance of the superconducting gap. The same arguments are valid for the diagram inclluding the Higgs mode, which involves the susceptibility χ jjf . Therefore, in the non-resonant eV-light Raman processes the impurity contribution [], which comes mainly from the processes described which are off-resonant in this regime, can be fairly neglected.

V. SPIN SUSCEPTIBILITY, SPIN-CHARGE AND SPIN-HIGGS COUPLING
The spin susceptibility for cuprates is given by the following expression [8,9]: with the bare spin susceptibility and g(q) = g 0 [1 − 0.1(cos q x + cos q y )]. However, we now show that the coupling of the spin channel (and therefore the spin susceptibility or the resonance peak) to either the amplitude of the superconductor or to the charge/Raman vertex is not possible in light experiments, i.e. : in the limit of vanishing transferred momentum, q → 0. Indeed, within this limit, we first calculate the susceptiblity which couples the spin channel (Ψ † k τ 0 Ψ k ) to the Higgs mode (τ 1 ): Analogously, for the coupling between the spin and the charge channel we get As shown by Venturini et al. in [10], it is possible to obtain some contribution of the spin resonance peak even in the Raman response only including a two-magnon process. In this case, however, it has been shown that the spin response is not affecting significantly the intensity of the A 1g Raman response, but it can rather change the peak position.

VI. EXPERIMENTAL DATA
The Raman spectra were measured in high vacuum (p < 10 −6 mbar) using standard equipment: they represent pure symmetry projections which were derived from the raw data by linear combinations [11]. The spectra are measured with reduced resolution to avoid the accumulation of surface layers at low temperature. The raw data for two different temperatures, T = 30 K < T c and T = 120 K > T c are shown in Fig. 4. The spectra contain contributions from electronic and phononic excitations. Due to the large step width and the poor resolution of approximately 20 cm −1 the phonon lines are not well resolved. Only the A 1g lines at 450 and 650 cm −1 are originating from phonons. In general, the phonons may be removed by subtracting the normal-state spectra from those in the superconducting state. Here, we put the plot difference R T =30K − R T =120K . Since the phonons depend only weakly on temperature no anomalies are seen in the spectra of Fig. 5. The negative intensities in Fig. 5 result from the suppression of the scattering inside the gap opening below the critical temperatureT c = 94 K. They are suppressed in the main text.

VII. NUMERICAL RESULTS AND PARAMETERS USED
We have first calculated the Raman response with the full Brillouin zone summation of the k-points, using the parameters ∆ max = 31 meV, broadening δ = 4 meV, t = 125 meV, B = 0.65, µ = −1.2 in the band structure dispersion in Eq. (8): the result is shown in Fig. 6. In particular, considering also the effect of the A 1g phonon peaked at ω ph = 80 meV ( [12]) on the Raman response, the full A 1g Raman response (which also includes quasiparticles and many-body Higgs). The bare phonon contribution is negligible, but when the many-body Higgs interactions are included, it contributes to increase the Raman peak and to shift it towards ω = 350 cm −1 , as well as to screen the response at higher frequencies revealing a shoulder at ω = 650 cm −1 . These features have to be included to reach a higher qualitative agreement with the experimental data. For the approach of the integration over the Fermi surface, as derived in Section III and whose results are shown in Fig.2 in the main text, we used the following values for the parameters which appear in Section III: ∆ max = 29 meV, broadening δ = 4 meV, α 1 = 136 meV, α 1 /β 1 = 1.09 and α 2 = −0.15 α 1 . They are all reasonably within the range used in the previous literature [5,13]. The parameters chosen for the results in the main text correspond to those which provide the best fit of all the three symmetries A 1g , B 1g , B 2g symultaneously. In fact,we notice that the results in peak position and intensity of each symmetry Raman response is not independent from the others for changing the parameters involved. Therefore, we also show systematic calculations for different values of the parameters in Fig. 7, where in each panel the parameter which differs from the optimal one is explicitly written. We notice that the actual intensities of the peak are within 10 − 15% with respect to the best fit, therefore showing a minor sensitivity to the parameters when changed within a reasonable interval. In addition, we still emphasize that changing a parameter affects the three symmetries simultaneously, and the best fit is not chosen by fine tuning of the A 1g symmetry, but on the most reasonable accordance of all the three symmetries at the same time, obtained with the same set of parameters. This is consistent with the robustness and the unambiguous presence of the Higgs oscillations' contribution.

VIII. TRANSIENT REFLECTIVITY CHANGE
For linear processes, the time-dependent response r(t) of a system which undergoes a time-dependent perturbation p(t ) is given by where f (t−t ) is the response function. In our case the response function is the real part of the first-order susceptibility, For nonlinear processes like those occurring in terahertz pump-optical probe (TPOP) experiments, we need to use higher order susceptibilities: thus, the nonlinear polarization is given by where the response function χ (3) (t − t ) is the third-order susceptibility, P (t) is the response of the system (induced nonlinear polarization), while E pump and E probe are the pump and probe electrical fields, respectively. The timedependent response function can be calculated with the inverse Fourier transform The transient reflectivity change due to a pump electric pulse E pump (t) is given by where χ (3) (t−t ) is the real part of the third order susceptibility in real time, t being the time at which the reflectivity change is detected, t < t the previous time at which the pump pulse interacts with the condensate.

Calculations and results
Recently, Katsumi et al. [14,15] demonstrated that the relative intensity between the transient reflectivity change projected onto the symmetries A 1g and the B 1g in a TPOP experiment on Bi2212 cannot be explained only with the QP contribution, thus deducing that the Higgs should play an important role. Here we considered an adiabatic electric pulse (in analogy to [14]) with the form E pump (t) = sin (2πω p t + ϕ) e −(t−t0)/(2σ 2 ) , with frequency ω p = 0.53 THz, ϕ = 1.65 rad, t 0 = 2.55 ps, σ = 0.63 ps. Then, we used Eq. (81) to evaluate the time-dependent transient reflectivity change on a d-wave superconductor, inserting the third-order susceptibility for a fixed symmetry containing the full many-body Higgs response (in analogy with the Raman calculation of the previous sections). The result, plotted in Fig. 8 a, clearly explains the maximum value of ∆R/R in the A 1g symmetry higher than the B 1g , despite the hypothesis raised in Katsumi's paper of the necessity to consider also strong correlation effects, such as the retarded phonon-mediated interaction in the s-wave case. Moreover, we predict the appearance of characteristic oscillations when a shorter THz pump (ω p = 5.57 THz, σ = 0.022 ps, see the inset of Fig. 8b) is used, simulating a quench experiment [6]. As shown in Fig. 8b, including the full many-body Higgs contribution to the third-order susceptibility, in this non-adiabatic regime oscillations with a period T ∼ 0.08 ps appear in the decaying region of ∆R(t)/R, while a lower intensity and no oscillations are present considering only the quasiparticles and charge fluctuations contributions to the A 1g response.