Mode-resolved, non-local electron–phonon coupling in two-dimensional spectroscopy (2025)

Interactions among constituent particles in many-body systems govern the physics of condensed-matter systems. Electron–phonon couplings (EPCs), in particular, affect the temperature dependence of carrier mobility and resistance, electronic energy band structures, phonon dispersion relations and lifetimes of those particles in solids¹. The roles EPCs have in ordered phases, such as conventional superconductivity, charge-density wave order and metal–insulator transitions are the subject of vast research fields². Understanding the role of phonons in high-temperature superconductors³ requires a careful modelling of the microscopic nature of EPCs, including non-local effects that lead to a variation in EPCs with electron momentum k (ref. ⁴). In fact, a strong coupling of phonons to the kinetic energy of electrons rather than potential energy has been predicted to lead to considerably higher critical temperatures for phonon-mediated superconductors^5,6. It is, therefore, crucial to identify materials that go beyond the paradigm of k-independent EPC in Holstein or Fröhlich models⁷. Moreover, via specific phonon modes, one can modulate the band structures^8,9,10,11, transport properties^12,13,14 and topology^15,16 on (sub-)picosecond timescales and with high precision in atomic motions. Such an ultrafast phonon control allows for achieving electronic states not accessible in static materials^17,18.

To advance our fundamental understanding of emergent phenomena and capability to exploit new functionalities, it is essential to differentiate contributions from individual phonon modes and different electronic energies to the overall EPC strength. Most experimental approaches, however, provide rather indirect information on EPCs without an explicit phonon-mode resolution. Conventional techniques, such as transport measurements, heat capacity or angle-resolved photoelectron spectroscopy, measure the renormalization of response functions or the band structure due to phonons. Such quantities typically include the average effect of many phonon modes and are also subject to other phenomena such as electron–electron interaction, rendering the phonon-mode-resolved EPC strength extraction impossible or unreliable^19,20. Although inelastic X-ray scattering^21,22 and time- and angle-resolved photoelectron spectroscopy have been proposed^19,23 for this purpose, these approaches have practical limitations coming from their intrinsic surface sensitivity, resolution at low energies and additional material-specific requirements^3,21,24,25.

Here we develop a two-dimensional (2D) spectroscopy protocol that allows for the direct extraction of the EPC matrix elements for individual optical phonon modes and different electron energies (or electron momenta for a known dispersion) simultaneously. In this 2D electron–phonon coupling spectroscopy (2D EPC), the EPC strengths of a solid manifest themselves as signal peaks in a measured 2D spectrum, where each of the two axes corresponds to the zone-centre phonon frequencies and the electron energies. We demonstrate how the EPC matrix elements are extracted from the measured 2D EPC signal without requiring extensive modelling. The interpretation of the measurement is facilitated by the fact that the signal vanishes in the absence of EPC. This is a direct manifestation of a general principle in 2D spectroscopy: cross-peaks can only arise because of correlations between the excitation and detection frequencies²⁶. Since our approach is all-optical and bulk sensitive, it is applicable to almost every kind of gapped solid. It is distinctly different from techniques that rely on the phonon-induced relaxation of electronic states like the recently demonstrated time- and angle-resolved photoemission experiments²⁴, which are limited to electronic states with a restricted scattering phase space and a sufficiently strong coupling to phonons. By contrast, our signal exclusively arises from electron–phonon correlations, making our experiment sensitive even to weak couplings and avoiding any averaging over electronic states.

To showcase our technique, we measure the EPC strength of two pronounced phonon modes in methylammonium lead iodide perovskite (MAPI) as a function of the electron energy. We find a surprisingly strong dependence on the electron energy, including a vanishing of the EPC strength at the electronic Γ point. This result indicates a pronounced coupling of phonons to the electronic kinetic energy in the sample akin to a Su–Schrieffer–Heeger-type (SSH-type) EPC. Hence, our technique is able to distinguish different microscopic EPC models by detecting the characteristic electron momentum dependence of the EPC strength, which is beyond the reach of conventional methods. In addition, we perform temperature- and polarization-dependent measurements to highlight the potential of our 2D EPC technique for studying the evolution of mode-resolved EPC as a function of various phase-space parameters (for example, temperature, pressure, strain, doping and chemical composition).

The basic concept of 2D EPC spectroscopy is illustrated in the schematic shown in Fig. 1a, where a series of three laser pulses interacts with the material. These interactions induce a polarization whose contribution at the third order in the electric field of the light waves is subsequently measured. To generate correlated excitations of specific electronic states and phonon modes, we use one pulse with broadband THz frequency (denoted by THz in Fig. 1a) and two pulses with optical frequencies. The individual frequencies of these two optical pulses are chosen to be non-resonant with any elementary excitation, whereas their sum is resonant with an electronic interband transition, as illustrated by the energy diagram shown in Fig. 1b. One of the two optical pulses has a broad bandwidth (denoted by BB), and the other, a narrow bandwidth (denoted by NB); therefore, the NB spectrum can be approximated by a delta function (for its detailed use in analysis, see Methods).

Figure 1c displays the time evolution of the initial state (denoted by the density matrix |a〉〈a|) subject to successive radiation–matter interactions (coloured arrows). The final state can then be obtained by summing over all possible Liouville-space pathways²⁷ (Fig. 1c). The radiation–matter interaction is described by the semiclassical Hamiltonian under the dipole approximation as H_int(t′) = –E(r, t′) ⋅ μ, where E(r, t′) is the position- and time-dependent electric field and μ is the dipole operator. The induced polarization at the third order in the electric field, P⁽³⁾(r, t′), is determined by the expectation value of the third-order density matrix, ρ⁽³⁾(t′), and can be expressed in terms of the third-order response function S_EPC(t₃, t₂, t₁) (Methods and ref. ²⁷). After a 2D Fourier transformation over times t and τ defined in Fig. 1a (the Fourier frequencies are denoted by ω_t/τ), the 2D EPC signal field E_EPC(ω_t, ω_τ) is obtained in terms of the response function S_EPC(ω_t, ω_τ) and the electric fields of the incident THz (BB) pulses, E_THz (E_BB) (Methods):

