Quench-drive spectroscopy and high-harmonic generation in BCS superconductors

In pump-probe spectroscopies, THz pulses are used to quench a system, which is subsequently probed by either a THz or optical pulse. In contrast, third-harmonic generation experiments employ a single multicycle driving pulse and measure the induced third harmonic. In this work, we analyze a spectroscopy setup where both a quench and a drive are applied and two-dimensional spectra as a function of time and quench-drive delay are recorded. We calculate the time evolution of the nonlinear current generated in the superconductor within an Anderson-pseudospin framework and characterize all experimental signatures using a quasiequilibrium approach. We analyze the superconducting response in Fourier space with respect to both the frequencies corresponding to the real time and the quench-drive delay time. In particular, we show the presence of a transient modulation of higher harmonics, induced by a wave mixing process of the drive with the quench pulse, which probes both quasiparticle and collective excitations of the superconducting condensate.


I. INTRODUCTION
The superconducting state of matter is characterized by a zoo of collective modes: These include, among others, Higgs, Leggett, bi-plasmon, and Bardasis-Schrieffer modes [1][2][3][4][5][6][7][8][9][10][11].The study of these modes is currently being established as a new field of collective mode spectroscopy, in the sense that bosonic excitations of the condensate reveal information about the underlying superconducting ground state and symmetry properties of the condensate itself [7,10,12,13].The Higgs mode, for instance, can be used as a spectroscopic tool to distinguish between different gap symmetries of unconventional superconductors [10].
The experimental study of collective superconducting modes poses significant challenges.Due to particle-hole symmetry, they generally cannot couple linearly to electromagnetic fields in the spatially homogeneous limit [3,14].Instead, they are activated in a two-photon Raman-like process [15,16].Thus, the main signature of collective modes consists of a renormalization of the nonlinear susceptibility, which can be probed by nonlinear spectroscopic techniques, such as high-harmonic generation, optical Kerr effect, and nonlinear optical conductivity measurements [17][18][19].
Generally, two distinct approaches exist to excite collective oscillations of a superconductor: The first is to apply a short-duration quench pulse, τ ∆ < 1 (τ being the pulse duration and ∆ the energy gap), to suddenly shrink the superconducting gap and excite the system into an out-of-equilibrium state [10].The second approach uses a single longer pulse, τ ∆ ≫ 1, to drive the material into a quasi-steady excited state [10,18,20].
Higher-dimensional spectroscopy techniques have been used to study the nonlinear response of various materials [21][22][23][24], but they have rarely been applied to superconductors [25].Instead, for superconductors, most time-resolved spectroscopies and theoretical studies have focused only on either short quench pulses or a multicycle driving pulse [5,26], without mixing them simultaneously.In the present work, we study the full evolution of the nonlinear current in the BCS superconducting state subjected to a spectroscopic setup where both a quench and a drive pulse are applied (Fig. 1(a)).We show how the spectroscopic data can be clearly analyzed in 2D Fourier space (ω, ω ∆t ), where the two frequency variables correspond to conjugates of real time t and the pump probe delay ∆t, respectively.While these spectra reduce to the aforementioned pump-probe and THG experiments in certain limits, we argue that quench-drive spectroscopy provides a comprehensive way to experimentally extract information of the nonlinear optical susceptibility and the spectrum of superconducting collective modes.We stress that one of the experimental advantages of the proposed quench-drive spectroscopy is that pulse frequencies need not be continuously scanned across a range of frequencies to probe the system, in contrast to simple driving protocols [27].Instead, to achieve frequency resolution and collective mode resonance, only the time delay between quench and drive pulses needs to be swept.
To solve the equations of motion for a superconductor we employ a pseudospin model, extended to describe the new quench-drive setup.This allows for the simulation of the evolution of the order parameter as well a calculation of the current induced in the superconductor [3,7].In addition, we present a diagrammatic approach to systematically interpret the two-dimensional spectra.
The paper is organized as follows: In Section II we introduce the quench-drive spectroscopy mechanism and explain its features diagrammatically.In Section III we describe the theoretical background.We use the pseudospin model to solve the Heisenberg equation of motion in the presence of an external field for the timedependent order parameter, and then calculate the generated nonlinear current.In Section IV we show the numerical results of the nonlinear current for the quenchdrive setup in a BCS superconductor.The discussion of results is provided in Section V. Finally, we give a summary and outlook on future applications and perspectives of quench-drive spectroscopy in Section VI.

