Theory of symmetry-resolved quench-drive spectroscopy: Nonlinear response of phase-fluctuating superconductors

Recent experiments on cuprates have shown the possibility of opening a gap above the superconducting critical temperature, in the so-called phase-fluctuating state, by enhancing the phase coherence of preformed Cooper pairs. Quench-drive spectroscopy, an implementation of 2D coherent spectroscopy, has emerged as a powerful tool for investigating out-of-equilibrium superconductors and their collective modes. In this work, we enrich the quench-drive scheme by developing a systematic generalization to study the nonlinear response of $d$-wave fully incoherent Cooper pairs in a symmetry resolved manner. In particular, we do not only show that it is possible to obtain a third harmonic signal from fully incoherent pairs with an equilibrium vanishing order parameter, but we also characterize the full flourishing 2D spectrum of the generated nonlinear response. The results provide a deeper theoretical insight on recent experimental results, opening the door to a new symmetry-driven design of future experiments on unconventional and enhanced superconductors.


I. INTRODUCTION
Since their discovery, high-temperature superconductors have been intensely studied because of their properties and rich phase diagram [1][2][3].These unconventional superconductors are characterized by a complex order parameter whose value depends on the quasiparticles' crystal momentum: it can assume both positive and negative values with maximum absolute value at the antinodal points of the Brillouin zone, while vanishing at the nodal points [4].This character is a result of the B 1g symmetry of the d x 2 −y 2 superconducting pairing, descending from their D 4h crystal structure [5,6].However, some features of this class of materials are still under debate, such as the conditions and possibility to induce and experimentally observe collective modes [7][8][9][10][11][12][13][14], or the origin of the pseudogap phase [15][16][17].
In particular, various attempts have been made to study and detect the amplitude Higgs mode even in unconventional superconductors, both investigating the nonequilibrium nonlinear behavior of these materials [12,18], and characterizing the symmetries of their response [19,20].Recent advances have shown that the nonlinear behavior of unconventional superconductors when probed by light emerges as the blending of different contributions, depending on electron-hole doping and impurity concentration, among others [21][22][23].
Beside this, it has been suggested that the pseudogap phase is a precursor of the superconducting state, characterized by finite pairing strength and pre-formed Cooper pairs with phase incoherence [15,16,24].Even if this picture is controversial and has been disproved to some extent, in cuprates, in a region of the phase diagram above the superconducting critical temperature, the superconducting phase is incoherent [25][26][27].When an electromagnetic field interacts with a phase- * matteo.puviani@mpl.mpg.defluctuating superconductor, for some values of intensity and frequency of the incident radiation, it is possible to induce phase coherence among pre-existing phaseincoherent Cooper pairs: this process is responsible for the transient enhancement of the order parameter, or the appearance of a finite superconducting gap in the case of complete phase-incoherence [24,28].If this transition from an incoherent to a partially-coherent phase is fast enough, such as when induced by a short-time quench pulse, then oscillations of the order parameter (quasiparticles' and amplitude mode's excitations) can be produced as well, similarly to what happens in light-induced superconductors [29].
Moreover, in the last years it has been theoretically shown and experimentally observed the generation of odd higher harmonics from driven superconductors: this result originates from the nonlinear behavior of the optical kernel in the superconducting state [8,18].2D coherent spectroscopy (2DCS) on superconductors [30] has developed as a systematic generalization of pump-probe [31] in the context of the broader concept of high-dimensional spectroscopy [32][33][34][35][36].In addition, quench-drive spectroscopy (Fig. 1(a)) has been proposed by Puviani et al. [37,38] as a specific scenario of THz 2DCS to study superconductors, combining a few-cycle short-time quench pulse and a long-time multi-cycle driving field.This allows to obtain a complex 2D nonlinear response embedding many nonlinear contributions, providing useful information on the optical kernel of the superconductor.In fact, this technique can be used on superconductors to study high-harmonic generation, and to address quasiparticles' excitations as well as collective states, as shown by a recent experimental realization [39].Since then, these two dimensional spectroscopies have widely developed [40][41][42][43], proving to be suitable for extracting novel information and details on the superconducting order parameter and its collective modes [23,44].
In this work we combine the quench-drive spectroscopy technique, which allows to investigate the non- equilibrium behavior of materials, with the symmetry selection allowed by pulses' polarization typically used in other spectroscopic techniques, such as Raman [5,45,46] or birefringence [28] spectroscopy.Here, we study the nonlinear response to quench-drive pulses of fully phaseincoherent Cooper pairs with d x 2 −y 2 -wave pairing symmetry (Fig. 1(b)-(c)), as in unconventional superconductors.In particular, we systematically investigate the order parameter's dynamics and the high-harmonic generation process as a function of the real time and the quench-drive delay time, as well as their Fourier spectra.Interestingly, we find that it is possible to observe a B 1g third harmonic signal in the nonlinear current, mediated by the transient increase of phase coherence, depending on the polarization and frequency of the quench and drive pulses.We also confirm previous results obtained only with a quench pulse [28], showing that the B 2g symmetry is the main responsible for the increase of phase coherence.Furthermore, we analyze the full rich 2D Fourier spectrum of the generated current, identifying mixed quench-drive signals with different intensities in different part of the spectrum, classifying them according to their symmetry.The paper is organized as follows: in Section II we provide a brief theoretical overview of nonlinear current generation by incoherent Cooper pairs with d x 2 −y 2 symmetry.In Section III we describe the foundations of symmetry-resolved quench-drive spectroscopy.In Section IV we show and analyze the numerical results, studying the symmetry-resolved nonlinear response obtained for different configurations of quench and drive pulses.Eventually, in Section V we conclude summarizing the work and providing an outlook for possible extensions and future research.In the Appendices A and B we propose the full theoretical calculations of the pseudospin model and the quench-drive nonlinear response generation, while in Appendix C we provide more results, obtained with a different choice of the quench and drive frequencies.