We will see below that the response function S_EPC(ω_t, ω_τ) is a measure of the coupling strength of phonons and electronic excitations at frequencies ω_τ and ω_t.

The starting point of our theoretical description of the material is the Hamiltonian H = H₀ + H_ep, which contains a non-interacting part \(\scriptstyle{H}_{0}={\sum }_{k,n}{\varepsilon }_{n,{\bf{k}}}{c}_{n,{\bf{k}}}^{\dagger }{c}_{n,{\bf{k}}}+{\sum }_{{\bf{q}},\lambda }{\omega }_{{\bf{q}},\lambda }({b}_{{\bf{q}},\lambda }^{\dagger }{b}_{{\bf{q}},\lambda }+1/2)\), as well as a generic electron–phonon interaction to the first order in nuclear displacements^1,7:

Here M_k,q,n,λ is the dimensionless coupling strength, the operator c_n,k (\({c}_{n,{\bf{k}}}^{\dagger }\)) annihilates (creates) an electron with band index n and energy ε_n,k, the operator \({b}_{{\bf{q}},\lambda }\,({b}_{{\bf{q}},\lambda }^{\dagger })\) annihilates (creates) a phonon of mode λ with occupation number v and energy ω_q,λ, and we have set ℏ = 1. We specifically consider an electronic two-band model with n = g, e, labelling the ground- and excited-state bands, respectively. The energy-level diagram of this model at fixed electron momentum k and q = 0 corresponding to the single-particle eigenstates \(| {n}_{{\bf{k}}},{\chi }_{n,\lambda }^{v}\left.\right\rangle\) is shown in Fig. 1b using the short-hand notation \(| n,v\left.\right\rangle \equiv | {n}_{{\bf{k}}},{\chi }_{n,\lambda }^{v}\left.\right\rangle\).

Then, we proceed by absorbing the EPC term of the Hamiltonian into a dressed phonon, redefining the phonon operator b_q,λ → B_q,λ as

The Hamiltonian becomes diagonal in the transformed phonon basis, \(H={\sum }_{{\bf{q}},\lambda }{\omega }_{{\bf{q}},\lambda }({B}_{{\bf{q}},\lambda }^{\dagger }{B}_{{\bf{q}},\lambda }+\frac{1}{2})+{H}_{{\rm{el}}}\), where the electronic part of the Hamiltonian H_el now also contains phonon-mediated electron–electron interactions. Also, we have used the identity \({M}_{{\bf{k}}+{\bf{q}},-{\bf{q}},n,\lambda }^{* }={M}_{{\bf{k}},{\bf{q}},n,\lambda }\), which follows from the hermiticity of equation (2). We denote electrons in the many-body ground state (that is, forming a Fermi sea) by |g〉; and those in the optically excited states with a single electronic interband transition at momentum k by \(| {e}_{{\bf{k}}}\left.\right\rangle ={c}_{e,{\bf{k}}}^{\dagger }{c}_{g,{\bf{k}}}| g\left.\right\rangle\). In the following, we drop the index q, focusing on the q = 0 sector relevant to this optical approach. In this sector, the dressed phonon operator B_q=0,λ is simply the original phonon operator shifted by a constant that depends on the occupation of the electronic bands. Henceforth, we refer to this shifted operator as the ‘phonon’ with the implication that its eigenstates depend on the electronic state. For instance, the phonon vacuum at q = 0 for the electronic ground state is defined by

By contrast, the phonon vacuum for an electronic excited state |e_k〉 is defined by

which means that \(\scriptstyle| {\chi }_{e,{\bf{k}},\lambda }^{\,v = 0}\left.\right\rangle\) takes the form of a coherent state when expressed in the basis \(\scriptstyle| {\chi }_{g,\lambda }^{\,v = 0}\left.\right\rangle\).

The phonon-mode- and electron-momentum-dependent EPC strengths M_k,n,λ can be associated with the measured response function S_EPC(t₃, t₂, t₁). To this end, we consider all the possible Liouville-space pathways representing different time orderings of overlapping pulses and different intermediate states of the multiphoton process (Supplementary Texts 1.3 and 1.4). The magnitude of the response function is proportional to all the transition matrix elements, \({\mu }_{{n}^{{\prime} }{v}^{{\prime} },nv}\)\(=\left\langle \right.{n}_{{\bf{k}}}^{{\prime} },{\chi }_{{n}^{{\prime} },\lambda }^{\,{v}^{{\prime} }}| \mu | {n}_{{\bf{k}}},{\chi }_{n,\lambda }^{\,v}\left.\right\rangle\), involved in a given pathway. For the ‘non-rephasing’ contribution, where the phase of the induced polarization is acquired continuously during t₁ and t₃, the three diagrams \({R}_{2}^{* }\), \({R}_{3}^{* }\) and R₄ (Fig. 1c) need to be taken into account (Supplementary Fig. 1 and ref. ²⁷). The nonlinear response function is given by the sum over all possible pathways (meaning all possible diagrams and intermediate states), and we find