II. QUENCH-DRIVE SPECTROSCOPY
We consider here a clean BCS superconductor without impurity scattering.We are interested in the nonlinear response of the superconductor, which is of third order in the external light field in materials with inversion symmetry.In particular, the quasi-equilibrium third order current is determined by the diamagnetic coupling to light, and its time-dependent expression reads [7,26,28] where (B * C)(x) = dyB(y)C(x − y) denotes a convolution integral.The function χ ρρ is the effective densitydensity response of the operator ρ = kσ For ease of notation, we have assumed that all applied electromagnetic pulses are polarized along the x-direction, i.e., A = A x.In the general case, where A can have arbitrary polarization, the density-density response becomes a tensor whose structure encodes additional information about material properties.A specific case of cross-polarized pulses is discussed in Apdx.C.
Within the BCS approximation, the gauge-invariant response for a single-band model has been computed in various references [3,15,26,28] and a general framework for pump-probe experiments based on a quasiequilibrium effective action formalism has been developed in Refs.[5,26].The density-density response is found to be peaked at the resonance frequency 2∆ of the Higgs mode, where 2∆ is the single-particle superconducting spectral gap.It was pointed out, however, that the resonance peak in the density-density response is the result of both single-particle contributions, stemming from quasiparticle excitations, and collective mode excitations.Importantly, it was shown that quasiparticles generally give the dominating contribution to the 2∆-peak in the clean limit, making observation of the Higgs mode difficult [15].Other collective modes of the condensate, such as Leggett [5,29,30], Bardasis-Schrieffer [9,31], or other relative phase modes [13] do contribute significantly to the density-density response and may even persist below the gap.Additionally, the Higgs modes may achieve a sizable signal in the presence of impurities [30,[32][33][34][35] or due to additional processes [36,37].In the present work, we will not focus on the attribution of spectral weight of the nonlinear response to their various origins and instead discuss spectroscopic measurement of the density-density response χ ρρ as a whole.
In Fourier space Eq.(1) can be expressed as where the δ-function is a manifestation of energy conservation, i.e., the three photon frequencies ω i have to sum up to the frequency ω of the induced current.The susceptibility χ ρρ in Eq. ( 3) enters with a functional dependence on the frequency variables ω 2 + ω 3 .In general, the integration over these variables scrambles the resonance spectrum of χ ρρ and no direct signature of collective modes can be recovered from j(ω).
Two approaches to circumvent this problem exist.First, one may choose A as a multicycle THz pulse of the form with τ d ω d ≫ 1 and ω d ∼ ∆ such that it has a narrow frequency spectrum of width τ −1 d centered around ±ω d .Then, the integration variables are constrained to ω i ≈ ±ω d and the susceptibilities are mostly evaluated at χ ρρ (0) and χ ρρ (±2ω d ) to yield the first or third harmonic, j(±ω d ), j(±3ω d ).To map out the functional dependence of χ ρρ one has to vary the driving frequency ω d .This, however, is not easily achievable experimentally.Instead, most current experiments fix the driving frequency ω d and instead attempt to shift the resonance energies contained in χ ρρ .For a superconducting mode, this is simply achieved by varying the temperature in the window (0, T c ).The clear disadvantages of this method are that (1) knowledge of the temperature dependence of the resonances of χ ρρ is required, (2) only modes above 2ω d are visible, and (3) thermal broadening effects are substantial.
The second approach consists of pump-probe setup.Here, we consider a novel pump-probe setup where in addition to a broadband quench pulse A q , the multi-cycle drive pulse A d is utilized.The quench A q has the same form as Eq. ( 4), but with τ q ≪ 1/∆.Further details of the pulse shapes are given in Appendix B. The two pulses are delayed with respect to each other by ∆t yielding the total field A(t) = A q (t + ∆t) + A d (t) (see Fig. 1 and Fig. 9).In frequency space, this results in a phase shift, where we have introduced the notation δ α,q = 1 for α = q and zero for α = d.
In a nonlinear THz experiment, the current j(t) is electrooptically sampled as a function of t and Fourier transformed numerically to obtain j(ω).Multiple such traces are recorded for varied ∆t to assemble the 2D spectrum j(ω; ∆t).Inserting Eq. ( 5) into Eq.( 3), we obtain where we are summing over all combinations {α i } of quench and drive pulse, A q and A d , respectively.We perform a Fourier transform in the parametric time delay ∆t to obtain We can represent the various terms in the sum over {α i } diagrammatically as depicted in Fig. 2. Here, the current is represented by a red wiggly line, and quench and probe pulses are depicted by blue and black wiggly photon lines, respectively.The density-density susceptibility χ ρρ is represented by a fermionic bubble whose internal frequency we have labeled ω n .
All external photon lines carry frequencies with profiles determined by the experimental bandwidth of the pulses A αi (ω i ), where the directionality is marked by arrows of the photon lines.The drive pulse will constrain the external frequencies to ∼ ±ω d .In fact, we will be assuming a sufficiently narrow-band drive pulse, τ d ≫ 1/∆, such that we can approximate The key advantage of the pump-probe geometry lies in the appearance of the second δ-function in Eq. ( 7) that introduces the experimentally accessible variable ω ∆t .We can see this as follows.When α 1 = d and α 2 = α 3 = q, the δ-function presents the constraint δ(ω ∆t +ω 2 +ω 3 ) and the density-density correlation function is evaluated at χ(ω 2 + ω 3 ) = χ(−ω ∆t ).Thus, χ can be pulled out of the integral in Eq. ( 7) and the measured current is directly proportional to the χ ρρresponse, whose frequency dependence can be mapped out by sweeping ω ∆t .The diagrammatic representation of this process is shown in Fig. 2(a).Making use of approximation (8), it follows from energy conservation that the current is non-zero only along the lines ω = ±ω d −ω ∆t in 2D frequency space (ω, ω ∆t ), where j is given by As similar discussion applies to the process depicted in Fig. 2(b).Here, α 1 = α 3 = d and α 2 = q.The susceptibility is evaluated as χ(±ω d − ω ∆t ) and determined the current along the lines ω = ±2ω d − ω ∆t and ω = −ω ∆t .Explicitly, the current is, Figure 2(c) describes the usual THG process that is independent of the quench.In 2D Fourier space it yields a signal at ω ∆t = 0 at the first and third harmonic frequencies of the drive, ω = ±ω d , ±3ω d , where the densitydensity susceptibility is evaluated at fixed values 0, ±2ω d .The wiggly lines represent photons: the black ones correspond from the driving field A d , the blue lines correspond to the quench pulse Aq, the red line represents the generated current.Solid lines denote fermionic bubbles that represent the nonlinear susceptibility χρρ.Due to the δ-function constraint in Eq. ( 7), all blue quench frequencies have to add up to −ω∆t.Energy conservation demands that all incoming frequencies add up to the frequency ω of the induced current.
The remaining diagrams of Fig. 2 can be separated into two classes.In Fig. 2(d) χ ρρ only depends on ω d and the discussion of the THG case applies.In Fig. 2(e-f) the dependence of χ ρρ on integration variables ω i cannot be removed.As a consequence, the resonance structure of χ ρρ is washed out by integration and one obtains a mostly constant signal for frequencies much smaller than the bandwidth of the quench pulse.
In summary, we have shown that the signal j(ω, ω ∆t ) falls onto diagonal lines ω = ±nω d − ω ∆t with n ∈ {0, 1, 2} in 2D Fourier space.From the 2D spectra, one can extract the density-density response according to where m = {0, 1} and the c i denote the background signal that is mostly constant in the limit of a broadband quench pulse.

III. MICROSCOPICAL DESCRIPTION
Having discussed the phenomenological structure of the nonlinear response in a quench-drive experiment, let us now microscopically investigate the response current of a conventional clean superconductor subject to quench and drive pulses.The solution is obtained solving the Bloch equations derived from the pseudospin model of the BCS Hamiltonian.

A. Equations of motion with quench-drive pulses
We write the BCS Hamiltonian using the pseudospin formalism [3,7,38,39] as with the pseudospin vector σk = 1 2 which is defined in Nambu-Gor'kov space, with spinor Ψk = (ĉ † k,↑ ĉ−k,↓ ) and the Pauli matrices τ = (τ 1 , τ 2 , τ 3 ).The pseudo-magnetic field is defined by the vector where ε k = ξ k − µ, ξ k being the fermionic band dispersion, µ the chemical potential.The superconducting order parameter Here, V is the pairing strength, and f k the form factor of the superconducting order parameter.For s-wave pairing one has f k = 1.
In the presence of an external field represented by the vector potential A(t), the pseudospin changes in time according to with σ k = ⟨σ k ⟩ and δσ k (t) = (x k (t), y k (t), z k (t)).The external electromagnetic field is included in the pseudomagnetic 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 Bloch form, ∂ t σ k = 2b k ×σ k , providing the set of differential equations where δ∆(t) = δ∆ ′ (t) + i δ∆ ′′ (t) is the time-dependent variation of the order parameter induced by the external field, such that ∆(t) = ∆ + δ∆(t).Here, we assumed a real order parameter at the initial time t = 0, so that y(0) = 0.The solution of Eq. ( 18) provides the time-dependent evolution of pseudospins, from which the time-dependent order parameter ∆(t) and the generated current j(t) can be calculated.A detailed derivation of the equations of motion is given in Appendix A.

B. Nonlinear current
The current generated by the superconductor in this quench-drive setup is given by the general expression where the velocity is calculated as , and the charge density is defined as We can expand the velocity as a function of the vector potential A(t, ∆t), and expand the current in powers of the external field.The first non-vanishing term of the nonlinear current j N L (t, ∆t) 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), that contains the information of the state of the system perturbed by the quench pulse.The unit vector r i , i = x, y, represents the two directions along which the output current is measured.