II. NONLINEAR RESPONSE OF INCOHERENT COOPER PAIRS
In this section we theoretically investigate the current generated by a clean high-temperature superconductor subject to quench and drive pulses: the solution is obtained by solving the Bloch equations derived from the pseudospin model of the BCS Hamiltonian described in Appendix A and B. We want to model a state with phasefluctuating superconductivity characterized by the presence of pre-formed incoherent pairs: for this purpose, as in Ref. [28], we use a new artificial equilibrium superconducting state obtained by adding a random momentumdependent phase ϕ k to the original Cooper pairs' state.Therefore, the strength of the pairing potential remains unchanged, as well as the number of total Cooper pairs, while the superconducting order parameter decreases due to the reduced coherence.According to the maximum angle ϕ max which defines the range of the random phase ϕ k , with ϕ k ∈ [−ϕ max , +ϕ max ], we are able to describe different conditions of the material, from the pure superconducting phase for ϕ max = 0, to the complete loss of coherence for ϕ max = π.We define the gap of the pure superconducting state ∆ (0) k = ∆ (0) 0 f k , and the superconducting order parameter in the presence of incoherent pairs as where V is the same pairing strength of the original state, and The superconducting gap in the new equilibrium state can be written in the pseudospin formalism as [47] ∆(ϕ) where we have introduced the equilibrium pseudospin components In order to describe the dynamics of the system, we use the Heisenberg's equation of motion with the new pseudomagnetic field defined as In the presence of an external gauge field represented by the vector potential A(t) coupling to the electrons, the pseudospin changes in time according to The external electromagnetic field is included in the pseudo-magnetic field by means of the Peierls' minimal substitution k → k − eA(t) in the fermionic energy, resulting in The equation of motion in Eq. ( 4) can be decomposed into a set of differential equations, whose solution provides the time-dependent value of the pseudospin σk (t).Once this term is known, we can obtain the value of the time-dependent order parameter ∆(ϕ) (t) = ∆(ϕ) (0) + δ ∆(ϕ) (t), as well as the generated nonlinear current (see Appendix B).However, we notice that the complex order parameter can be written as where θ is the global phase of the superconducting gap, which differs from the local phase of the Cooper pairs in momentum space, ϕ k .As a consequence, the gap equation is not self-consistent anymore (see Eq. ( 1)) and the value of the gap is subject to some time-dependent noise due to the phase incoherence of the preformed pairs.In the full generated current, we can distinguish two non-vanishing contributions: namely, a linear component with the same oscillating behavior of the driving field and a non-linear term including all higher orders Since the third pseudospin component in equilibrium is independent on the phase coherence (Eq.( 3)), the linear current in Eq. ( 9) is always nonzero, even for fully incoherent Cooper pairs and vanishing gap.More details on the solution of the equation of motion and the derivation of the generated current for the quench-drive setup are provided in Appendices A and B.