where the three sums only run over states |a〉 = |n, v〉 allowed by the respective diagrams in Fig. 1c. Here we use the short-hand notation m_abcd = μ_baμ_adμ_cbμ_dc and ω_ba denotes the energy difference between states |b〉 and |a〉.

In the absence of EPC, the transition matrix elements are simply \({\mu }_{{n}^{{\prime} }{v}^{{\prime} },nv}={\delta }_{{v}^{{\prime} }v}\left\langle \right.{n}_{{\bf{k}}}^{{\prime} }| {\mu }_{e}| {n}_{{\bf{k}}}\left.\right\rangle +{\delta }_{{n}^{{\prime} }n}\left\langle \right.{\chi }_{n,\lambda }^{{v}^{{\prime} }}| {\mu }_{p}| {\chi }_{n,\lambda }^{v}\left.\right\rangle\), with \(\left\langle \right.{\chi }_{n,\lambda }^{{v}^{{\prime} }}| {\mu }_{p}| {\chi }_{n,\lambda }^{v}\left.\right\rangle\)\(\propto {\delta }_{{v}^{{\prime} },v\pm 1}\). Hence, a single photon either leaves the phonon number unchanged (interband transition) or changes it by ±1 (intraband transition). The first photon in Fig. 1c resonantly excites a single phonon, and so, we can fix |a〉 = |g, 0〉 and |b〉 = |g, 1〉 for all three diagrams. The second and third photons do not need to be resonant, but if we focus on the non-rephasing contribution (that is, the same sign between ω_τ and ω_t), we only find one possible set of intermediate states, namely, |c〉 = |e, 1〉 and |d〉 = |g, 1〉, for the first two diagrams. For the third diagram, there are two possibilities, namely, |c〉 = |g, 0〉 or |c〉 = |e, 1〉 and |d〉 = |e, 0〉. The detailed diagrams for both non-rephasing and rephasing contributions at all the possible phase-matching directions are provided in Supplementary Fig. 1. In the absence of EPC, the electronic transitions do not depend on the phonon occupation and vice versa, that is, μ_g0,e0 = μ_g1,e1 and μ_g0,g1 = μ_e0,e1. We can then readily convince ourselves that the product of the matrix elements in equation (6) is identical for all three diagrams. Evaluating the sum over the four remaining terms yields

which equals zero. The response function considering all the rephasing contributions also undergoes total cancellation in a similar manner. This means the response function and, therefore, the signal vanishes when the phonons do not couple to the electronic transition. This cancellation has the same origin as the vanishing of the cross-peaks in the absence of, for example, intermolecular interactions in 2D infrared (IR) spectroscopy²⁸.

The presence of EPC, on the other hand, has two important consequences for the 2D EPC spectrum. First, the total spectroscopic signal is non-zero and is, in fact, a function of EPC strength as it depends on the difference in the matrix elements \({\mu }_{e0,g0}-{\mu }_{e1,g1}\propto {M}_{{\bf{k}},\lambda }^{2}{{\rm{e}}}^{-{M}_{{\bf{k}},\lambda }^{2}/2}\), with the effective EPC strength M_k,λ ≡ M_k,e,λ – M_k,g,λ. Therefore, although the intraband phonon transitions are independent of the EPC, the interband transitions only depend on the effective EPC strength, M_k,λ. The fact that the 2D EPC spectrum is sensitive to the difference in the EPC strength from the conduction and valence bands is a common feature of experimental methods based on optical transitions²⁹. Second, part of the spectral weight of the interband transitions is transferred to higher phonon numbers and, therefore, additional pathways start to contribute (see the ‘Discussion’ section). It is noteworthy that EPC can also result in a phonon frequency shift with band index, which is not included here. However, we believe our results below to remain qualitatively true in that case too, as such effects will similarly destroy the perfect cancellation explained above.

For the experimental demonstration, the pulse characteristics are chosen to match the energy scales of MAPI at room temperature. The THz pulse with a spectral range up to 60 meV contains sufficient coverage to excite all inorganic-sublattice-based phonon modes reported in MAPI (Supplementary Fig. 2a)³⁰. To gain access to the broad range of electronic transition energies above the bandgap, we tune the NB pulse from 1,300 nm to 2,000 nm and produce a BB pulse with a spectral width of 1,000–1,250 nm (Supplementary Fig. 2b). Both these photon energies are away from any elementary excitations in MAPI. These two optical (that is, BB and NB) pulses collinearly propagate through the sample simultaneously, and the delay time between the optical pulse pair and the THz pulse, τ, is scanned. The phase of the signal field E_EPC(t, τ) is resolved by interference with a local oscillator (LO) field (Supplementary Fig. 2c). The LO field is generated after the sample, and both fields are dispersed at a spectrograph and detected by a camera. The thus-measured intensity I_EPC(ω_t, τ) containing the squared magnitude of the total field, E_tot(t, τ) = E_EPC(t, τ) + E_LO(t) + c.c., is then Fourier transformed along the delay time τ, producing the frequency-domain 2D EPC spectrum I_EPC(ω_t, ω_τ) (Supplementary Fig. 3 and Methods).