IV. NUMERICAL RESULTS
We now present the results obtained from the numerical implementation of time-dependent Bloch equations described in the previous section, solved by means of a Runge-Kutta-4 algorithm without linearization or further analytical approximations.We used the fermionic band dispersion ε k = −2t(cos k x + cos k y ) at half-filling, setting the lattice constant a = 1.The point µ = 0 is special in the sense that it has perfect particle-hole symmetry as well as a van Hove singularity at the Fermi level.But from the viewpoint of collective modes, this symmetry point does not bear any special significance other than presenting a local minimum of the Higgs contribution compared to single-particle excitations in the nonlinear response [15].
We used the values of t = 125 meV for the nearestneighbor hopping energy, s-wave order parameter ∆ 0 = 15.8 meV (corresponding to a frequency of 3.82 THz), and a summation over the full Brillouin zone with a square sampling and a total number of N k = 10 6 points.For the time-dependent evolution we used a time-step of δt = 3 • 10 −4 ps, and for the quench-drive delay δ∆t = 2.5 • 10 −2 ps.For the quench we used a few cycles pulse with central frequency ω q = 4.77 THz, while for the driving an asymmetric pulse with central fre- quency ω d = 4.3 THz, so that both satisfy the condition ω q(d) < 2∆.Both pulses were considered linearly polarized along the x-direction, while the maximum amplitude of the electric field of quench and drive pulse was taken E q = 10.5 kV/cm and E d = 4.7 kV/cm, respectively, assuming a value for the lattice constant of a = 3 Å.For more details on the pulses used, see Appendix B. In Appendix C we analyze the 2D spectra for the case of cross polarization of the two pulses, while in Appendix D we show additional results that were computed using a symmetric Gaussian driving pulse instead of an asymmetric one.