III. SYMMETRY-RESOLVED NONLINEAR 2D SPECTROSCOPY
In this section we propose the theoretical foundations for the symmetry-resolved quench-drive spectroscopy, identifying the main nonlinear components for different configurations and their corresponding symmetry [48].First, we can conveniently write the frequency spectrum of gap oscillations, as obtained from the solution of the Bloch equations after transforming into Fourier space, using the convolution operation ( * ) as follows [38]: with i, j ∈ {x, y}, and the vector potential A i,j including both the quench and the driving fields.The Raman-like factor where ĵ is the unitary vector along the direction of j) represents the second-order lightmatter coupling and includes the overall symmetry of the gap oscillations (Fig. 2).
ρρ γijγ kl (with {i, j, k, l} ∈ {m, q, d}, which are the measurement, quench and drive axis, respectively) in a 2D quench-drive spectrum on a D 4h crystal.We considered three (I)-(III) given input (quench and drive) directions as well as two (a),(b) measurement axes.
of the D 4h point group as follows [45]: with x′ = (x+ ŷ)/ √ 2 and ŷ′ = (x− ŷ)/ √ 2, corresponding to an angle with respect to the x axis of π/4 and −π/4, respectively.The general rule given the angles α and θ with respect to the x axis reads [28]: Similarly, the third-order nonlinear current, which represents the lowest-order non-vanishing nonlinear contribution, can be written as with i, j, k, l ∈ {x, y}.We notice that in this expression the order parameter's oscillations of Eq.( 11) are embedded into the time-dependent pseudospin component δσ z k (t, ∆t).Considering its spectrum in Fourier space, it can be conveniently re-written as where we introduced the third-order nonlinear susceptibility χ ρρ .Here we omitted the sum over k and the frequency dependencies for convenience of notation.
As an example, we can derive the symmetry of the xyx ′ y response, which enters the nonlinear current along the x direction with interaction of pulses along x′ and ŷ , as follows: with χ (3) which is the only nonvanishing term in Eq.( 16) after summing over the full Brillouin zone.When analyzing the quench-drive spectra, we can substitute the subscripts m, q and d representing the measurement, quench and drive axis, respectively, to the Eq.s ( 11), ( 14) and (15).
As shown by Puviani et al. [38], there are six contributions of the third-order nonlinear susceptibility in a quench-drive spectroscopy setup, which sum up to pro-   vide the full nonlinear response, namely (Fig. 3): Each of them can be decomposed into the D 4h symmetry irreps as shown before.For example, selecting the output along the x axis parallel to the driving field, and with a quench pulse along the xy diagonal, i.e. m = x, q = x ′ , d = x, for the purely driving response we get: where the first Raman factor represents the measurement-driving vertex, while the second corresponds to the driving-driving one.Analogously, for the mixed quench-drive response quadratic in the quench amplitude field we have and A1g + χ (3) which will appear at ν ̸ = 0 in the two-dimensional quench-drive Fourier spectrum of the nonlinear response.
Interestingly, this example shows in practice how the nature of the two-dimensional spectroscopy allows to extract an A 1g symmetry response (and similarly the B 1g and B 2g ) from only one susceptibility component.Moreover, the presence of multiple contributions for different values of the 2D frequency components (ω, ν) allows to measure and selectively address all the symmetries response with only one experiment.