From the measured I_EPC(ω_t, ω_τ) spectrum, we extract the 2D EPC response function S_EPC(ω_t, ω_τ) by the following analysis protocol. First, we isolate the desired phase-matching condition. To do so, we introduce an additional time delay between the LO and the signal field that contributes to opposite responses between the rephasing and non-rephasing pathways, and separate it temporally (equations (15) and (16), Supplementary Figs. 4 and 5 and Supplementary Text 1.1). Second, we correct the relative phase differences among three incident pulses at the sample by using a SiN membrane as a reference sample. We, thus, obtain the real and imaginary spectra for both non-rephasing and rephasing contributions, separately (Supplementary Fig. 6 and Supplementary Text 1.2). In fact, such a full separation could be useful for analysing the lineshape of isolated resonances³¹. Here we focus on the slightly larger non-rephasing contribution. Finally, we normalize the isolated 2D spectrum \({I}_{{\rm{EPC}},{{\bf{k}}}_{{\rm{s}}}}({\omega }_{t},{\omega }_{\tau })\) by the spectra of E_THz(ω_τ) and E_LO(ω_t). The LO spectrum is, in fact, the upconverted BB spectrum by the NB frequency (Supplementary Fig. 2c and Methods).

The non-rephasing 2D EPC spectrum of the sample obtained from this analysis, S_EPC(ω_t, ω_τ), is shown in Fig. 2 for a range of electron energies (~1.6–2.2 eV). Evidently, both 1-THz and 2-THz modes couple to electrons with energy above the bandgap to some extent as they contribute to the 2D EPC signal. The signals from both modes monotonically increase as the electron energy increases up to 0.6 eV above the bandgap, with the contribution of the 1-THz mode consistently larger than that of the 2-THz mode.

For quantifying the EPC strength at a given phonon mode λ and electronic energy ε_k, it is important to remember that the 2D EPC spectrum S_EPC(ω_t, ω_τ) (Fig. 2c) depends on four transition matrix elements, namely, two phononic and two electronic transition elements (for example, Fig. 3a (top right) shows μ_g1,g0 and μ_em,e(m+1) for the phononic transitions and μ_e(m+1),g1 and μ_g0,em for the electronic transitions). Moreover, the 2D EPC signal field contains the phonon and electron densities of states. Therefore, to isolate the pure contribution of the EPC strength to the 2D EPC response function for a given phonon mode λ, we normalize the measured response function S_EPC(ω_t, ω_τ) (Fig. 2c) by both absorbance spectra along the electronic (Fig. 2a) and phononic (Fig. 2b) energies, resulting in \({S}_{{\rm{EPC,norm}}}^{\lambda }({\omega }_{t})\) (Fig. 2d (closed circles); the theoretical analysis is given in the ‘Discussion’ section). We note here that our approach is applicable to generic nuclear displacements, but the quantitative comparison of EPC strengths can only be performed for phonon modes with non-zero oscillator strengths. The result indicates a consistently stronger coupling of electrons in this energy range to the 1-THz mode than to the 2-THz mode. Both normalized signals monotonically increase from zero near the bandgap and the growth slows down as the electron energy approaches 2 eV.

Using our theoretical model, we can relate the measured 2D EPC signal to the EPC strength M_k,λ. Our starting point is an expression for the induced polarization P_EPC,k(ω_t, ω_τ) at fixed electronic momentum k (Methods). Because transitions can involve multiple phonon excitations in the presence of EPC, the signal field is emitted at a series of discrete frequencies spaced by the phonon frequency, that is, ε_k + mω_λ, where ε_k = ε_e,k – ε_g,k is the momentum-dependent electronic transition energy. Considering the non-rephasing contribution (that is, ω_τ > 0 and ω_t > 0), such emissions take place via four different pathways (Fig. 3a). We evaluate these diagrams explicitly in Supplementary Text 1.4 and the resulting induced polarization is given by \({P}_{{\rm{EPC}},{\bf{k}}}({\omega }_{t},{\omega }_{\tau })\propto {\sum }_{m}\delta\)\(({\epsilon }_{{\bf{k}}}+m{\omega }_{\lambda }-{\omega }_{t}){C}_{m,{\bf{k}}}{{\rm{e}}}^{-{M}_{{\bf{k}},\lambda }^{2}}{({M}_{{\bf{k}},\lambda })}^{2m}({M}_{{\bf{k}},\lambda }^{2}-m)/m!\) (Fig. 3b). Here C_m,k is a constant that weakly depends on k (Methods). The polarization changes sign as a function of the detection frequency ω_t, crossing zero at a frequency corresponding to the typical number of additional phonon excitations, \(m={M}_{{\bf{k}},\lambda }^{2}\). Interestingly, the sum over all frequencies, ∫dω_tP_EPC,k(ω_t), vanishes regardless of the EPC strength. This cancellation originates from the destructively interfering pathways discussed above. However, in the presence of EPC, the negative contributions to the induced polarization acquire a stronger blueshift than the positive ones, leading to an imperfect cancellation. This leads to the characteristic dispersive lineshape (Fig. 3b) expected for the 2D spectra of discrete levels³¹. Our measured spectrum, however, has a different lineshape as the electronic excitations form a continuum and the full signal involves a convolution of the density of states with P_EPC,k.

The total induced polarization at a particular frequency ω_t is given by a sum over all electronic momenta:

where ν is the density of states corresponding to the transition energy ε_k. The density of states is non-zero for arguments above some minimal energy ε_Γ, which is the smallest energy separation between the electronic ground and excited states with the same momentum. Equation (8) can be used to numerically evaluate the 2D EPC spectrum for a given momentum-dependent coupling strength M_k,λ. To determine the EPC strength from a measured spectrum, we use an analytical approximation. Specifically, we assume the coupling strength \({M}_{{\varepsilon }_{{\bf{k}}},\lambda }\equiv {M}_{{\bf{k}},\lambda }\) to vary little with ε_k on the scale of the phonon frequency ω_λ, and we obtain

where f(ω_τ) = Γω_λ/[(ω_τ – ω_λ)² + Γ²] is the lineshape along ω_τ and C(ε_k) is given in the Methods. This equation can be inverted and we find