A. Two-dimensional quench-drive spectroscopy
We solved the time-dependent equation of motion in the quench-drive spectroscopy setup using the pseudospin model and calculated the nonlinear current generated by the condensate as described in the previous section.The result plotted as a function of the real time evolution t and the quench-drive delay time ∆t, which is measured as the interval between the maximal peaks of the envelopes of the two pulses (see Appendix B for more details on the pulse shapes), is shown in Fig. 3.We notice that the diagonal line signifies the arrival of the quench pulse, which overlaps with the drive for times t = ∆t in the range t, ∆t ∈ [0, 2] ps.In this region the response is modulated as a function of the delay time ∆t.Next, we Fourier transform the real time variable t into the frequency ω.In Fig. 4 we show the nonlinear gener- ated current intensity as a function of the quench-drive delay time ∆t, namely |j N L (ω, ∆t)|.We notice that both the first and the third harmonic of the fundamental driving frequency ω d are modulated in delay time ∆t, with maximum intensity in the interval 0 ps ≤ ∆t ≤ 2 ps, which corresponds to the range of interference between the quench and the drive pulses, as shown in Fig. 3. Additionally, the signal intensity does not vanish away from ω d and 3ω d , where we instead observe a striped pattern, with each intensity line tilted towards the central time ∆t = 1 ps.These features can be more readily interpreted by plotting the 2D Fourier transform of the current, i.e., as a function of the frequency ω and the delay-time frequency ω ∆t , respectively, shown in Fig. 5.In particular, we notice the first and third harmonic as strong peaks in the central vertical line, at ω ∆t = 0, which corresponds to the equilibrium response of the driving pulse in absence of the quench.Note that it is sufficient to plot two quadrants of the nonlinear current in 2D frequency space, since it follows from j(t, ∆t) ∈ R that j(ω, ω ∆t ) = j(−ω, −ω ∆t ).The modulations in ∆t appear here as broad diagonal lines in 2D frequency space as expected from Eq. ( 11).These features correspond to a dynamically generated four-wave mixing signal due to both the quench and the drive pulses.

