Non-linear Terahertz driving of plasma waves in layered cuprates

Francesco Gabriele, Mattia Udina, Lara Benfatto

https://doi.org/10.1038/s41467-021-21041-6, Volume: 12, Pages: null

Article Type: research-article Article History

Publisher: Nature Publishing Group UK

- Facebook
- Twitter
- Linkedin
- Whatsapp
Altmetric

Table of Contents

Introduction
Results
Methods
Supplementary information
Supplementary information

Abstract

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 cuprates. So far, THz light pulses have been used to achieve non-linear control of the out-of-plane Josephson plasma mode, whose frequency lies in the THz range. However, the high-energy in-plane plasma mode has been considered 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 markedly different thermal effects for the out-of-plane and in-plane response, linking in both cases the emergence of non-linear photonics across T_c to the superfluid stiffness. Our results show that THz light pulses represent a preferential knob to selectively drive phase excitations in unconventional superconductors.

Josephson coupling determines the superconducting phase stiffness and sets the energy scale of plasma waves. Here, the authors show that THz light can induce two-plasmon excitations of both out-of-plane and in-plane phase modes, leading however to markedly different resonant and thermal effects due to the strong anisotropy of the Josephson couplings.

Keywords

Terahertz optics, Superconducting properties and materials

Gabriele,Udina,and Benfatto: Non-linear Terahertz driving of plasma waves in layered cuprates

Introduction

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¹. 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, as 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. Although for ordinary superconductors, this energy scale is far above the THz range, in layered cuprates the existence of a weak Josephson coupling among neighboring layers^3–5 provides a natural mechanism to push down to the THz range the frequency of the interlayer Josephson plasma mode (JPM), as it was proposed long ago in order to account for the soft plasma edge appearing below T_c in standard reflectivity experiments^6–10. More recently, the possibility to manipulate such interlayer JPM by intense THz pulses has been experimentally proven^11,12, and theoretically discussed within the context of the non-linear equation of motion for the phase variable^11–16. This approach turned out to successfully capture the main features of a series of recent experiments^17,18, 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^19–21 have been discussed so far only within the context of the SC amplitude (Higgs) mode or BCS response, that consists in lattice-modulated charge fluctuations in the clean limit^22,23. Indeed, as their excitation energy scales in both cases as 2Δ, which range from 5 to 10 THz in cuprates, they appear in principle a better candidate than high-energy in-plane plasma waves. As it has been recently discussed by several authors^24–27, even small disorder affects significantly the non-linear response by triggering in general all processes mediated by the paramagnetic electronic current, that is no more conserved. This affects the relative strength of the various processes, making ultimately the Higgs mode^24–27 as well as charge/phase modes²⁷ dominant at strong disorder. The various processes can be further distinguished by their dependence on the pump polarization, and for cuprates the Higgs response is strongly isotropic at all disorder levels, whereas the BCS one has a shallow maximum for field polarized along the diagonal of square lattice unit cell²⁷. Nonetheless, the experiments show at least two features, which do not easily match our current expectation for both the Higgs and the BCS response: (i) a monotonic temperature dependence as T increases^20,21, with a persistence above T_c²¹ and (ii) a finite and doping-dependent polarization dependence with a minimum for field polarized along the diagonal¹⁹.

Here, we provide a complete theoretical description of the JPM contribution to the non-linear response of layered cuprate superconductors, focusing both on third-harmonic generation (THG) and pump-probe protocols for pump fields applied both out-of-plane and in plane. We first address the out-of-plane response and we show that the basic mechanism behind non-linear 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^22–33. 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 opposite momenta. As a consequence the excitation of out-of-plane JPMs strongly depends on the thermal probability to populate excited states and on the matching condition between the pump frequency and the JPM frequency scale, resulting in a non-monotonic dependence of the THG in temperature. We then turn our attention on the in-plane response. In this case, as the frequency scale of the in-plane JPMs is much larger than T_c and of the THz pump frequency, the THG monotonically scales in temperature with the in-plane superfluid stiffness. In addition, in contrast to the Higgs mode^22,26,27, for a light pulse polarized in the planes the signal coming from JPMs is in general anisotropic, as 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^17–21, but they also offer a perspective to design future experiments aimed at selectively tune non-linear photonic of Josephson plasma waves in layered cuprates.

Fig. 1