where \({S}_{{\rm{EPC,norm}}}^{\lambda }({\omega }_{t})\) corresponds to the normalized 2D EPC spectrum for a given phonon mode (Fig. 2d) and we have only omitted prefactors that do not depend on ω_λ and ω_t (Methods).

To determine the phonon-mode- and electron-energy-specific coupling strength from this expression, we fit the data in Fig. 2d with a polynomial function for each mode with the same zero crossing and use the transition energy ε_k = ε_Γ + Δcos(ka). Here the bandgap ε_Γ = 1.68 eV has been obtained from the zero crossing of the fit functions, and the bandwidth of MAPI, Δ = 1.2 eV, was taken from another work³². The resulting effective EPC strength \({M}_{{\varepsilon }_{{\bf{k}}},\lambda }\) is shown in Fig. 2e. Note that the intensity of the signal drops to zero at the band edge due to the vanishing density of electronic states below the gap. Independently, as a result of our fit, we find that the EPC strength for both modes approaches zero at the band edge.

With the thus-obtained EPC strength \({M}_{{\varepsilon }_{{\bf{k}}},\lambda }\), we can then reproduce the full 2D EPC spectrum. To model the two pronounced phonon modes (namely, the ~1-THz and ~2-THz modes), we choose the frequency, width and oscillator strength to match the experimental THz spectrum shown in Fig. 2b. With this, we obtain the 2D EPC spectrum shown in Fig. 3c by numerically evaluating equation (8). The calculated normalized spectrum (Fig. 2d, solid lines) is in good agreement with the experimental data.

The phonon-mode-resolved EPC strength as a function of electron energy (Fig. 2e) is the main result of this study. We observe a roughly linear drop of the EPC at the band edge ε_k = ε_Γ for both modes, which corresponds to a vanishing of the EPC strength as ~k² at k = 0. Such a strong dependence on the electron wavevector is in marked contrast to a standard Holstein coupling term, which assumes a coupling of phonons to the on-site density of electrons and does not vary with k. The behaviour shown in Fig. 2e suggests instead an EPC that affects the kinetic energy rather than the potential energy of electrons.

To gain further insight, we consider the SSH model as a paradigmatic model with a k-dependent EPC. In Supplementary Text 1.5.1, we show that the EPC strength for a q = 0 optical phonon mode in the SSH model indeed vanishes at k = 0 as M_k,q=0 ∝ k. The SSH model, however, is not directly applicable to our case as it describes longitudinal phonons. We have, therefore, developed another toy model for transverse phonons with an EPC that vanishes at q = k = 0 (Supplementary Text 1.5.2). Although the microscopic coupling mechanism in both models is rather distinct, we show that the strong electron momentum dependence of the EPC in the latter can be understood as an effective renormalization of the hopping strengths in the spirit of the SSH model. In contrast to the SSH model, however, the EPC strength of transverse phonons vanishes quadratically, M_k,q=0 ∝ k², as required by symmetry and in agreement with our experimental observations. Moreover, the characteristic behaviour at q = k = 0 of these two toy models agrees with the more realistic electron–phonon models for longitudinal and transverse modes discussed elsewhere^4,33 in the context of cuprates. We conclude that measuring the electron momentum dependence of the EPC can be a powerful indicator of the microscopic structure of coupled electron–phonon systems and help distinguish the couplings of phonons to the kinetic and potential energies of electrons.

Since we are able to directly measure and extract the phonon-mode- and electron-energy-resolved EPC strength at a fixed condition (for example, temperature), we are now in a position to further explore how the EPC evolves across a phase transition. MAPI is known to undergo a structural phase transition from a tetragonal (I4/mcm) to an orthorhombic (Pnma) phase below T_c ≈ 160 K (ref. ³⁴) and to have optical and electronic properties that are highly sensitive to temperature changes^35,36,37. Due to the symmetry reduction in the low-temperature phase, the degeneracy of both 1- and 2-THz modes in the high-temperature phase is lifted and both modes split into two modes below T_c (ref. ³⁸). Both splittings occur in such a way that additional modes with slightly lower frequencies (denoted by ~0.7 THz and ~1.5 THz; Fig. 4a) appear besides the original modes (denoted by ~1 THz and ~2 THz; Fig. 4a). Note that all the mode frequencies increase slightly with temperature. In 2D EPC (Fig. 4a; the 2D EPC spectra are shown in Supplementary Fig. 10), the EPC of each of these modes can be clearly distinguished. However, the original modes keep dominating the coupling strengths (with the 1-THz mode having a consistently stronger EPC than the 2-THz mode). Given the mode degeneracy above T_c, the combined contributions, for example, from the ~1.5-THz and ~2-THz phonons to the EPC strength at ~2 THz is overall reduced across the transition from orthorhombic (below-T_c) to tetragonal (above-T_c) phase.

Finally, another intriguing attribute of the phonon-mode- and electron-energy-specific EPC strength arises when we choose different polarization combinations of the three incoming pulses. With the given laboratory axis X and the propagation direction Z, the polarization directions of THz, NB, BB and signal fields are written in order, for example, XYYX (Supplementary Fig. 11). As evident from Fig. 4b, the black (all parallel) and red (cross-polarization between the THz and the two optical pulses) curves reveal a contrast between the coupling of the 1-THz mode and that of the 2-THz mode. With the cross-polarization condition, the signal of the 1-THz mode is suppressed largely compared with that of the 2-THz mode, which highlights a clear anisotropic nature of EPC only for the 1-THz case at ~2 eV. Further theoretical studies will be beneficial to gain microscopic insights into the geometrical nature of the phonon-mode-specific EPC and its sensitivity to the electronic momentum in terms of relevant symmetries.