B. High-harmonic generation and transient excitation of the superconductor
To clearly distinguish the shape and position of the peaks observed in the frequency 2D plot, we performed one dimensional cuts of Fig. 5 along various lines.Fig. 6 shows the plot along the vertical line at ω ∆t = 0.This corresponds to an equilibrium high-harmonic generation  due to the driving pulse only.We observe a first harmonic peak at ω d and a third harmonic signal at 3ω d .Measurement of the temperature dependence of the THG peak would correspond to a usual THz THG experiment.The intensity of the fundamental and third harmonics (continuous line) are of the same order, since only the third order is plotted.The total current has a dominant linear first harmonic response (dashed line).In addition to the first and third harmonic, however, we notice the presence of an additional shoulder peaks at a frequency ω = 2∆ + ω d and at ω = 2∆ − ω d .These are the result of the intrinsic Higgs and quasiparticle resonance at 2∆ [40] in conjunction with a wave-mixing process with a driving photon of frequency ±ω d .In Fig. 7 we show the prototypical case of a horizontal cut in Fig. 5 along ω = 3ω d .The peak at ω ∆t /ω d = 0 is the equilibrium third harmonic as visible in Fig. 6, while the smaller peaks at ω ∆t /ω d = −1, −2, −3 stem from the modulation of the third harmonic due to the quench.Additional smaller peaks appear in Fig. 4 as a modulation in the delay time ∆t of the third harmonic, giving rise to the characteristic striped pattern.

