Non-linear Terahertz Driving of Plasma Waves in Layered Cuprates

The hallmark of superconductivity is the rigidity of the quantum-mechanical phase of electrons, responsible for superfluid behavior and Meissner effect. The strength of the phase stiffness is set by the Josephson coupling, which is strongly anisotropic in layered superconducting cuprates. So far, THz light pulses have been efficiently used to achieve non-linear control of the out-of-plane Josephson plasma mode, whose frequency scale lies in the THz range. However, the high-energy in-plane plasma mode has been assumed to be insensitive to THz pumping. Here, we show that THz driving of both low-frequency and high-frequency plasma waves is possible via a general two-plasmon excitation mechanism. The anisotropy of the Josephson couplings leads to marked differences in the thermal effects among the out-of-plane and in-plane response, consistently with the experiments. Our results link the observed survival of the in-plane THz non-linear driving above $T_c$ to enhanced fluctuating effects in the phase stiffness in cuprates, paving the way to THz impulsive control of phase rigidity in unconventional superconductors.

Order and rigidity are the essential ingredients of any phase transition. In a superconductor the order is connected to the amplitude of the complex order parameter, related to the opening of a gap ∆ in the single-particle excitation spectrum. The rigidity manifests instead in the quantum-mechanical phase of the electronic wave function, associated with the phase of the order parameter 1 . Twisting the phase is equivalent to an elastic deformation in a solid, meaning that its energetic cost is vanishing for sufficiently slow spatial variations. On the other hand, since phase fluctuations come along with charge fluctuations, long-range Coulomb forces push the energetic cost of a phase gradient to the plasma energy ω J 1,2 . While for ordinary superconductors this energy scale is far above the THz range, in layered cuprates the weak Josephson coupling among neighboring layers [3][4][5] pushes down the frequency of the inter-layer Josephson plasma mode (JPM) to the THz range. 6,7 The possibility to manipulate the inter-layer JPM by intense THz pulses has been theoretically discussed long ago within the context of the non-linear equation of motion [6][7][8][9][10][11] . This approached turned out to successfully capture the main features of a series of recent experiments 12,13 , even though a full quantum treatment of the JPM able to capture thermal effects across T c is still lacking. On the other hand, non-linear effects induced by strong THz pulses polarized in the planes [14][15][16] have been discussed so far only within the context of the SC amplitude (Higgs) mode, whose excitation energy ω H = 2∆, ranging from 5 to 10 THz in cuprates, appears as a better candidate than high-energy in-plane plasma waves. Nonetheless, the observed monotonic temperature dependence of the nonlinear response 15,16 , its persistence above T c 16 and its polarization dependence 14 do not easily match the expectations for the Higgs mode. The same problem holds considering lattice-modulated charge fluctuations, which are expected to dominate in the clean limit 17,18 but become less relevant [19][20][21][22] and strongly isotropic 22 when disorder is considered.
Here we provide a complete theoretical description of Figure 1. (a) Schematic view of the mexican-hat potential for the free energy F (ψ), with ψ the complex order parameter of a superconductor below Tc. A phase-gradient excitation corresponds to a shift along the minima, while a Higgs excitation moves the system away from the minimun. An intense light pulse with almost zero momentum can excite simultaneously two plasma waves with opposite momenta (in red) or a single Higgs fluctuation (in blue). (b)-(c) Feynman-diagrams representation of the (b) plasma-waves or (c) Higgs contribution to the non-linear optical response. Here wavy lines represent the e.m. field, solid/dashed lines the plasmon/Higgs field, respectively.
the JPM contribution to the non-linear response of layered cuprate superconductors, focusing both on thirdharmonic generation (THG) and pump-probe protocols. We first show that the basic mechanism behind nonlinear photonic of Josephson plasma waves is intrinsically different from the one of the Higgs mode, see Fig.  1. By pursuing the analogy with lattice vibrations in a solid, the Higgs mode is like a Raman-active optical phonon mode. It has a finite frequency at zero momentum, and its symmetry allows for a finite quadratic coupling to light [17][18][19][20][21][22][23][24][25][26] . The phase mode behaves instead like an acoustic phonon mode, pushed to the plasma energy by Coulomb interaction, carrying out a finite momentum at nonzero frequency. As such, zero-momentum light pulses can only excite simultaneously two JPMs with op-posite momenta, making this process strongly dependent on the thermal probability to populate excited states. This feature differentiates drastically the temperature dependence of the THG associated with out-of-plane or inplane JPMs, since the frequency scale of the former is comparable to T c , while it is much larger for the latter. In addition, in contrast to the Higgs mode 17,21 , for a light pulse polarized in the planes the signal coming from JPMs is in general anisotropic, since the momenta carried out by the two plasmons can be along different crystallographic axes. All these features not only contribute to the understanding of the existing experimental measurements, 12-16 but they also offer a perspective to design future experiments aimed at selectively tune non-linear photonic of Josephson plasma waves in layered cuprates. Let us first focus on the out-of-plane JPM. We take a layered model with planes stacked along z. In the SC state the Josephson coupling J ⊥ of the SC phase φ n between neighboring planes sets an effective XY model: An electric field polarized along z enters the Hamiltonian via the minimal-coupling substitution 1 θ n → θ n − (2π/Φ 0 )dA z , with θ n = φ n − φ n+1 , d interlayer distance and Φ 0 = hc/(2e). The corresponding out-of-plane current density I z = −∂H/∂(cA z ) is given by: where J c = 2eJ ⊥ / S, with S surface of each plane. The Josephson current (2) naturally admits an expansion in powers of A z to all orders: where the explicit time convolution of Eq. (3) has been omitted for compactness. Here, following the same approach used so far to investigate the Higgs response 17,18,23,24 , we rely on a quasi-equilibrium description, where the leading effect of the intense THz pump field is to trigger a third-order χ (3) response mediated by plasma waves. The quantum generalization of the model (1) has been widely discussed within several contexts 6,9,10,27,28 . Here we follow the approach of Ref.s 27,28 where long-range Coulomb interactions are introduced within a layered model appropriate for cuprates (see Ref. 29 ). The Gaussian quantum action for the phase mode at long wavelength has the usual form: where ω 2 J = c 2 /ελ 2 c = 8πedJ c / ε is the energy scale of the out-of-plane JPM, iω m = 2πmT are Matsubara frequencies and v = 2 ε/(16πe 2 ), with ε the background dielectric constant. In the classical limit only ω m = 0 is relevant and one recovers the leading term of Eq. (1), i.e. a discrete phase gradient along z, as expected for the Goldstone mode.
To compute the third-order contribution in Eq. (3) we need to derive the effective action S (4) A for the gauge field up to terms of order O(A 4 z ) (see Ref. 29 ). By coupling the gauge field A z to the phase mode via the minimalcoupling substitution in Eq. (2) and by expanding the cosine term, one finds that: where dots denote additional terms not relevant for the χ (3) response. The second term in Eq. (5) can be treated as a perturbation with respect to S G 29 , so that integrating out the JPM one obtains: where I 0 is an overall constant. The vanishing of the denominator in Eq. (7) identifies the resonance of the non-linear kernel. Since the physical mechanism behind the THG is the excitation of two plasma waves, the largest I T HG in Eq. (8) occurs when the pump frequency matches the plasma frequency, i.e. ω = ω J . This has to be contrasted e.g. to the case of the THG from the Higgs mode. In this case the e.m. field excites nonlinearly a single amplitude fluctuation δ∆, via a term like A 2 δ∆. 17,18,23,24 As a consequence the non-linear kernel, identified by the dashed line in Fig. 1c, is proportional to a single Higgs fluctuation, and the THG (8) is resonant when the pump frequency matches half the mode energy, i.e. ω = ω H /2 = ∆, as observed in conventional superconductors 23,26 .
The temperature dependence of the JPM non-linear kernel (7) and the corresponding THG (8) for a narrowband pulse are shown in Fig. 2a-d for different values of the pump frequency ω. Here we modelled J ⊥ (T ) and the corresponding ω J (T ) according to the out-of-plane [����] superfluid stiffness measured in Ref. 13 . As we explained, the possibility to match the resonance condition in the THG (8) depends on the relative value of the pump frequency ω with respect to ω J (T = 0) ≡ ω J,0 . In the case where ω < ω J,0 , as for ω = ω 3 in Fig. 2a, the temperature-dependence of I T HG (ω 3 ) is dominated by the maximum at the temperature where ω J (T ) = ω 3 . On the other hand, when ω ≥ ω J,0 , as it is the case for ω = ω 1 , ω 2 , the resonant excitation of the plasma mode cannot occur. However, I T HG is still non-monotonic, see Fig. 2b, due to the fact that by increasing temper- prefactor decreases, while the coth(βω J (T )/2) term increases, accounting for the thermal excitation of plasma modes. This thermal effect is particularly pronounced for the out-of-plane JPM since ω J,0 is of the same order of the critical temperature T c . The absolute value of I T HG depends also on the damping γ present in Eq. (7), which plays the same role 29 of a linear damping term in the equationsof-motion approach. In Fig. 2c,d we show the results for a temperature-dependent γ(T ) = γ 0 + r(T ), where r(T ) = r 0 e −∆/T has been taken in analogy with previous work 11 to mimics dissipative effects from normal quasiparticles. In this case the plasma resonance is progressively smeared out by increasing temperature, and for out-of-resonance conditions the THG signal looses rapidly intensity as the system is warmed up.
The THG for a field polarized along z has been measured so far only by means of a broadband pump. 13 To make a closer connection with this experimental setup we then simulated (see Ref. 29 ) the THG for a short (τ = 0.85 ps) pump pulse E p (t) with central frequency Ω/2π = 0.45 THz, as shown in Fig. 2g. The frequency spectrum of the resulting non-linear current I N L z presents then a broad peak around 3Ω, as shown in Fig. 2e. The integrated spectral weight of the 3Ω peak is shown in Fig. 2f at several temperatures. Following Ref. 13 we used Ω ω J,0 , so the narrow-band response should corresponds to the case ω = ω 2 of Fig. 2d. However, the broadband spectrum of the pump pulse enhances the response at intermediate temperatures and apart from a small deep around T = 0.2T c the signal scales with the superfluid stiffness, in good agreement with the available experimental data. In the broad-band case the nature of the non-linear kernel can also be probed via a typical pump-probe experimental setup, schematically summarized in Fig. 2g. As it has been theoretically described in Ref. 18,30 for the transmission geometry, the oscillations of the differential probe field with and without the pump δE pr (t pp ) as a function of the pump-probe time delay can be directly linked to the resonant non-linear optical kernel. In the case of the out-of-plane response (7) one then obtains (see Ref. 29 ): When the pump pulse is short enough one can approximate A 2 z (t) δ(t) and Eq. (9) shows that the differential field δE pr (t pp ) oscillates at twice the JPM frequency, and not at the frequency of the mode, as it occurs for the Higgs mode observed in conventional superconductors 31 . This prediction is confirmed when a realistic pump pulse is used in Eq. (9), as shown in Fig. 2h, which reproduces very well the 2ω J oscillations reported at low-temperature in pump-probe experiments in reflection geometry 12 .
Let us consider now the effects of a strong THz pulse polarized within the plane. In this case we can generalize the model (4) by taking into account both the twodimensional nature of the phase fluctuations in the plane and the anisotropy of penetration depth measured experimentally in cuprates, [3][4][5] where λ c 10−100λ ac depending on the material and the doping, and λ ac 2000 Å, so that ω J = c/ √ ελ ac is much larger than the out-of-plane one. Following again the microscopic derivation outlined e.g. in Ref. 27,28 we obtain where k = (k x , k y ) and we promoted the phase difference to a continuum gradient for the in-plane phase mode. To describe the non-linear coupling to the e.m. field we rely again on a quantum XY model, whose coupling constant is the effective in-plane stiffness J = 2 c 2 d/16πe 2 λ 2 . Even though the microscopically-derived phase-only action is not in general equivalent to the XY model 28 , for cuprates this can still represents a reasonable starting point 27 . By minimal-coupling substitution ∇φ(r) − (2π/Φ 0 )A we then obtain, in full analogy with Eq. (5), that: (11) By following the same steps as before we obtain a quartic action of the form (6), but the non-linear kernel becomes a tensor which admits two different K xx;xx and K xx;yy components (see Ref. 29 ): where K has the same structure of Eq. (7), provided that J ⊥ and ω J are replaced by J and ω J . The frequency and temperature dependence of K is shown in Fig. 3a. The in-plane stiffness J is taken as linearly decreasing, in analogy with experiments 3-5 . Since ω J (T = 0) ≡ ω J,0 is of the order of the eV, we only considered the case of THz pump frequencies ω i < ω J,0 . As one can see, when ω i is a fraction of ω J,0 the resonance condition ω i = ω J (T ) is still attained at temperatures where the kernel is large enough to give rise to a pronounced maximum in the THG intensity. However, when ω i ω J,0 the resonance is only attained near to T c where the prefactor has already washed out the two-plasmon resonace, and the THG scales with the superfluid stiffness. This can be easily seen from Eq. (7), since by putting ω 0 in the denominator, and considering that coth(βω J ) 1 at all relevant temperatures, from ω J ∝ J one finds The scaling of the THG intensity in the THz regime with J has several consequences. First, I T HG monotonically increases below T c , in striking contrast with the pronounced maximum one would expect for a resonance at ω i = 2∆(T ), due to the Higgs 23,24 or charge fluctuations 17,18 . Second, the superfluid stiffness appearing in the THG response is the one measured at THz frequencies. As such, due to pronounced fluctuations effects at this frequency scale, it vanishes in cuprates well above T c 16,32,33 . The I T HG at pump frequencies significantly smaller than ω J,0 closely follows the same behavior, as we exemplify in Fig. 3c where we report a simulation of the superfluid stiffness with a fluctuation tail above T c . Interestingly, both the monotonic suppression 15 and the persistence of non-linear effects above T c 15, 16 have been recently reported in THG and THz Kerr measurements in cuprate superconductors. Finally, due to the tensor structure of the in-plane kernel (12), the non-linear current associated with JPM for a pump field with a polarization angle θ with respect to the x crystallographic direction scales with: where K A1g/B1g = (K xx;xx ± K xx;yy )/2. The resulting I T HG (θ) ∝ |K(θ)| 2 is shown in Fig. 3d. According to Eq. (12), for JPM is K A1g /K B1g = 2. This result is rather different from the theoretical expectations for other collective modes. Indeed, in a single-band model, as appropriate for cuprates, the Higgs signal has only a A 1g component 17,24 . The density-fluctuations response has a largely dominant B 1g symmetry in the clean case 17 , but it becomes predominantly isotropic in the presence of even a weak disorder 22 . As a consequence, the recent observation 14 of a sizeable B 1g component at optimal doping in Bi2212 compounds cannot be simply ascribed to these collective excitations. On the other hand, it is worth noting that the ratio K A1g /K B1g = 2 for JPMs only holds within the phenomenological approach based on the quantum XY model, where the tensor admits the structure (12). Indeed, within a microscopically-derived phase-only model the interacting terms in the phase can differ from the one obtained within the XY model, as discussed for the clean case in Ref. 28 . On this view, while one expects in general an anisotropy of the non-linear JPM response, the exact value of the K A1g /K B1g ratio has to be determines within a microscopic approach. Our work establishes the theoretical framework to manipulate and detect JPMs in layered cuprates across the superconducting phase transition. The basic underlying mechanism relies on the excitation of two plasma waves with opposite momenta by an intense field. For the out-of-plane response, we support the well-established approach based on non-linear sine-Gordon equations, 6,9,10,12,13 adding a complete description of thermal effects and highlighting the possibility to tune the resonant excitation of JPMs by changing the temperature. For the in-plane response we suggest the possible relevance of JPMs to explain several puzzling aspects emerging in recent measurements in different families of cuprates. [14][15][16] An open question remains a quantitative estimate of the signal coming from the JPMs, as compared to the one due to Higgs or charge-modulated density fluctuations. Indeed, as the recent theoretical work done in the context of conventional superconductors demonstrated [19][20][21][22] , even weak disorder becomes crucial to make such a quantitative estimate, and to establish the polarization dependence of the response 22 . Here we notice that the large value of the in-plane plasma frequency comes along with a large value for the in-plane stiffness J , which controls the non-linear coupling of the JPM to the e.m. field. This suggest that especially near optimal doping, where J attains its maximum value, a two-plasmon THG signal can be comparable to other effects. On this perspective, the theoretical and experimental investigation of non-linear phenomena induced by intense THz pulses represents a privileged knob to probe relative strength of pairing and phase degrees of freedom in unconventional superconducting cuprates. The derivation of the quantum action for the phase degrees of freedom can be done following a rather standard approach, see e.g. Ref.s 1,27,28 and references therein. The basic formalism relies on the quantum-action representation of a microscopic superconducting model in the presence of long-range Coulomb interactions. The collective variables corresponding to the amplitude, phase and density degrees of freedom are introduced via an Hubbard-Stratonovich decoupling of the interacting superconducting and Coulomb term. This allows one to integrate out explicitly the fermionic degrees of freedom in order to obtain a quantum action in the collective-variables only, whose coefficients are expressed in terms of fermionic susceptibilities, computed on the SC ground state. The result for the Gaussian phase-only action in the isotropic three-dimensional case reads: (S1) Here D s = 2 c 2 /4πe 2 λ 2 andχ ρρ is the density-density susceptibility dressed at RPA level by the Coulomb interaction where χ 0 ρρ represents the bare charge susceptibility, which reduces in the static limit to the compressibility of the electron gas, i.e. χ 0 ρρ (ω m = 0, q → 0) ≡ κ. The nature of the Goldstone phase mode is dictated by the form of the charge susceptibility. For the neutral system Coulomb interactions are absent andχ ρρ in Eq. (S1) can be replaced by the bare one χ 0 ρρ . Thus, in the long-wavelength limit the pole of the Gaussian phase propagator defines, after analytical continuation to real frequencies iω m → ω + iδ, a sound-like Goldstone mode: ω 2 = (D s /κ)q 2 . On the other hand, in the presence of Coulomb interaction the long-wavelength limit of the charge compressibility (S2) scales as χ ρρ → 1/V (q). In the usual isotropic three-dimensional case V (q) = 4πe 2 /q 2 , where ε is the background dielectric constant, and one easily recovers from Eq. (S1) that where ω 2 P ≡ 4πe 2 D s / 2 ε = c 2 /λ 2 ε coincides with the usual 3D plasma frequency. In the case of cuprates one should start from a layered model where the in-plane and out-of-plane superfluid densities are anisotropic, so that the D s q 2 term in Eq. (S1) is replaced by (4D ⊥ /d 2 ) sin 2 (k z d/2) + D k 2 , with D ⊥/ = 2 c 2 /4πe 2 λ 2 c/ac . In addition, one can also introduce an anisotropic expression for the Coulomb interaction, to account for the discretization along the z direction 27 . Following e.g. the derivation of Ref. 27 one then recovers in the long-wavelength limit the two expressions (4) and (10) in the main text.
An alternative but equivalent approach is instead the one followed e.g. in Ref.s 6,9,10,12,13 , where one deals with the equations of motion for the plasmon, coupled to the electromagnetic fields. The connection between the two approaches has been derived in details in Ref. 9,10 . Once more, the authors start from a microscopic SC layered model, and integrate out the fermionic degrees of freedom in order to build up an effective action for the phase field. The effective quantum action then reads: where r i is the two-dimensional in-plane coordinate running over each SC layer and δ α is the in-plane versor along the direction α = x, y. In Eq. (S4) the quantum term accounts for the capacitive coupling C 0 = s/4πR 2 D between the planes, s and R D being, respectively, the layer thickness and the Debye length. By retaining leading orders in φ in the cosine terms the Gaussian action of Eq. (S4) describes a sound mode, in full analogy with Eq. (S1) in the absence of RPA resummation of the density response. Indeed, as emphasized in the main text, the presence of long-range interactions is crucial in order to lift the sound mode to a plasmon. In Ref. 9,10 this is achieved by adding explicitly the electric and magnetic fields, and the corresponding scalar and vector potentials. To describe the out-of-plane JPM one needs an electric fieldẼ z n,n+1 polarized perpendicularly to the planes. The magnetic field will then lie in the plane, and we can take without loss of generalityB y n,n+1 along the y in-plane direction. Hence the Gaussian action becomes: where ∆ x φ n ≡ φ n,ri (τ ) − φ n,ri+x , a and D = d + s d are, respectively, the in-plane and out-of-plane lattice spacings. By means of the Maxwell equations one can replaceẼ z n,n+1 = (Ṽ n −Ṽ n+1 )/D andB y n,n+1 = (Ã x n+1 −Ã x n )/D into Eq. (S5). The explicit integration of the e.m. potentials then leads: where α = C/C 0 and describes the full dispersion of the plasma mode as a function of q = (k z , k), with k laying in the in-plane propagation direction. The pole equation for the Gaussian phase mode, i.e. ω 2 = ω 2 P (q), is completely equivalent to the solution of the linearized sine-Gordon equation for Josephson plasma waves previously addressed in the literature 6,9,10,12,13 . In cuprates the constant α is usually very small, so the main dispersion of the plasmon comes from the last term of Eq. (S7), which accounts for the inductive coupling between planes. In this approximation, the Gaussian phase fluctuations identify a collective mode whose energy dispersion is obtained as the pole of the Guassian propagator for phase fluctuations: .