Non-linear excitation of phase and Higgs modes.

a Schematic view of the mexican-hat potential for the free energy F(ψ), with ψ the complex order parameter of a superconductor below T_c. A phase-gradient excitation corresponds to a shift along the minima, whereas a Higgs excitation moves the system away from the minimum. An intense light pulse with almost zero momentum can excite simultaneously two plasma waves with frequency ω_J and opposite momenta (in red) or a single Higgs fluctuation with frequency ω_H = 2Δ (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.

Results

Two-plasmon non-linear response

To elucidate the basic mechanism behind the two-plasmon non-linear response we first discuss the case of 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¹ θ_n → θ_n − (2π/Φ₀)dA_z, with θ_n = ϕ_n − ϕ_n+1, d interlayer distance and Φ₀ = 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 () naturally admits an expansion in powers of A_z to all orders:

where the explicit time convolution of Eq. () has been omitted for compactness. Here, following the same approach used so far to investigate the Higgs response^22,23,28,29, we rely on a quasi-equilibrium description, where the leading effect of the intense THz pump field is to trigger a third-order χ⁽³⁾ response mediated by plasma waves. The quantum generalization of the model () has been widely discussed within several contexts^{11,14,15,32,34,35}. Here we follow the approach of refs. ^34,35 where long-range Coulomb interactions are introduced within a layered model appropriate for cuprates (see Methods). The Gaussian quantum action for the phase mode at long wavelength has the usual form:

where

ω_{J}^{2} = c^{2} / ε λ_{c}^{2} = 8 π e d J_{c} / ℏ ε

is the energy scale of the out-of-plane JPM, iω_m = 2πmT are Matsubara frequencies and

ν

= ℏ²ε/(16πe²), with ε the background dielectric constant. In the classical limit only ω_m = 0 is relevant and one recovers the leading term of Eq. (), 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⁽⁴⁾ for the gauge field up to terms of order $O (A_{z}^{4})$ (see Methods). By coupling the gauge field A_z to the phase mode via the minimal-coupling substitution in Eq. (2) and by expanding the cosine term, one finds that:

where dots denote additional terms not relevant for the χ⁽³⁾ response. The second term in Eq. () can be treated as a perturbation with respect to

S_{⊥}^{G}

, see Supplementary Note , so that integrating out the JPM one obtains:

where

A_{z}^{2} (i ω_{m})

is defined as the Fourier transform of

A_{z}^{2} (τ)

and K_⊥(iω_m) is the non-linear optical kernel of the system, given by the convolution of two JPM propagators, as represented diagrammatically in Fig. 1b. After analytical continuation to real frequency we get:

with K₀ a constant prefactor and γ accounts for the plasmon dissipation (see Supplementary Note and ). From Eq. (), it immediately follows (see Methods) that

⟨ I_{z}^{N L} (t) ⟩ = 4 \int d t^{'} A_{z} (t) K_{⊥} (t - t^{'}) A_{z}^{2} (t^{'})

. Therefore, for a monochromatic incident field

A_{z} = A_{0} \cos (ω t)

the non-linear current admits both a term oscillating at ω, which gets mixed with the linear response, and one oscillating at 3ω, whose intensity is given by^22,23,28,29

where I₀ is an overall constant. The vanishing of the denominator in Eq. () identifies the resonance of the non-linear kernel. As the physical mechanism behind the THG is the excitation of two plasma waves, the largest I^THG in Eq. () occurs when twice the pump frequency matches the 2ω_J kernel resonance, i.e., ω = ω_J. This has to be contrasted, e.g., to the case of the THG from the Higgs mode. In this case, the electromagnetic (e.m.) field excites non-linearly a single amplitude fluctuation δΔ, via a term like A²δΔ^22,23,28,29. 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 () is resonant when the pump frequency matches half the mode energy, i.e., when ω = ω_H/2 = Δ, as observed in conventional superconductors^28,31 for strong (up to ~100 kV/cm) but not too intense fields³⁶.

Out-of-plane THG

Once derived the two-plasmon contribution to the non-linear optical kernel, let us compute the THG for a field polarized in the out-of-plane direction. The temperature dependence of the JPM non-linear kernel (7) and the corresponding THG (8) for a narrow-band pulse are shown in Fig. 2a–d for different values of the pump frequency ω. Here we modeled J_⊥(T) and the corresponding ω_J(T) according to the out-of-plane superfluid stiffness measured in ref. ¹⁸. In general, the THG for the out-of-plane response is not monotonic, as one has to face with three different temperature effects in K(2ω): (i) the suppression of J_⊥(T) and ω_J(T) with temperature; (ii) the increase of $\coth (β ω_{J})$ with temperature, owing to thermal activation of the plasmon population; (iii) the resonance condition 2ω = 2ω_J(T), that is achieved at the temperature where the (fixed) pump frequency matches the value of ω_J, and 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 ω = ω₃ in Fig. 2a, the temperature dependence of I^THG(ω₃) is dominated by the maximum at the temperature where ω_J(T) = ω₃. On the other hand, when ω ≥ ω_J,0, as it is the case for ω = ω₁, ω₂, the resonant excitation of the plasma mode cannot occur, and the temperature dependence of the THG is controlled by the opposite effects (i–ii), which lead to a non-monotonic dependence of I^THG(T). The thermal effect (ii) is particularly pronounced for the out-of-plane JPM, as ω_J,0 is of the same order of the critical temperature T_c. The absolute value of I^THG depends also on the damping γ present in Eq. (7), which has the same role of a linear damping term in the equations-of-motion approach, see Supplementary Note 1. In Fig. 2c, d, we show the results for a temperature-dependent γ(T) = γ₀ + r(T), where r(T) = r₀e^−Δ/T has been taken in analogy with previous work¹⁶ 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 rapidly loses intensity as the system is warmed up.

Fig. 2

Non-linear excitation of out-of-plane JPM.

a–d Narrow-band pulse. Temperature and frequency dependence of the non-linear kernel (7) ∣K(2ω, T)∣, normalized to its T = 0 value for the pumping frequency ω = ω_J,0, for constant a and temperature-varying c damping γ. The dashed line denotes ω_J(T)/ω_J,0. b, d show the corresponding I^THG(ω_i, T) for three values ω_i of the pump frequency, normalized to its T = 0 value for the pumping frequency ω₂ = ω_J,0. The dashed line represents J_⊥(T)/J_⊥(0). e–h Broadband pulse. e Spectrum of the non-linear current $I_{z}^{N L}$ as a function of frequency, normalized to the central frequency Ω of the pump pulse, shown explicitly in g. The intensity of the THG signal is now obtained by integrating the peak around 3Ω (gray region in e). Its temperature dependence, normalized to its T = 0 value, is shown in f, with the same color code of the curves of e. g Schematic of the pump-probe setup: a weak probe field E_pr(t) impinges on the sample with a variable time delay t_pp with respect to the intense pump pulse E_p(t). g Time-dependence of the differential probe field δE_pr(t_pp) measured with and without the pump, at different temperatures. The periodicity of the oscillations matches the 2ω_J(T) value at each temperature. Here we set ω_J,0/2π = 0.47 THz in accordance with the experiments¹⁷.

The THG for a field polarized along z has been measured so far only by means of a broadband pump¹⁸. To make a closer connection with this experimental setup we then simulated (see Methods) 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_{z}^{N L}$ 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. ¹⁸ we used Ω ≃ ω_J,0, so the narrow-band response should corresponds to the case ω = ω₂ 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.

Out-of-plane pump-probe oscillations

In the broadband 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. ^23,37 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 Methods):

where

F (T) \equiv J_{⊥}^{2} \coth (β ω_{J} / 2) / ω_{J}^{2}

. When the pump pulse is short enough one can approximate

A_{z}^{2} (t) ≃ δ (t)

and Eq. () 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³⁸. This prediction is confirmed when a realistic pump pulse is used in Eq. (), as shown in Fig. 2h, which reproduces very well the 2ω_J oscillations reported at low-temperature in pump-probe experiments in reflection geometry¹⁷.

In-plane THG

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 two-dimensional nature of the phase fluctuations in the plane and the anisotropy of penetration depth measured experimentally in cuprates^3–5, where λ_c ≃ 10−100λ_ab depending on the material and the doping, and λ_ab ≃ 2000 Å, so that $ω_{J}^{∥} = c / \sqrt{ε} λ_{a b}$ is much larger than the out-of-plane one. Following again the microscopic derivation outlined, e.g., in ref. ^34,35 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} λ_{a b}^{2}