V. DISCUSSION
We can understand all features in the spectrum shown in Fig. 5 by considering each of the diagrams in Fig. 2 that represent the induced current expanded to third order in various combinations of powers of A d and A q .The equilibrium THG signal proportional to A 3 d has to be independent of ∆t and therefore falls onto the vertical line ω ∆t = 0 in the 2D spectrum j(ω, ω ∆t ).This is represented by the diagram Fig. 2(c), where only the driving field acts on the condensate, and the current spectrum can be described as a function of real-time evolution as in Fig. 6.Interestingly, we also notice that in addition to the aforementioned fundamental and third harmonic, a shoulder peak of the third harmonic at a frequency ω d + 2∆ appears when the driving pulse is asymmetric and not Gaussian-shaped (see Appendix D for the data with the symmetric envelope).This is the direct consequence of the effective quench induced by the driving, which launches free Higgs oscillations alongside the quasiparticle contribution and enhances the intensity of the nonlinear susceptibility at ω H = 2∆ [3,7,36].In principle, a THG experiment for a single temperature would suffice to identify the collective mode resonance.However, this approach strictly relies on the condition 2∆ ≈ 2ω d and is specific to the asymmetric pulse shape [36] (See also Appendix D).The processes described by diagrams (a,b,d-f) in Fig. 2, which involve at least one photon of the quench, are responsible for the signals along diagonal lines ω = ±nω d − ω ∆t with n ∈ {0, ±1, ±2}.The spectral window in which these lines can be observed are related to the bandwidth of the quench pulse.Here, we differentiate between even and odd n.Odd diagonal lines show a peak at ω ∆t = −2∆ and even diagonals are peaked at ω ∆t = ω d − 2∆.This is expected from Eq. ( 11) since the susceptibility χ ρρ is peaked at 2∆ for the modelled single-band superconductor.From Eq. 11, we would additionally expect a peak at ω ∆t = −ω d − 2∆ for the line at n = −2.However, while the diagonal n = −2 line in principle is present, its spectral weight is negligibly small within the corresponding frequency range.Note that diagonal lines at ±n are related by frequency inversion j(ω, ω ∆t ) = j(−ω, −ω ∆t ).
By inspecting the 2D spectrum in Fig. 5, it is now straightforward to extract the resonances of the nonlinear susceptibility χ ρρ .In our case, we observe the four peaks on the diagonal lines from which we extract the value 2∆.If the superconducting condensate supports additional collective modes that nonlinearly couple in the electromagnetic response, their mode frequencies can be readily extracted as well.

VI. CONCLUSION AND OUTLOOK
In this work we proposed and analyzed from a new perspective a pump-probe spectroscopy setup on conventional clean superconductors with a combination of a single-cylce THz quench pulse and a multi-cycle driving THz probe field.We used a numerical approach based on the Anderson-pseudospin model to solve the equations of motion and to calculate the generated nonlinear current.In addition, we investigated the nonlinear optical processes by means of a diagrammatic approach to interpret and explain the obtained results.
In particular, we showed that, in addition to the usual third harmonic generation measured in driving experiments, new features are obtained in a two-dimensional spectrum of the nonlinear generated current.These features are manifest as diagonal lines in 2D frequency space of the time and pump-probe delay and allow for a direct extraction of resonances in the nonlinear susceptibility.
The susceptibility encodes the intrinsic superconducting response of the quasiparticles and resonances of transient excitation of the Higgs mode.
The advantage of a two-dimensional analysis of quench-drive spectroscopy is manifest in the possibility to scan a wider frequency spectrum at once with fixed parameters of quench and drive pulses, by scanning the quench-drive delay time.In addition, with the present setup the quench pulse allows to push the system out of equilibrium, quenching and shrinking the superconducting gap, allowing the driving pulse to probe different states of the superconductor, resulting in different peak profiles and positions in the 2D frequency spectra.
It is also interesting to examine the possibility to extend the quench-drive spectroscopy framework to the case of cuprates, which exhibit a different symmetry of the order parameter in momentum space, and preformed phase-incoherent Cooper pairs [41], which can reveal more information on the competing orders and their symmetries.
All in all, we believe that this work can pave the way towards coherent time-dependent multi-dimensional spectroscopy on superconductors in the THz regime.A full two-dimensional pump-pump-probe spectroscopy with coherent pulses will be the focus of a future work.Its possibilities range from coherent control of superconductors to the study of competing orders, such as superconductivity, charge-density wave, and bi-plasmon among others [21,42].Other systems where the Higgs response is known to be enhanced, such as cuprates, could be interesting to investigate with this spectroscopic approach, to efficiently study the transient non-equilibrium response of quasiparticles and the Higgs mode and to unveil the features of their rich phase diagram.Here, we repeated the same calculations of the nonlinear current generation in the quench-drive spectroscopy setup as in the main text, but using a symmetric Gaussian envelope for the drive pulse.The results are shown in Figs.12-13 and correspond to those in Figs.5-8 in the main text, where the asymmetric drive was used.In the top panel of Fig. 13, the equilibrium high-harmonic generation does not include the shoulder peak at ω = 2∆ + ω d , ω ∆t = 0, since it was generated by the initial effective quench of the asymmetric driving field.However, all other features of high-harmonic modulation and transient excitation at ω ∆t = 2∆ are still present, since they originate from the wave mixing of the quench and the drive pulses, independently of their shape.

