Recently developed terahertz (THz) two-dimensional coherent spectroscopy (2DCS) is a powerful technique to obtain materials information in a fashion qualitatively different from other spectroscopies. Here, we utilized THz 2DCS to investigate the THz nonlinear response of conventional superconductor NbN. Using broadband THz pulses as light sources, we observed a third-order nonlinear signal whose spectral components are peaked at twice the superconducting gap energy $2\mathrm{\Delta }$. With narrow-band THz pulses, a THz nonlinear signal was identified at the driving frequency $\mathrm{\Omega }$ and exhibited a resonant enhancement at temperature when $\mathrm{\Omega }=2\mathrm{\Delta }$. General theoretical considerations show that such a resonance can arise only from a disorder-activated paramagnetic coupling between the light and the electronic current. This proves that the nonlinear THz response can access processes distinct from the diamagnetic Raman-like density fluctuations, which are believed to dominate the nonlinear response at optical frequencies in metals. Our numerical simulations reveal that, even for a small amount of disorder, the $\mathrm{\Omega }=2\mathrm{\Delta }$ resonance is dominated by the superconducting amplitude mode over the entire investigated disorder range. This is in contrast to other resonances, whose amplitude-mode contribution depends on disorder. Our findings demonstrate the unique ability of THz 2DCS to explore collective excitations inaccessible in other spectroscopies.

We compute the nonlinear current J^NL induced by an intense THz field in superconducting cuprates. We find that J^NL is dominated by the BCS response, showing comparable para- and diamagnetic contributions, in agreement with θ-resolved THG experiments.

Multiband superconductors can host collective excitations with marked differences with respect to their single-band counterpart. We first study the spectrum of collective amplitude fluctuations in a clean two-band superconductor, showing that the spectral weight of the Higgs mode rapidly deviates from the naive extension of the single band case as the interband coupling is turned on. These results are then used to critically analyze the nonlinear optical response in ${\mathrm{MgB}}_2$, providing an explanation for the apparently contradictory results of recent experiments, pointing towards a selective relevance either of the Leggett mode or of the amplitude fluctuations at twice the lower gap. By using exact numerical simulations and realistic estimates of disorder we compute the relative contribution of the quasiparticle, amplitude, and phase fluctuations to the nonlinear optical response. We show that at low pumping frequency only the resonance at twice the smaller gap emerges, as due to the BCS response, while the Leggett mode dominates only in a narrow range of higher pumping frequencies matching its low-temperature value. Our findings provide a fresh perspective on the potential of nonlinear THz spectroscopy to detect collective modes in other multiband systems, such as iron-based superconductors.

Josephson plasmons in layered superconductors represent a natural source of optical nonlinearity, thanks to their intrinsically anharmonic nature. Here we derive the selection rules behind nonlinear plasmonics showing its dependence on plasmonic branches hidden to other spectroscopies, such as RIXS. We benchmark our results for the case of layered cuprates, showing how in a layered system the combined effect of plasmon dispersion and light polarization can move the resonance of the bilayer system away from the plasma edge measured in linear spectroscopy. Our results demonstrate the dependence of the nonlinear THz response on the convoluted plasmon dispersion in a momentum region complementary to RIXS, and offer a possible perspective for the generation of THz pulses by artificially designed Josephson heterostructures.

Electromagnetic waves propagating in a layered superconductor with arbitrary momentum, with respect to the main crystallographic directions, exhibit an unavoidable mixing between longitudinal and transverse degrees of freedom. Here we show that this basic physical mechanism explains the emergence of a well-defined absorption peak in the in-plane optical conductivity when light propagates at small tilting angles relative to the stacking direction in layered cuprates. More specifically, we show that this peak, often interpreted as a spurious leakage of the c-axis Josephson plasmon, is instead a signature of the true longitudinal plasma mode occurring at larger momenta. By combining a classical approach based on Maxwell’s equations with a full quantum derivation of the plasma modes based on modeling the superconducting phase degrees of freedom, we provide an analytical expression for the absorption peak as a function of the tilting angle and light polarization. We suggest that an all-optical measurement in tilted geometry can be used as an alternative way to access plasma-wave dispersion, usually measured by means of large-momenta scattering techniques like resonant inelastic X-ray scattering (RIXS) or electron energy loss spectroscopy (EELS).

The term two-dimensional coherent spectroscopy (2DCS) usually refers to experimental setups where a coherently generated electric field in a sample is recorded over many runs as a function of two time variables: the delay $τ$ between two consequent excitation pulses and the time $t$ over which the signal is emitted. While its implementation in the femtosecond time domain for studying vibrational molecular states has been developed for over two decades, its experimental application in the THz domain to interacting electronic systems remains in its infancy. This work provides a general theoretical framework for describing and interpreting 2DCS using a many-body language based on a perturbative diagrammatic expansion, as widely applied in linear spectroscopy. Focusing on centrosymmetric systems, we show that interpreting the 2D maps can be recast into two complementary problems. The first is the evaluation of a third-order response function to the gauge field. In the velocity gauge, this leads to semi-analytical expressions that both reduce computational complexity and assist in assigning spectral features to microscopic processes, as shown using a toy model of electrons undergoing a charge-density wave transition. The second is a careful treatment of multi-wave propagation effects, which, in bulk systems, can obscure the intrinsic nonlinear response, demonstrated here for soft superconducting Josephson plasmons. Our results provide a solid foundation for extending 2DCS to complex interacting systems and offer a flexible method to realistically model nonlinear responses across arbitrary spectral widths.