. Even though the microscopically derived phase-only action is not in general equivalent to the XY model³⁵, for cuprates this can still represent a reasonable starting point³⁴. By minimal-coupling substitution ∇ ϕ(r) − (2π/Φ₀)A_∥ we then obtain, in full analogy with Eq. (), that:

By following the same steps as before we obtain a quartic action of the form (), but the non-linear kernel becomes a tensor, which admits two different K_xx;xx and K_xx;yy components (see Methods):

where K_∥ has the same structure of Eq. (), 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. As

ω_{J}^{∥} (T = 0) \equiv ω_{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 resonance, and the THG scales with the superfluid stiffness. This is easily seen from Eq. (), since by putting ω ≃ 0 in the denominator, and considering that

\coth (β ω_{J}^{∥}) ≃ 1

at all relevant temperatures, from

ω_{J}^{∥} \propto \sqrt{J_{∥}}

one finds

The scaling of the THG intensity in the THz regime with J_∥ has several consequences. First, I^THG monotonically increases below T_c, in striking contrast with the pronounced maximum one would expect for resonance at ω_i = Δ(T), owing to the Higgs^28,29 or BCS response^22–25,27. Second, the superfluid stiffness appearing in the THG response is the one measured at THz frequencies. As such, owing to both fluctuations effects and inhomogeneity it vanishes in cuprates well above T_c^21,39,40. Interestingly, whatever is the origin of persistence of the finite-frequency stiffness above T_c, it directly implies a persistence of the I^THG above T_c, as we exemplify in Fig. 3c, where we report a simulation of the superfluid stiffness with a finite tail above T_c (see also Supplementary Note for more details). Both the monotonic suppression²⁰ and the persistence of non-linear effects above T_c^20,21 have been recently reported in THG and THz Kerr measurements in cuprate superconductors. As explained above, they can hardly be reconciled with the typical 2Δ resonance expected for the BCS response or for the Higgs mode, both within clean^22,23 and disordered^24,25,27 models for superconductors. A second experimental finding that does not properly fit the BCS and Higgs scenario for the in-plane pump field is the polarization dependence of the response, i.e., the dependence of I^THG on the angle θ, which the pump field forms with the x crystallographic axis. Indeed, within a disordered superconducting model with a realistic band structure, the Higgs signal has an isotropic contribution, while the BCS one has a relative maximum for a field applied along the diagonal²⁷. As a consequence, the recent observation¹⁹ of a sizeable response with a minimum along the diagonal direction in optimal and over-doped Bi2212 compounds cannot be simply ascribed to these collective excitations. It is then worth exploring the polarization dependence of the JPM signal. Owing to the tensor structure of the in-plane kernel (), the non-linear current owing to JPM for a pump field with a polarization angle θ scales with:

where we introduce the standard decomposition of the kernel by means of the irreducible representation of the square lattice, i.e.,

K_{A_{1 g} / B_{1 g}} = (K_{x x; x x} \pm K_{x x; y y}) / 2

. The resulting I^THG(θ) ∝ ∣K(θ)∣² is shown in Fig. 3d. According to Eq. (), for JPM is

K_{A_{1 g}} / K_{B_{1 g}} = 2

. As mentioned above, so far JPM are the only candidate to give a B_1g contribution to the THG. In this view, the doping dependence of the anisotropy observed experimentally within the Bi2212 family¹⁹, and the reported cuprate-family dependence²⁰, both offer a potentially privileged knob to explore the relative importance of phase-fluctuation effects in cuprates. It is worth noting that the anisotropy of the kernel for JPMs follows form the intrinsics anisotropy of the two-plasmon excitation process, and the ratio

K_{A_{1 g}} / K_{B_{1 g}} = 2

only holds within the phenomenological approach based on the quantum XY model, where the overall coupling constant of the action () is isotropic. However, 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. ³⁵. As a consequence, although one expects, in general, an anisotropy of the non-linear JPM response, the value of the

K_{A_{1 g}} / K_{B_{1 g}}

ratio could also be influenced by microscopic details.

Fig. 3

Non-linear excitation of in-plane JPM.

a Temperature and frequency dependence of the non-linear kernel ∣K_∥(2ω, T)∣, normalized to its T = 0 value for the pumping frequency ω = ω_J,0, for constant damping γ. The dashed line denotes $ω_{J}^{∥} (T) / ω_{J, 0}^{∥}$ . b I^THG(ω_i, T) for three ω_i values of the pump frequency, marked in a, normalized to I^THG(ω_i, 0). The dashed line represents J_∥(T)/J_∥(0). As one can see, when $ω_{J, 0}^{∥} / ω_{i}$ is increased the THG intensity progressively approaches the temperature dependence of the stiffness. c Effect of superconducting fluctuations on the THG. Here the dashed line simulates the experimental behavior of the J_∥(T) measured at THz frequencies, with a pronounced tail above T_c. When $ω_{J, 0}^{∥} / ω_{i}$ is increased also the THG signal survives above T_c, following the fluctuating stiffness. d Angular dependence of I^THG(θ) = ∣K(θ)∣², where K(θ) is given by Eq. (14).

Discussion

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^{11,14,15,17,18}, 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^19–21. Although for the out-of-plane response the strong incoherent quasiparticle transport automatically suppresses all electronic mechanisms, leaving the JPM as the only plausible candidate to explain non-linear effects, for the in-plane case an open question remains a quantitative estimate of the signal coming from the JPMs, as compared with the one owing to the Higgs or to BCS quasiparticle excitations. Indeed, as the recent theoretical work on disordered superconducting models demonstrated^24–27, even weak disorder becomes crucial to estimate the relative strength of the various possible processes, and to establish the polarization dependence of the response²⁷. Nonetheless, as we have shown, even in the absence of a quantitative estimate of the hierarchy of the various effects the temperature and polarization dependence of the non-linear response can be used to discriminate different contributions. In our modeling, 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 suggests that especially near optimal doping, where J_∥ attains its maximum value, a two-plasmon THG signal can be comparable to other effects. An interesting additional question is a possibility that a finite supercurrent triggered by a very strong THz field, as the one recently discussed within the context of second-harmonic generation in conventional superconductors^31,41, could also allow for single-plasmon excitations processes. From this perspective, the theoretical and experimental investigation of non-linear phenomena induced by intense THz pulses represents a privileged knob to probe the relative strength of pairing and phase degrees of freedom in unconventional superconducting cuprates.

Methods

Effective quantum action

The derivation of the quantum action for the phase degrees of freedom can be done following a rather standard approach, see, e.g., refs. ^1,34,35 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:

Here D_s = ℏ²c²/4πe²λ² and

{\tilde{χ}}_{ρ ρ}

is the density–density susceptibility dressed at RPA level by the Coulomb interaction V(q):

where

χ_{ρ ρ}^{0}

represents the bare charge susceptibility, which reduces in the static limit to the compressibility of the electron gas, i.e.

χ_{ρ ρ}^{0} (ω_{n} = 0, q \to 0) \equiv κ

. 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

{\tilde{χ}}_{ρ ρ}

in Eq. () 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ω_n → ω + iδ, a sound-like Goldstone mode: ω² = (D_s/κ)q². On the other hand, in the presence of Coulomb interaction the long-wavelength limit of the charge compressibility () scales as

{\tilde{χ}}_{ρ ρ} \to 1 / V (q)

. In the usual isotropic three-dimensional case V(q) = 4πe²/q², where ε is the background dielectric constant, and one easily recovers from Eq. () that

where

ω_{P}^{2} \equiv 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_sq² term in Eq. () is replaced by

(4 D_{⊥} / d^{2}) \sin^{2} (k_{z} d / 2) + D_{∥} k_{∥}^{2}

, with

D_{⊥ / ∥} = ℏ^{2} c^{2} / 4 π e^{2} λ_{c / a c}^{2}

. In addition, one can also introduce an anisotropic expression for the Coulomb interaction, to account for the discretization along the z direction³⁴. Following, e.g., the derivation of ref. ³⁴ one then recovers in the long-wavelength limit the two expressions () and (). Notice that at long-wavelengths the result () coincides also with the one based on the non-linear sine-Gordon equations, as shown in refs. ^11,14,15. In this case, however the effect of long-range forces is included via the coupling to the electromagnetic gauge and scalar potentials, which are eliminated to derive the equations of motion for the phase variables. Further technical details on this analogy are provided in Supplementary Note .

Computation of the non-linear kernel

The current I_α in the α = (x, y, z) direction is defined as usual as the functional derivative with respect to A_α of the action S_A. Thus, to compute the third-order contribution to I_z in Eq. (3) we need to expand the e.m. action up to terms of order $O (A_{z}^{4})$ . The coupling term of the JPM to $A_{z}^{2}$ in Eq. (5) leads to a $A_{z}^{4}$ contribution after integrating out the plasmon, see Supplementary Note 2. This is represented by the Feynmann diagram of Fig. 1b. Here each solid line denotes the Gaussian phase mode, obtained by Eq. (4) as $⟨ ∣ ϕ (i ω_{m}, k_{z}) ∣^{2} ⟩ = {\{4 \sin^{2} (k_{z} d / 2) [ω_{m}^{2} + ω_{J}^{2}]\}}^{- 1}$ . With straightforward calculations one gets:

where

A_{z}^{2} (i ω_{m}) = \sum_{i ω_{m}^{'}} A_{z} (i ω_{m}^{'}) A_{z} (i ω_{m} - i ω_{m}^{'})

is the Fourier transform on

A_{z}^{2} (τ)

. Eq. () coincides with Eq. (), once defined

K_{⊥} (i ω_{m}) = K_{0} \frac{J_{⊥}^{2}}{ω_{J}} \frac{\coth (β ω_{J} / 2)}{4 ω_{J}^{2} + ω_{m}^{2}}

. After analytical continuation iω_m → ω + iδ to real frequencies one then recovers Eq. (). From Eq. () one directly derives

I_{z} (i ω_{n}) = - \partial S_{A}^{(4)} / \partial A_{z} (- i ω_{n}) = - 4 \sum_{i ω_{m}^{'}} A_{z} (i ω_{n} - i ω_{m}^{'}) K (i ω_{m}^{'}) A_{z}^{2} (i ω_{m}^{'})

, we used the parity of the kernel to write the four possible derivatives in the same way. After analytical continuation to real frequency one has:

whose Fourier transform to real time gives

⟨ I_{z}^{N L} (t) ⟩ = - 4 \int d t^{'} A_{z} (t) K_{⊥} (t - t^{'}) A_{z}^{2} (t^{'})

, as stated in the main text. For a monochromatic field, A(ω) is proportional to a delta function peaked at the pump frequency, and one recovers Eq. (). Notice that the infinitesimal positive δ in the analytical continuation of the kernel is promoted here to a finite and temperature-dependent value γ(T) = γ₀ + r₀e^−Δ/T to account for plasmon dissipative effects, as explained in Supplementary Note . To better reproduce the pump-probe experimental findings¹⁷, in Fig. 2 we fixed γ₀/2π = 0.08 THz, while r₀ = 0.3ω_J,0 in panels c,d and r₀ = 0.6ω_J,0 in panels e–h. Here ω_J,0/2π = 0.47 THz is the out-of-plane plasma frequency at T = 0. In Fig. 3, instead, we set γ₀ = 0.1ω_J,0, where now ω_J,0/2π = 240 THz is the T = 0 value of the in-plane plasma frequency.

For what concerns the in-plane JPM, we follow the same procedure starting from the interaction term of Eq. (11). In this case, the quartic action has a structure similar to Eq. (18) provided that K_⊥ is replaced by a two-component tensor:

where K_ii;jj(iω_m) = M_ijK_∥, and

M_{i j} = \sum_{k} \frac{k_{i}^{2} k_{j}^{2}}{k^{4}}

. As a consequence, up to an overall normalization factor, one has that

M_{x x} = M_{y y} \propto \int_{0}^{2 π} d ϕ \cos^{4} ϕ = 3 π / 4

while

M_{x y} = M_{y x} \propto \int_{0}^{2 π} d ϕ \cos^{2} ϕ \sin^{2} ϕ = π / 4

, leading to Eq. (). If θ is the angle of the pump field in the plane, i.e.,

A = (A \cos θ, A \sin θ, 0)

, then the polarization angle-dependence of the in-plane non-linear optical kernel is²²:

that can be rewritten in the form of Eq. ().

Broadband 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 broadband pulse with central frequency Ω, the THG is more generally associated with the 3Ω component in the non-linear current^22,23. We then computed the non-linear current from Eq. (19) by using a realistic pump spectrum A(ω), obtained by Fourier transform of $A_{z} (t) = A_{0} e^{- t^{2} / τ^{2}} \sin (Ω t)$ . The result for I_z is shown in Fig. 2e at different temperatures. 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 reproduces accurately the experimental pulse profile of ref. ¹⁸.

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_pr(t) + E_p(t). Here, we will refer for simplicity to the transmission configuration, as discussed in ref. ^23,37, where one measures the variation δE_pr(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_{pr} (t) \propto \int d t^{'} A_{z}^{pr} (t) K (t - t^{'}) {(A_{z}^{p})}^{2} (t^{'})$ . By considering a fixed t_g acquisition time and implementing the time delay t_pp between the pump and the probe, δE_pr(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) = \int \frac{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. ¹⁷ 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 simulations in Fig. 2e–h, we used the broadband pump pulse described above. For the in-plane response measured in ref. ²¹, 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 consequence, the oscillations are absent and δE_pr(t_pp) simply scales as the square of the pump field, modulated by F(T) and by the polarization encoded in the kernel (12). Indeed, if the pump field forms an angle θ with the x axis and the probe is applied, e.g., along the x axis, from Eq. (9), properly generalized for the planar configuration, one easily sees that $δ E_{x} ~ K_{x x; x x} \cos^{2} θ + K_{x x; y y} \sin^{2} θ = K_{A_{1 g}} + K_{B_{1 g}} \cos (2 θ)$ . This is exactly the decomposition used to analyzed the transient reflectivity measured in ref. ¹⁹.

Peer review information: Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

The online version contains supplementary material available at 10.1038/s41467-021-21041-6.

Acknowledgements

We acknowledge useful discussions with C. Castellani, A. Cavalleri, and D. Nicoletti. We thank the authors of ref. ¹⁸ for providing us with the experimental data used to estimate J_⊥(T) in Fig. 2. This work has been supported by the Italian MAECI under the Italian-India collaborative project SUPERTOP-PGR04879, by the Italian MIUR project PRIN 2017 No. 2017Z8TS5B, by Regione Lazio (L.R. 13/08) under project SIMAP and by Sapienza University under project Ateneo 2019 (Grant No. RM11916B56802AFE).

Author contributions

F.G. and M.U. contributed equally to this work. L.B. conceived the project and supervised its development. F.G., M.U. and L.B. performed the analytical calculations. F.G. and M.U. performed the numerical simulations. L.B. wrote the manuscript with inputs from all the authors.

Data availability

All data generated during this study are included in this published article (and its supplementary information files).

Competing interests

The authors declare no competing interests.

References

Nagaosa, N. & Heusler, S. Quantum Field Theory in Condensed Matter Physics, Texts and monographs in physics (Springer, 1999).

Anderson

. Random-phase approximation in the theory of superconductivity.

doi: 10.1103/PhysRev.112.1900

Shibauchi

. Anisotropic penetration depth in La_2−xSr_xCuO₄.

doi: 10.1103/PhysRevLett.72.2263

Panagopoulos

. Anisotropic magnetic penetration depth of grain-aligned HgBa₂Ca₂Cu₃O_8+δ.

doi: 10.1103/PhysRevB.53.R2999

Hosseini

. Survival of the d -wave superconducting state near the edge of antiferromagnetism in the cuprate phase diagram.

doi: 10.1103/PhysRevLett.93.107003

Tamasaku

, Nakamura

, Uchida

. Charge dynamics across the cuo₂ planes in la_2−xsr_xcuo₄.

doi: 10.1103/PhysRevLett.69.1455

Homes

, Timusk

, Liang

, Bonn

, Hardy

. Optical conductivity of c axis oriented yba₂cu₃o_6.70: evidence for a pseudogap.

doi: 10.1103/PhysRevLett.71.1645

Kim

. Strong damping of the c-axis plasmon in high-tc cuprate superconductors.

doi: 10.1016/0921-4534(95)00198-0

Basov

, Timusk

, Dabrowski

, Jorgensen

. c-axis response of yba₂cu₄o₈: a pseudogap and possibility of josephson coupling of cuo₂ planes.

doi: 10.1103/PhysRevB.50.3511

10.

van der Marel

, Tsvetkov

. Transverse optical plasmons in layered superconductors.

doi: 10.1007/BF02548125

11.

Savel’ev

, Yampol’skii

, Rakhmanov

, Nori

. Terahertz josephson plasma waves in layered superconductors: spectrum, generation, nonlinear and quantum phenomena.

doi: 10.1088/0034-4885/73/2/026501

12.

Laplace

, Cavalleri

. Josephson plasmonics in layered superconductors.

13.

Koyama

, Tachiki

. I − V characteristics of josephson-coupled layered superconductors with longitudinal plasma excitations.

doi: 10.1103/PhysRevB.54.16183

14.

Machida

, Koyama

, Tachiki

. Dynamical breaking of charge neutrality in intrinsic josephson junctions: common origin for microwave resonant absorptions and multiple-branch structures in the I − V characteristics.

doi: 10.1103/PhysRevLett.83.4618

15.

Machida

, Koyama

, Tanaka

, Tachiki

. Theory of the superconducting phase and charge dynamics in intrinsic josephson-junction systems: microscopic foundation for longitudinal josephson plasma and phenomenological dynamical equations.

doi: 10.1016/S0921-4534(99)00612-7

16.

Savel’ev

, Rakhmanov

, Yampol’skii

, Nori

. Analogues of nonlinear optics using terahertz josephson plasma waves in layered superconductors.

doi: 10.1038/nphys358

17.

Rajasekaran

. Parametric amplification of a superconducting plasma wave.

doi: 10.1038/nphys3819

18.

Rajasekaran

. Probing optically silent superfluid stripes in cuprates.

doi: 10.1126/science.aan3438

19.

Katsumi

. Higgs mode in the d -wave superconductor Bi₂Sr₂CaCu₂O_8+x driven by an intense terahertz pulse.

doi: 10.1103/PhysRevLett.120.117001

20.

Chu

. Phase-resolved Higgs response in superconducting cuprates.

doi: 10.1038/s41467-020-15613-1

21.

Katsumi, K., Li, Z. Z., Raffy, H., Gallais, Y. & Shimano, R. Superconducting fluctuations probed by the Higgs mode in Bi₂Sr₂CaCu₂O_8+x thin films, Phys. Rev. B102, 054510 (2020)

22.

Cea

, Castellani

, Benfatto

. Nonlinear optical effects and third-harmonic generation in superconductors: cooper pairs versus higgs mode contribution.

doi: 10.1103/PhysRevB.93.180507

23.

Udina

, Cea

, Benfatto

. Theory of coherent-oscillations generation in terahertz pump-probe spectroscopy: from phonons to electronic collective modes.

doi: 10.1103/PhysRevB.100.165131

24.

Silaev

. Nonlinear electromagnetic response and higgs-mode excitation in bcs superconductors with impurities.

doi: 10.1103/PhysRevB.99.224511

25.

Murotani

, Shimano

. Nonlinear optical response of collective modes in multiband superconductors assisted by nonmagnetic impurities.

doi: 10.1103/PhysRevB.99.224510

26.

Shimano

, Tsuji

. Higgs mode in superconductors.

doi: 10.1146/annurev-conmatphys-031119-050813

27.

Seibold, G., Udina, M., Castellani, C. & Benfatto, L. Third harmonic generation from collective modes in disordered superconductors. Phys. Rev. B. 103, 014512 (2021).

28.

Matsunaga

. Light-induced collective pseudospin precession resonating with higgs mode in a superconductor.

doi: 10.1126/science.1254697

29.

Tsuji

, Aoki

. Theory of anderson pseudospin resonance with higgs mode in superconductors.

doi: 10.1103/PhysRevB.92.064508

30.

Schwarz

. Classification and characterization of nonequilibrium higgs modes in unconventional superconductors.

doi: 10.1038/s41467-019-13763-5

31.

Yang

. Lightwave-driven gapless superconductivity and forbidden quantum beats by terahertz symmetry breaking.

doi: 10.1038/s41566-019-0470-y

32.

Sun

, Fogler

, Basov

, Millis

. Collective modes and terahertz near-field response of superconductors.

doi: 10.1103/PhysRevResearch.2.023413

33.

Vaswani, C. et al. Light quantum control of persisting Higgs modes in iron-based superconductors. Preprint at https://arxiv.org/abs/2011.13036 (2020a).

34.

Benfatto

, Caprara

, Castellani

, Paramekanti

, Randeria

. Phase fluctuations, dissipation, and superfluid stiffness in d-wave superconductors.

doi: 10.1103/PhysRevB.63.174513

35.

Benfatto

, Toschi

, Caprara