In conclusion, we have demonstrated how 2D spectroscopy can be used to obtain detailed information about the strength and microscopic origin of EPC in solids. This approach extends conventional 2D spectroscopy of discrete modes like molecular excitations or excitons^28,31,39 to band structures in solids with a continuous spectrum, which enables us to determine the electron momentum dependence of the EPC coupling strength.

Compared with established experimental techniques that are sensitive to EPC, 2D EPC spectroscopy has certain limitations and unique strengths. First, 2D EPC spectroscopy probes the differences in EPC strengths for different electronic bands. This is tied to the fact that our technique directly probes the bare EPC matrix elements rather than indirect consequences of the coupling like the lifetime of electronic quasiparticles, which require an average over the final electronic scattering states. Because 2D EPC spectroscopy directly probes the correlated electron–phonon excitations, it remains unaffected by spurious effects such as electron–electron interactions, which may lead to ambiguity in all-electronic measurements. The band-resolved EPC strength, therefore, can be obtained by comparing 2D EPC with other experiments or ab initio calculations. Moreover, the EPC difference is itself a relevant physical observable, for instance, it determines the energy renormalization and non-radiative lifetime of excitons due to phonons⁴⁰. Second, 2D EPC spectroscopy probes optically active phonons at zero momentum. This may be particularly relevant for determining the EPC strength near instabilities driven by q = 0 phonons, such as structural phase transitions, charge-density wave systems or excitonic insulators (see, for example, ref. ⁴¹). Third, our 2D EPC spectroscopy measures the EPC strength as a function of electron energy. For a known band structure, we can, therefore, extract the electron momentum dependence of the EPC strength averaged over isoenergy surfaces.

Most importantly, the unique combination of electron momentum (k) and phonon-mode resolution of our technique enables us to detect signatures of the microscopic origin of the EPC. The strong k dependence along with the vanishing of the EPC strength at k = 0, which we identify in MAPI, are powerful evidence for a coupling of phonons to the kinetic energy rather than potential energy. This makes our technique ideally suited for experimentally distinguishing non-local SSH-type couplings from local Holstein-type couplings. Incidentally, the contrast between both models is the most pronounced at phonon momentum q = 0, where our technique operates (Supplementary Text 1.5). This is tied to the fact that at finite phonon momentum, the interference of electron and lattice waves leads to a partial cancellation of k dependence.

Hence, 2D EPC spectroscopy yields complementary information compared with traditional techniques sensitive to EPC such as angle-resolved photoelectron spectroscopy, heat capacity or transport measurements, and can, therefore, play an important role in improving our understanding of EPC. Our work paves the way for studying the evolution of mode-specific EPC under changes in thermodynamic or mechanical conditions, as well as phonon-mediated light-induced ultrafast control of condensed materials. The combination of material-specific calculation and this experimental approach is expected to provide deeper insights into the nature of EPCs.

Perovskite film preparation

Methylammonium iodide was purchased from Greatcell Energy. Lead acetate (98% purity) was purchased from Tokyo Chemical Industry. Anhydrous dimethylformamide (>99.8% purity) was purchased from Sigma-Aldrich. Substrates (0.3-mm-thick glass) were cleaned by sonication in water with soap for 30 min, rinsed afterwards with deionized water three times, and again sonicated in ethanol and acetone (both for 30 min). Afterwards, the substrates were dried under a nitrogen flow and subjected to ultraviolet–ozone treatment (FHR UVOH 150 LAB, 250 W) for 20 min with an oxygen-feeding rate of 1 l min^–1 right before spin coating. Film preparation was carried out in a nitrogen-purged glovebox. The precursor solution consists of 477 mg (3 mmol) methylammonium iodide and 325.3 mg lead acetate (1 mol) in 1 ml dimethylformamide. The film preparation was done by depositing 50 μl of the precursor solution onto a cleaned fused silica substrate and spin coating it at 4,000 rpm (ramp, ±1,000 rpm s^–1) for 1 min. Afterwards, the film was left to dry at room temperature for 10 min and subsequently annealed on a 100 °C hotplate^42,43. The films were stored in the glovebox before their use.

Experimental details

The output of an amplified femtosecond Ti:sapphire laser (1 kHz repetition rate, 50-fs pulse duration) is split into ~1 mJ, 1 mJ and 0.1 mJ for generating the THz, NB and BB near-infrared (NIR) pulses, respectively. The broadband THz pulses are generated via a femtosecond two-colour laser focusing in long air-plasma filaments⁴⁴. The NB NIR pulses with a bandwidth of ~20 cm⁻¹ are obtained by sending the output of a commercial optical parametric amplifier into a custom-built pulse shaper, which spatially selects spectral components by the combination of a ruled reflective diffraction grating and a slit. The BB NIR pulses are produced via white light generation at 4-mm-thick yttrium aluminium garnet crystal⁴⁵ (Supplementary Fig. 2).

We combine the NB and BB pulses with a dichroic beam combiner in such a way that they collinearly propagate through the sample. The pulse energies of NB and BB pulses at the sample position are of ~1–2 μJ and ~0.02 μJ, respectively, and are focused by an off-axis parabolic mirror with a focal length of 15 cm. The THz pulses are focused by another off-axis parabolic mirror with a focal length of 10 cm. To measure the signal in a phase-resolved fashion, we adopt the heterodyne detection. The LO pulses are, therefore, generated after the sample position by sum-frequency generation of NB and BB pulses at a 10-μm-thick β-barium borate crystal. Since we set NB and BB to overlap in time, when they pass through the sample and reach the β-barium borate crystal, new pulses with a frequency of ω_NB + ω_BB are generated. The thus-generated LO and signal fields are collimated via a CaF₂ lens and collinearly propagate through an analyser polarizer. The intensity of the LO field interfering with and without the signal field are spatially separated by a vibrating motorized mirror (galvanometer), spectrally resolved and imaged separately onto a charge-coupled device camera. The area of the setup around the THz-beam pathway was enclosed in a box and purged with dry nitrogen. The sample is kept in a cryostat.