(S8)
By analytical continuation iω m → ω + iδ in Eq. (S8) we then get The relation (S9) is the same that one obtains by using the equation of motion approach discussed in Ref.s 6,9,10,12,13 . In this case, one introduces directly the variable θ n ≡ φ n − φ n+1 which represents the phase difference between nearest-neighbour layers. It is then shown to satisfy the equation of motion 6 where ∂ 2 n f n ≡ f n+1 +f n−1 −2f n is the second-order discrete differential operator along the z direction, and analogously ∂ 2 x for the x direction. As one can easily check, when sin θ n ≈ θ n Eq. (S10) admits a wave solution θ n (x = mξ 0 , t) ∝ exp[i(kx + k z nd − ωt)] where the frequency ω and the momentum q = (k, k z ) satisfy Eq. (S9). In the approach of Ref.s 6,9,10,12,13 , based on the study of the equation of motions, the electromagnetic field is completely eliminated and the non-linear effects are included by retaining the full sin θ n term in the sine-Gordon model (S10). In this case, a real q solution for a propagating waves is only possible if one retains the full momentum dispersion in Eq. (S7). In contrast, in our approach the plasma mode is first computed at Gaussian level, and then non-linear effects originate by retaining in Eq. (S5) the full cosine term appearing in Eq. (S4), which is responsible for non-linear coupling to the gauge potential. The phase mode is then integrated out in order to obtain the complete electromagnetic response, as required to describe non-linear effects in the currents, see Eq. (3) and (6) in the main text. This is indeed the same Let us consider now the case of an in-plane polarized external e.m. field. Following the same scheme adopted for the out-of-plane case we find that the in-plane fourth-order effective action is:   is the polarization-dependent tensor, whose components read: Hence the tensor components of the in-plane non-linear optical kernel are those enlisted in Eq. (12) of the main text. As a final step, let us show how the additional σ reg |ω m | term of Eq. (S12) can be added to the non-linear kernel in order to account for the effect of dissipation. In this case, the calculation is done by introducing a finite spectral function to the phase mode A(z) ≡ zσ (z 2 −ω 2 J ) 2 +z 2 σ 2 reg . One then finds that, in general, the kernel becomes: where b(z) = 1 e βz −1 is the Bose function. If σ << ω J it can be shown that Eq. (S25) can be approximated, after analytical continuation, as: Eq. (S26) is the expression used, indeed, to compute all the quantities of interest in the main text. In analogy with Ref. 6 we also assumed that where r(T ) = r 0 e −∆/T 6 , and γ 0 is a small regularization constant, which prevents the non-linear optical kernel K (diss) to be ill-defined at T = 0. Both γ 0 and r 0 parameters are fixed by looking at the number of time-resolved oscillations observed experimentally in the pump-probe set up of Ref. 12 at low temperatures. To better reproduce the experimental findings, in Fig. 2 of the main text we fixed γ 0 /2π = 0.08 THz, while r 0 = 0.3ω J,0 in panels c,d and r 0 = 0.6ω J,0 in panels e-h. There ω J,0 /2π = 0.47 THz is the out-of-plane plasma frequency at T = 0. In Fig. 3, instead, we set γ 0 = 0.1ω J,0 , where now ω J,0 /2π = 240 THz is the T = 0 value of the in-plane plasma frequency.