Figure 1 .
Figure 1.Quench-drive spectroscopy.(a) Quench-drive spectroscopy setup.A single-cylce quench pulse and a multicycle drive pulse excite both Higgs mode and quasiparticles, resulting in a third-harmonic generation (THG) and a dynamical modulation of higher harmonics.In addition, the driving pulse effectively quenches the system, launching Higgs oscillations.In this illustration we show the asymmetric driving pulse.(b) Representation of the frequency spectrum of the quench pulse Aq(ω) (orange), centered at ω = ωq, and the sum (SFG) and difference frequency generated (DFG) pulses A 2 (ω) (green), centered at ω = 0 and ω = 2ωq, respectively.The grey vertical line represents the position of the critical value 2∆.In the inset the real-time quench pulse is shown.

Figure 2 .
Figure 2. Diagrammatic representation of the nonlinear processes.The wiggly lines represent photons: the black ones correspond from the driving field A d , the blue lines correspond to the quench pulse Aq, the red line represents the generated current.Solid lines denote fermionic bubbles that represent the nonlinear susceptibility χρρ.Due to the δ-function constraint in Eq. (7), all blue quench frequencies have to add up to −ω∆t.Energy conservation demands that all incoming frequencies add up to the frequency ω of the induced current.

Figure 3 .
Figure 3. Time-delay two-dimensional plot of the generated nonlinear current.2D plot of the nonlinear output current j N L (t, ∆t) generated by the driven superconductor, as a function of real time t and the quench-drive delay time ∆t.The narrow diagonal stripes are generated by the quench, while the vertical ones are the response to the drive pulse.The intersection provides a wave mixing pattern for t, ∆t ∈ [0, 2] ps.Here we used the asymmetric drive pulse, with quench and drive pulse frequencies respectively ωq = 4.77 THz and ω d = 4.3 THz.

Figure 4 .
Figure 4. Frequency-time delay two-dimensional plot of the generated nonlinear current.2D plot of the nonlinear current j N L (ω, ∆t) generated by the driven superconductor, as a function of relative frequency ω/ω d (with driving frequency ω d = 4.3 THz) and the quench-drive delay time ∆t.This corresponds to a Fourier transform in the horizontal direction in Fig. 3.The fundamental (ω = ω d = 4.3 THz) and the third harmonic (ω = 3ω d = 12.9 THz) are both modulated in the delay time ∆t.

Figure 5 .
Figure 5. Two-dimensional Fourier-transformed plot of the nonlinear current.2D plot of the generated nonlinear output current intensity j N L (ω, ω∆t) as a function of the real frequency ω and the quench-drive delay frequency ω∆t.It corresponds to the two-dimensional Fourier transform of the data in Fig.3.The vertical response at ω∆t = 0 corresponds to the quench-free superconducting signal, namely the high-harmonic generation due to the driving field.The diagonal lines, instead, represent the transient modulation of the higher-harmonics due to the quench-drive wave mixing.

Figure 6 .
Figure 6.Driven high-harmonic generation.Plot of the generated nonlinear (continuous line) and total (dashed line) current as a function of frequency, j N L (ω) and jtot(ω), respectively.This plot corresponds to a vertical cut in Fig. 5 along ω for ω∆t = 0.The peaks at ω = ω d = 4.3 THz and ω = 3ω d = 12.9 THz correspond to the fundamental and third harmonic, respectively.The smaller peak at ω = 12 THz corresponds to the transient excitation of Higgs and quasiparticles with ω = ω d + 2∆, due to the asymmetric driving pulse.