IV. NUMERICAL RESULTS AND DISCUSSION
In this section we present the results obtained from the numerical implementation of the expressions and timedependent Bloch equations show in Sec.II and III.We modeled the electronic band dispersion of the unconventional superconductor as ϵ k = −2t(cos k x + cos k y ) − µ, where the quasimomentum components are expressed in units of the lattice constant a.We used the values of t = 125 meV for the nearest-neighbour hopping energy, chemical potential µ/t = −0.2,obtaining an electron occupation n = 0.9 as in Ref. [24].For the d x 2 −y 2 order parameter ∆ For the pulses we used a few-cycle quench and a gaussianshaped long-duration drive, with amplitudes A q = 0.8 and A d = 0.8, and gaussian envelopes 2σ 2 q = 0.01 ps 2 and 2σ 2 d = 5 ps 2 , respectively.The maximum intensity used for each pulse is provided for the corresponding vector potential in units of ℏ/(e a), where e is the electron charge and a the lattice constant.Moreover, we chose the frequency of quench and drive to be different but both in the THz spectrum, with Ω d = 11 THz and Ω q = 7 THz.In general, different choices of amplitude and frequencies can be made in order to suppress or enhance specific symmetry contributions.In this work, we focused on the fully phase-incoherent Cooper pairs, with ϕ max = π, for which ∆(ϕ) = 0.In our calculations and analysis we restricted ourselves to only three quench-drive symmetry configurations, namely studying the behavior of the superconducting gap (Sec.IV A), and the generated nonlinear current (Sec.IV B) along the (a) x and (b) y direction, as well as their corresponding 2D spectra.More results, obtained with drive and quench frequencies Ω d = 3.66 THz and Ω q = 7.48 THz, respectively, are provided in Appendix C: in this case the quench pulse is nearly resonant with the bare superconducting gap, while the drive is at a much lower energy.