S3. MODELLING OF THE BROAD-BAND PUMP PULSE
For a narrow-band multicycle pulse one can assume a monochromatic incident field, and the THG is simply related to the non-linear optical kernel via Eq. (8). However, for a broad-band pulse with central frequency Ω, the THG is more generally associated with the 3Ω component in the nonlinear current 17,18 : as shown e.g. in Fig. 2e in the main text (with i, j = z and K ij = K ⊥ ) at different temperatures. Here, A z (ω) is given by the Fourier transform of A z (t) = A 0 e −t 2 /τ 2 sin (2πΩt), while A 2 z (ω) is defined as the Fourier transform of A 2 z (t). The τ = 0.85 ps and Ω/2π = 0.45 THz parameters are set in such a way that the e.m. field E z (t) ∝ −∂A z (t)/∂t well reproduces the experimental pulse profile of Ref. 13 .

S4. PUMP-PROBE CONFIGURATION
In a pump-probe experiment designed to excite the out-of-plane JPM both the pump and probe fields are polarized along z, i.e. E z = E probe (t) + E pump (t). Here we will refer for simplicity to the transmission configuration, as discussed in Ref. 18,30 , where one measures the variation δE probe (t) of the transmitted probe field with and without the pump, so that terms not explicitly depending on the pump field cancel out. This allows one to express it as δE probe (t) ∝ dt A probe z (t)K(t − t )(A pump z ) 2 (t ). By considering a fixed t g acquisition time and implementing the time-delay t pp between the pump and the probe, δE probe (t g ; t pp ) becomes a function of t pp only, as given by the first line of Eq. (9). Finally, by computing from Eq. (7) the non-linear kernel in time domain, i.e. K(t) = dω 2π K(ω)e −iωt = F (T )e −γt sin(2ω J t), we derive the last line of Eq. (9).
For the reflection geometry used in Ref. 12 the basic mechanism is the same, so that one expects that the differential reflectivity signal scales with the convolution of the non-linear kernel times the pump field squared given in Eq. (9). For the calculation of Fig. 2h in the main text we used the simulation of the broad-band pump field explained above. For the in-plane response measured in Ref. 16 , the huge frequency mismatch between the spectral components of the gauge field and 2ω J implies that only the term with t = t pp survives in the integral (9). As a consequences the