Data analysis

The 2D EPC spectroscopy uses the third-order nonlinear optical process, and the generated 2D EPC signal field E_EPC is proportional to the induced polarization to the third order of the electric field of the incident light:

where t is the time variable of the signal field E_EPC and τ is a time delay of the two IR pulses relative to the THz pulse. The electric fields E_THz(ω₁), E_NB(ω₂) and E_BB(ω₃) correspond to the THz, NB and BB pulses at frequencies ω₁, ω₂ and ω₃, respectively. Here ω_s ≡ ω₁ ± ω₂ ± ω₃ and k_s ≡ k₁ ± k₂ ± k₃. The nonlinear response function S_EPC is expressed in the time domain in terms of correlation functions with dipole moment operators μ:

where θ(t) is the Heaviside function and ρ(–∞) is the equilibrium-reduced density matrix for the system eigenstates at t₀ (Fig. 1a)²⁷. The relevant diagrams for equation (12) are shown in Fig. 1c. Each term of the commutators indicates a distinctive Liouville-space pathway at a given phase-matching direction k_s. Then, the raw measured spectrum in the frequency domain, I_EPC(ω_t, ω_τ), contains convoluted contributions from each of the incoming fields themselves, the interference of signals from different phase-matching directions and the relative phase shift of the incoming fields at the interface, apart from the response function S_EPC(ω_t, ω_τ) that needs to be isolated.

Technically, the combined signal field and LO field are guided to the spectrometer (that is, heterodyne detection), where they are dispersed by the grating. This dispersion is mathematically described by a Fourier transform of the field, E_tot(t, τ) = E_EPC(t, τ) + E_LO(t), over time t. The signal intensity measured by the square-law detector (charge-coupled device camera) for the detection frequency ω_t ≥ 0 is equivalent to \({I}_{{\rm{EPC}}}({\omega }_{t},\tau )\propto {\int}_{T}{\rm{d}}t{[{E}_{{\rm{tot}}}({\omega }_{t},\tau ){{\rm{e}}}^{-{\rm{i}}{\omega }_{t}t}+{E}_{{\rm{tot}}}(-{\omega }_{t},\tau ){{\rm{e}}}^{{\rm{i}}{\omega }_{t}t}]}^{2}\) with the accumulation time period T. By squaring the term under the integral, neglecting terms oscillating at very high frequency (~2ω_t) and subtracting the homodyne terms of the form E_EPC(ω_t, τ)E_EPC(–ω_t, τ) and E_LO(ω_t, τ)E_LO(–ω_t, τ), we obtain the heterodyne-detected part of the measured signal (Supplementary Fig. 3a). Through further Fourier transform over the time delay τ of the THz pulse relative to the NB and BB pulses, the measured (heterodyne-detected) 2D EPC spectrum is given by

Meanwhile, from the Fourier transform of equation (11) along t, E_EPC(ω_t, ω_τ) can be written as

Since we use an NB NIR pulse with frequency ω_NB in the experiment, we approximate the NB NIR spectrum by a delta function: E_NB(ω₃) = δ(ω_NB – ω₃) + δ(ω_NB + ω₃). By performing the integration of equation (14) and use the experiment condition that E_BB(ω_t + ω_NB – ω_τ) = 0 (ω_τ ≪ ω_t + ω_NB), we obtain

Similarly, for E_EPC(–ω_t, ω_τ), we obtain

By substituting equations (16) and (15) into equation (13), we obtain the 2D EPC spectrum:

Due to the collinear geometry of our experiment, the two terms in equation (17), that is, signals from different phase-matching conditions (non-rephasing pathways (ω_τ > 0 and ω_t > 0) and rephasing pathways (ω_τ > 0 and ω_t < 0)) are not spatially separated and spectrally very close (ω_τ ≪ ω_t). This co-propagation results in the interference of the two signals in the raw 2D spectrum (Supplementary Fig. 3b). To separate these two types of contribution, we use the fact that the two terms in equation (17) respond to the time delay t′ between the LO and signal field in an opposite way:

where \({I}_{{\rm{EPC}}}^{\,{\rm{nr}}}\) (\({I}_{{\rm{EPC}}}^{\,{\rm{r}}}\)) indicates the signal intensity contributed from the non-rephasing (rephasing) pathway. Hence, in the inverse Fourier transform of the 2D spectrum along ω_t, contributions from \({I}_{{\rm{EPC}}}^{{\rm{nr}}}\) and \({I}_{{\rm{EPC}}}^{{\rm{r}}}\) are temporally separated and one could select one of the two phase-matching directions (Fig. 2c)³⁹.

Once we isolate the 2D spectrum of non-rephasing contributions, we could now utilize from equation (18) that the 2D response function S_EPC(ω_t, ω_τ) can be obtained by dividing the measured 2D spectrum \({I}_{{\rm{EPC}}}^{{\rm{r}}}({\omega }_{t},{\omega }_{\tau })\) directly by the incident field spectra, E_THz(ω_τ) and E_LO(ω_t) (each spectrum is shown in Supplementary Fig. 2a,c). Note that we do not scan the time between the NB and BB pulses; therefore, the spectrum of the BB pulse is not needed.