A. Emergent superconducting gap and oscillations
At first, we calculated the behavior of the order parameter, i.e. the superconducting gap, within the quenchdrive spectroscopy setup.In Fig. 4 (I)-(III) (a) we show the 2D time-dependent behavior of the absolute value of the gap.Since the initial system is formed by incoherent pairs, the initial superconducting gap is zero.However, when the quench and drive pulses perturb the incoherent state, they are able to induce coherence in the Cooper pairs, giving rise to a finite gap value, in accordance with Ref. [28].However, thanks to the quench-drive spectroscopic technique, exploiting the symmetry resolution for different quench and drive directions, we can analyze more in depth the gap behavior and the symmetry of its oscillations.Indeed, in Ref. [24] it was shown that a quench pulse along the x axis, i.e. with α q = 0, tends to reduce the superconducting gap decreasing coherence, while with α q = π/4 it is increased.Here, we go beyond that scheme observing that, with the given frequencies of the pulses, a long driving pulse with α d = 0 can also induce coherence in a fully incoherent setup, while a quench along the same direction keeps suppressing it (plots (I)(a) and (III)(a) of Fig. 4).On the other hand, a long driving pulse with α d = π/4 can also increase the gap coherence, but with less efficiency (plot (II)(a)).In order to understand this, we can use the symmetry table in Fig. 2: in fact, for both the (I)(a) and (III)(a) conditions, the pairs are excited mainly in the B 1g symmetry channel.On the contrary, in the (II)(a) scheme the gap is excited with a predominant B 2g symmetry.As a consequence, the B 1g symmetry enhances the gap if used in a driving, while it tends to suppress it if imposed by a short quench.Additional information can be extracted from the analysis of the 2D Fourier spectra of the complex gap, as shown in Fig. 4(b).On the one hand there is a 2Ω d oscillation for the (II) scheme, which results from the B 2g excitation, while no 2∆ peak (originating from quasiparticles' and amplitude mode excitations) appears here.On the other hand, in the schemes (I) and (III), where the B 1g symmetry is mainly excited as the relevant one, we notice dominant frequency components at ω ≈ 0.5 Ω d and ω ≈ 2.5 Ω d .The reason is that within this symmetry the dominant excitation of the superconducting gap is provided by the quasiparticles' excitation and amplitude mode, which have an intrinsic frequency of ω = 2 ∆(ϕ) and ω ≈ 0.4 ∆(ϕ) , as predicted in Ref. [20].

B. Nonlinear current generation
Since the order parameter is not easily accessible in a direct way in experiments, we analyze here the generated nonlinear current by the material: this is because the linear current contains a strong response from the incoherent pairs, while the interesting information is contained in the purely nonlinear part.In Fig. 5 we show the 2D current and the corresponding spectra for the (I)-(III), (a)-(b) configurations, indicating the main symmetry contributions to each term, obtained from Fig. 3.We first notice that the current measured along the x direction (sub-plots (a)) follows the behavior of the gap in Fig. 4, even though the intensity peak for the B 1g sym-metries (I),(III) is one order of magnitude larger than the one with B 2g (II).This can be partly ascribed to the pulse duration and the frequency difference between the quench and the drive, even though the corresponding gap intensities are in the opposite order.The calculations performed selecting the polarization along the y axis are particularly interesting: in fact the response of configuration (I)(b) vanishes (in accordance to the symmetry-resolved susceptibility in Fig. 3), and the response in (II)(b) is surprisingly lower than the one in (III)(b), even though the gap for t = ∆t in the latter case is smaller than in the former.We can also notice that the response in (II)(b) occurs only when quench and drive overlap and extends along the t axis, while the current in (III)(b) is visible only along the diagonal t = ∆t, starting when the driving overlaps with the quench.This means that in the former case the χ mdqq , χ mddq and χ mqqd are the most relevant contributions, while in the latter χ mqdd is, with B 2g and B 1g dominant symmetry, respectively.Overall, the B 2g symmetry is responsible for the gap enhancement from a short pulse, while the B 1g symmetry dominates when a long driving is applied, as well as in the nonlinear current generation.Additional information can be extracted from the 2D spectra, obtained with the Fourier transform of the timedependent plots (Fig. 6).In general, the signals at ν = 0 are independent of the quench pulse, while all the diagonal lines originate from at least a quench pulse component.The horizontal lines with ν = const., which appear in I(a) and III(a) in correspondence of the first harmonic signal, are also independent on the ω frequency and are generated by χ (3) mddd .We first notice that, while in the (II)(b) scheme the most prominent features are peaks at (ω = Ω d , ν = nΩ d ) followed by diagonal spectral lines, in (III)(b) the diagonal features peaked below ω = Ω d are more visible, at ω ≈ Ω q , ν ≈ Ω q .In particular, the diagonal signal starting from the origin and with ω = −ν is the sum of the contributions of nonlinear susceptibilities χ The third harmonic generated by the driving pulse, appearing along the vertical axis for ν = 0, is generated by the third-order nonlinear susceptibility χ (3) mddd , and appears in Fig. 6 I(a) and III(a).Its importance is twofold: firstly, this generally proves that it is possible to generate a third harmonic response even in fully incoherent Cooper pairs exhibiting an initially null gap, when properly quenched and driven.This feature has been experimentally shown in cuprate superconductors above their critical temperature, where a phase-fluctuating phase with vanishing gap is expected [10].Secondly, the third harmonic is generated only when the B 1g symmetry is explicitly excited (see also Fig. 3).However, we can also notice that in configuration (III)(b) there is a non-vanishing third-harmonic component at ν = 0, originating from a diagonal line which accidentally overlaps with ν = 0 due to a higher-order quench-drive mixing, of the kind χ (5) mq,qddd .

V. CONCLUSION
In this work we have calculated the nonlinear response of a phase-fluctuating superconductor with d x 2 −y 2 pairing symmetry without phase coherence, characterized by vanishing superconducting gap in equilibrium.We have adopted the recently proposed quench-drive spectroscopy scheme [37,38], with THz pulses, inducing a finite superconducting gap and analyzing the generated nonlinear current response.In particular, we have developed a symmetry-resolved analysis, which allows to selectively address symmetry components according to the quench and drive pulses and the measurement axis chosen.Our results confirm previous findings on the possibility to enhance the phase coherence with quench pulses, and show additional unique features such as third-harmonic generation and high-harmonic modulations for specific quench and drive polarizations and symmetry measurements.In the future, symmetry-resolved quench-drive spectroscopy can be applied in different scenarios, ranging from the study of novel superconductors, exploring the phase diagram, to the characterization of the nonlinear response and collective modes in superconductors.
with δσ k (t) = (x k (t), y k (t), z k (t)).The external electromagnetic field is included in the pseudo-magnetic field by means of the minimal substitution k → k − eA(t) in the fermionic energy, resulting in The Heisenberg equation of motion for the pseudospin can be written in the Bloch form, Here, for simplicity of calculations and without loss of generality, we assumed a real order parameter, ∆ ′′ (t = 0) = 0, at the initial time t = 0, so that y(0) = 0.
Appendix B: Quench-drive nonlinear response of a superconductor In order to describe a quench-drive experiment we have to choose the appropriate vector potential A(t) = A q (t) + A d (t) = A q (t − t q ) + A d (t − t d ), where A q (t) is the quench pulse centered at time t = t q , A d (t) is the driving field centered at t = t d .Introducing the timedelay ∆t = t d − t q and putting t d = 0 we can rewrite A(t) = A q (t + ∆t) + A d (t).Therefore the expressions in Eq. (A8) depend on both t and ∆t.The solution of Eq. (A8) provides the time-dependent pseudospin, from which the time-dependent order parameter ∆(t) and the generated current j(t) can be calculated.From the self-consistent gap equation we get The current generated by the superconductor in this quench-drive setup is given by the expression In order to separate the linear and the nonlinear contributions to the full generated output current, we first expand the velocity in series of powers of the vector potential A: Here we omitted the explicit time-dependence of A and κ from t and ∆t.Now we can rewrite Eq. (B3) as In particular, the equivalence ∇ κ v κ A=0 = ∇ k v k holds.Therefore we can simplify Eq. (B3) writing Additionally, we expand the electron number where we used the relation nk = 2σ z k + 1.Therefore, Eq. (B2) can be expanded in the lowest orders as We can decompose the generated current along a generic x ′ axis, in order to extract specific symmetry components: The contribution to the current in Eq. (B7) at the lowest order in the external field is given by which vanishes due to parity.At the next order, the linear term reads The first non-vanishing nonlinear term generated by the driving pulse is the third order component where z k (t, ∆t) is the third component of the pseudospin vector σ k (t, ∆t), containing the information of the state of the system perturbed by the quench pulse.The full nonlinear response, which is given by the sum of all the odd orders of the current expansion, can be conveniently calculated by  In general it is useful to extract the 2D frequency spectrum of such a response, in order to analyze the relevant high harmonics: for this reason, we compute the 2D Fourier transform with respect to the evolution time t and the quench-drive delay time ∆t, obtaining the reciprocal variables ω ≡ F(t) and ν ≡ F(∆t), respectively.
As an example, the 2D Fourier transform of Eq.(B11) to provide the third harmonic response of the driving frequency along the direction x′ is j x where F x ′ is an appropriate function independent of the frequency which contains information on the driving shape, the measurement axis and the quasiparticles' momentum.

Appendix C: Additional results
We present here additional results, obtained for the same quench and drive intensities as the ones in the main text, i.e.A d = A q = 0.8, as well as pulses' duration and shape, but with different frequencies: namely, Ω d = 3.66 THz and Ω q = 7.48 THz, respectively.As a consequence, the quench is nearly resonant with the maximum superconducting equilibrium gap, ∆ max = 31 meV = 7.5 THz, while the driving pulse is far from it.The generated nonlinear current is therefore affected by these conditions, and the response appears in some cases qualitatively and quantitatively different from the one obtained in the main text, even if the symmetries involved in the quench-drive spectra are the same.
We first analyze the behavior of the absolute value of the superconducting gap as a function of the real time t and the delay time ∆t.Interestingly, we realize that the gap is very poorly excited in configurations (I) and (III) due to the B 1g symmetry and the driving contribution, with a maximum amplitude of about 3 meV, and as low as 1 meV on the central peak of the driving field at t = 0. On the other hand, the scheme (II) has a higher gap excitation.The corresponding 2D Fourier spectra show that, for schemes (I) and (III) there are no proper gap oscillations, but rather an almost frequency-independent enhancement, plus quench-induced contributions (vertical lines in Fig. 7(I),(III)(b)).On the other hand, for the scheme (II) with diagonal quench and drive pulses, where the B 2g symmetry is excited, a gap oscillation at ω = 2Ω d due to the driving appears, as well as at ω = 2∆, which includes Higgs and quasiparticles' excitations at twice the induced gap amplitude, around 8.5 meV.
We now turn to the generated current: due to the different resonance conditions, we expect the current responses involving a quench pulse to be more intense, saturating the purely drive signals.In particular, for schemes  (I) and (III) where the gap excitation and oscillations are much smaller, we expect the susceptibility term independent of the frequency to be the most relevant [8].In response involves the driving pulse, and the corresponding susceptibility is now more far from resonance.On the other hand, the scheme (II), with quench and drive pulses along the xy diagonal axis, provides now a slightly stronger response, involving mainly the quench pulse.The 2D Fourier spectra in Fig. 9 are even more dense of information.In fact, the spectra of schemes (I) and (III) present much fewer features than with the choices of frequency in the main text: in particular, (I)(a) and (II)(a) have a weaker third harmonic generation and only one diagonal line, representing the nonequilibrium modulation due to the quench pulses.Moreover, the current mea-sured along the y axis in (III)(b) has no first harmonic contribution, and is saturated by the same nonequilibrium modulation of (I),(III) (a).On the other hand, the spectra of (II) are much more complex, exhibiting more and stronger frequency modulations and the emergence of a non-equilibrium third harmonic at ω = 3Ω d , ν = 0.All in all, we have observed how the nonlinear signal is still present in (I) and (III) configurations, the current intensity being higher than in scheme (II), even if the gap is less excited in the former.The reason of this behavior can be ascribed once again to the symmetries involved and here identified.

Figure 1 .
Figure 1.Quench-drive spectroscopy of unconventional superconductors.(a) The 2D quench-drive spectroscopy is performed with a short (quench) pulse, followed by a long gaussian-shaped (drive) pulse at a delayed time ∆t.The output signal is analyzed as a function of the real time t.(b) Band structure (on the left) and d x 2 −y 2 gap symmetry (on the right) of the unconventional superconductors studied in this work.The nodal points are identified with ∆ k = 0, for k = (kx, ky = ±kx), while the antinodal ones are at k = (0, ±π), (±π, 0).

FrequencyFigure 4 .
Figure 4. Gap oscillations and frequency spectra.(a) 2D oscillations in (t, ∆t)of the absolute value of the superconducting gap, |∆|, for the three scheme configurations (I)-(III) described in the main text and illustrated by the plots of quench and drive pulses.The main symmetry contributions are written for the strongest signals, according to the table inFig.2.(b) Absolute value of the 2D Fourier transform of the full complex gap, |F{∆(t, ∆t)}| = |∆(ω, ν)|.
Figure 4. Gap oscillations and frequency spectra.(a) 2D oscillations in (t, ∆t)of the absolute value of the superconducting gap, |∆|, for the three scheme configurations (I)-(III) described in the main text and illustrated by the plots of quench and drive pulses.The main symmetry contributions are written for the strongest signals, according to the table inFig.2.(b) Absolute value of the 2D Fourier transform of the full complex gap, |F{∆(t, ∆t)}| = |∆(ω, ν)|.
max (cos k x − cos k y )/2 we chose the value ∆ (0) max = 31 meV.The calculations were performed with a summation over the full Brillouin zone {k x , k y } ∈ {−π, π} with a homogeneous square sampling and a total number of k points N k = 10 6 .For the timedependent evolution we used a time-step of δt = 3 • 10 −4 ps, and for the quench-drive delay δ∆t = 2.5 • 10 −2 ps.

A 2 AFigure 5 .
Figure 5. 2D nonlinear current.Plots of the generated nonlinear current as a function of real time t and the quench-drive delay time ∆t, for three different schemes (I)-(III) described in the main text, and two polarized output measures along (a) x and (b) y axis, respectively.The symmetries written represent the main contribution, according to the table of the nonlinear susceptibilities in Fig. 3. Be aware of the different color scale for each plot.

Figure 6 .
Figure 6.2D nonlinear current spectra.Plots of the Fourier transform of the nonlinear current in Fig.5, for three different schemes (I)-(III) described in the main text, and two polarized output measures along (a) x and (b) y axis, respectively.Be aware of the different Log color scale for each plot.

Figure 7 .
Figure 7. Gap oscillations and frequency spectra.(a) 2D oscillations in (t, ∆t)of the absolute value of the superconducting gap, |∆|, for the three scheme configurations (I)-(III) described in the main text and illustrated by the plots of quench and drive pulses.(b) Absolute value of the 2D Fourier transform of the full complex gap, |F{∆(t, ∆t)}| = |∆(ω, ν)|.

Figure 8 .
Figure 8. 2D nonlinear current.Plots of the generated nonlinear current as a function of real time t and the quench-drive delay time ∆t, for three different schemes (I)-(III) described in the main text, and two polarized output measures along (a) x and (b) y axis, respectively.This figure corresponds to Fig. 5, here obtained with different frequencies of quench and drive pulses, as explained in the main text.Be aware of the different color scale for each plot and with respect to Fig. 5.

Figure 9 .
Figure 9. 2D nonlinear current spectra.Plots of the Fourier transform of the nonlinear current in Fig.8, for three different schemes (I)-(III) described in the main text, and two polarized output measures along (a) x and (b) y axis, respectively.This figure corresponds to Fig.6, here obtained with different frequencies of quench and drive pulses, as explained in the main text.Be aware of the different Log color scale for each plot.

Fig.s 8
and 9 the measured nonlinear responses in time and the corresponding 2D Fourier spectra are shown, re-spectively.We notice that the current measured along the x axis for schemes (I) and (III), involving mainly the B 1g symmetry, is quantitatively different from the one in Fig. 5.The lower intensity (the scales of Fig.s 5 and 8 are different), in fact, is explained by the fact that the main Table of symmetries of Raman factors, γrs, with r, s ∈ {q, p}, where the labels q and p represent the quench and drive pulses' directions, respectively.
Table of symmetry-resolved nonlinear spectra contributions.The table shows the six symmetry-resolved components of the third-order nonlinear susceptibility χ