Figure 7 .
Figure 7. Transient modulation of higher harmonics.Frequency-resolved spectral weight of the third harmonic current as a function of ω∆t, j N L (ω∆t, ω = 3ω d ).This corresponds to a horizontal cut in the 2D Fourier-transform plot in Fig. 5 along ω∆t at ω = 3 ω d .It is possible to identify a transient modulation at frequencies −ω d , −2ω d and −3ω d , respectively.

Fig. 8
photon of frequency ±ω d .In Fig.7we show the prototypical case of a horizontal cut in Fig.5along ω = 3ω d .The peak at ω ∆t /ω d = 0 is the equilibrium third harmonic as visible in Fig.6, while the smaller peaks at ω ∆t /ω d = −1, −2, −3 stem from the modulation of the third harmonic due to the quench.Additional smaller peaks appear in Fig.4as a modulation in the delay time ∆t of the third harmonic, giving rise to the characteristic striped pattern.Fig. 8 corresponds to a diagonal line in the 2D plot of Fig. 5 passing through the point ω ∆t = 0, ω = ω d , and projected along the ω ∆t axis.The peak at ω ∆t = 0 is the signal of the first harmonic.Of particular interest is the peak at ω D = −2∆, which is a direct consequence of the quasiparticle resonance at ±2∆, represented by the process in Fig. 2(d).Due to the wave mixing of the quench and the drive, we have here isolated the intrinsic superconducting response with the characteristic frequency of 2∆.Moreover, the peaks in Fig. 5 along the diagonals placed at ω ∆t = −2∆ + ω d are resulting from the process represented by the diagram in Fig. 2(b), and they disappear when quench and drive have perpendicular polarization, since the corresponding interaction vertex vanishes (see Appendix C).

Figure 9 .
Figure 9. Quench and drive pulses.Plot of the amplitude of the vector potential A corresponding to the quench (top left, frequency ωq = 4.77 THz) and drive pulse, respectively, used in the calculations: the asymmetric drive pulse (center left) and Gaussian shaped drive pulse (bottom left) have a frequency ω d = 4.3 THz.On the right their Fourier spectrum in frequency is shown: the asymmetric drive has a sharp peak in frequency with a slower decay, while the Gaussian-shaped drive is narrower in frequency.

Figure 11 .
Figure 11.Nonlinear current along the x axis: traces from the 2D Fourier transform.Left: Nonlinear current as a function of the frequency ω, obtained with a trace along the vertical axis from the left plot in Fig.10.Right: Nonlinear current as a function of ω∆t with the constraint ω = ω d − ω∆t, obtained from a diagonal cut passing by the fundamental harmonic of the left plot in Fig. 10.

Figure 12 .
Figure 12.Two-dimensional Fourier-transformed plot of the nonlinear current.2D plot of the generated nonlinear output current intensity j N L (ω, ω∆t) as a function of the real frequency ω and the quench-drive delay frequency ω∆t.It corresponds to the two-dimensional Fourier transform of the data in Fig.3.The vertical response at ω∆t = 0 corresponds to the quench-free superconducting signal, namely the high-harmonic generation due to the driving field.The diagonal lines, instead, represent the transient modulation of the higher-harmonics due to the quench-drive wave mixing.The symmetric Gaussian shaped driving pulse was used here.

Figure 13 .
Figure 13.High-harmonic generation and transient modulation with Gaussian drive.Top: nonlinear current j N L (ω) at ω∆t = 0.In contrast to Fig.6, no peak at ω = 2∆ + ω d is present in this case.Center: nonlinear current modulation as a function of the quench-drive frequency ω∆t, j N L (ω∆t, ω = 3ω d ).This corresponds to Fig.7using the Gaussian envelope for the drive pulse.Bottom: j N L (ω∆t) obtained with the condition ω = −ω∆t + ω d , equivalent to Fig.8.