Finally, the relative phase differences of the incident fields at the sample are compensated by using the 2D EPC spectrum of a reference sample measured in the experimentally identical condition. A reference sample must have no resonance; therefore, the response function spectrum \({S}_{{\rm{EPC}}}^{{\rm{ref}}}({\omega }_{t},{\omega }_{\tau })\) of the reference sample is a real constant over the spectral range (which ω_t and ω_τ cover). Then, the phase-shifted incident fields, that is, \({\tilde{E}}_{{\rm{THz}}}({\omega }_{\tau })={E}_{{\rm{THz}}}({\omega }_{\tau })\exp {\rm{i}}{\omega }_{\tau }{t}^{{\prime}{\prime}}\) and so on, where t″ is relative phase difference, can be cancelled out by⁴⁶

In this experiment, we use a 1-μm-thick SiN membrane as the reference sample (Supplementary Figs. 7–9).

Theoretical model for 2D EPC spectroscopy

Here we theoretically calculate the 2D spectrum for the generic electron–phonon model in equation (2). In Supplementary Text 1.3, we calculate the induced polarization P_EPC,tot(ω_t) by considering the contributions from all possible diagrams and intermediate states (Supplementary Fig. 12). We find that the contributions from diagrams \({R}_{2}^{* }\) and \({R}_{3}^{* }\) to the polarization cancel for our specific pulse sequence, and we only need to evaluate the contribution from diagram R₄. All the possible pathways relevant for emission at frequency ε_k + mω_λ of the non-rephasing condition contain the sequence of states |g, 0〉→|g, 1〉→|i〉→|e_k, m〉→|g, 0〉 and we have to sum over all the intermediate states |i〉 (Fig. 3a). Since the transition dipole matrix element of the conduction and valence bands is non-zero, the state |i〉 must also be in the conduction or valence band; otherwise, one of the transition matrix elements μ_iν,gν′ or μ_iν,eν′ would be zero by symmetry. The detailed evaluation of all the pathways is presented in Supplementary Text 1.4. Summing over all the intermediate states, we obtain \({P}_{{\rm{EPC}},{\bf{k}}}({\omega }_{t}={\varepsilon }_{{\bf{k}}}+m{\omega }_{\lambda })\propto {{\rm{e}}}^{-{M}_{{\bf{k}},\lambda }^{2}}{({M}_{{\bf{k}},\lambda })}^{2m}({M}_{{\bf{k}},\lambda }^{2}-m)/m!\), for a given electron momentum k and phonon number m of mode λ. The polarization at a detection frequency ω_t is then given by summing over all phonon numbers and all electronic momenta of the final state as

Here ν(ε) is the electronic density of states, Γ_λ is the intrinsic broadening of the phonon mode λ, ω_τ is the Fourier transform of the delay time τ between the first pulse and the NB and BB pulses, and C_m,k is given by

where ω_2/3 and σ_2/3 are the central frequencies and time-domain widths of the two IR pulses, respectively, which are assumed to be Gaussian. Writing this expression, we have assumed that the frequencies of the incoming pulses are chosen to be resonant with the respective transitions. We can perform the sum over m analytically, if we ignore the weak dependence of C_m,k on m, that is, C_m,k ≃ C_0,k, and we obtain

where δ′(ω) ≃ [δ(ω + ω_λ) – δ(ω)]/ω_λ denotes the derivative of a δ function. Substituting this into equation (22), we find

In this expression, only ν(ε_k) and \({M}_{{\bf{k}},\lambda }^{2}\) have a strong dependence on k, and, thus, we can further simplify the result as equation (9) with the definition

To determine the EPC strength from the experimental data, we fit the normalized 2D EPC spectrum for a given phonon mode, \({S}_{{\rm{EPC,norm}}}^{\lambda }({\omega }_{t})\) (Fig. 2d). In our theoretical model, the normalized spectrum \({S}_{{\rm{EPC,norm}}}^{\lambda }({\omega }_{t})\) corresponds to the intensity I_EPC ∝ ω_tP_EPC,tot divided by the THz absorbance \(\propto {\omega }_{\lambda }{\vert \langle {\chi }_{g,\lambda }^{0}| {\mu }_{p}| {\chi }_{g,\lambda }^{1}\rangle \vert }^{2}\varGamma /[{({\omega }_{\tau }-{\omega }_{\lambda })}^{2}+{\varGamma }^{2}]\) and IR absorbance ∝ω_tν(ω_t)|〈e₀|μ_e|g〉|² (ref. ²⁷). We then use equation (10) to determine the EPC strengths. We perform the integration by assuming the density of states to be a power law as ν(ω_t) ∝ (ω_t – ε_Γ)^α and by fitting the experimental data for \({S}_{{\rm{EPC,norm}}}^{\lambda }({\omega }_{t})\) in Fig. 2d with a polynomial. Specifically, we fit the spectrum of the 1-THz phonon mode by a quadratic function and that of the 2-THz mode by a linear function, assuming both functions cross zero at the gap edge, which we find to be ε_Γ = 1.68 eV. This allows us to obtain \({M}_{{\varepsilon }_{{\bf{k}}},\lambda }\) (Fig. 2e). We verify that the coupling strength obtained from the analytical expression in equation (10) indeed reproduces the experimental data by plotting the full numerical result based on equations (8) and (22) in Fig. 2d. The broadening of the phonon modes are obtained by fitting the THz absorbance in Fig. 2c with a Lorentz oscillator model of absorbance³⁸, from which we find Γ₁ = 0.11 THz and Γ₂ = 0.27 THz. The oscillator strength of the 2-THz mode is found to be approximately four times higher than that of the 1-THz mode.