Accessing the anisotropic non-thermal phonon populations in black phosphorus

We combine femtosecond electron diffuse scattering experiments and first-principles calculations of the coupled electron-phonon dynamics to provide a detailed momentum-resolved picture of the ultrafast lattice thermalization in a thin film of black phosphorus. The measurements reveal the emergence of highly anisotropic non-thermal phonon populations which persist for several picoseconds following excitation of the electrons with a light pulse. Combining ultrafast dynamics simulations based on the time-dependent Boltzmann formalism and calculations of the structure factor, we reproduce the experimental data and identify the vibrational modes primarily responsible for the carrier relaxation via electron-phonon coupling and the subsequent lattice thermalization via phonon-phonon scattering. In particular, we attribute the non-equilibrium lattice dynamics of black phosphorus to highly-anisotropic phonon-assisted scattering processes, which are primarily mediated by high-energy optical phonons. Our approach paves the way towards unravelling and controlling microscopic energy-flow pathways in two-dimensional materials and van der Waals heterostructures, and may also be extended to other non-equilibrium phenomena involving coupled electron-phonon dynamics such as superconductivity, phase transitions or polaron physics.

phonon-phonon scattering. In particular, we attribute the non-equilibrium lattice dynamics of black phosphorus to highly-anisotropic phonon-assisted scattering processes, which are primarily mediated by high-energy optical phonons. Our approach paves the way towards unravelling and controlling microscopic energy-flow pathways in twodimensional materials and van der Waals heterostructures, and may also be extended to other non-equilibrium phenomena involving coupled electron-phonon dynamics such as superconductivity, phase transitions or polaron physics.
Black phosphorus (BP) exhibits a tunable band gap in the mid-IR, 1-3 high carrier mobilities, 4-6 and a layered crystal structure. These features make it a versatile platform to explore novel device concepts, such as field-effect transistors, saturable absorbers, and polarizationsensitive photodetectors. 2,4,5,[7][8][9] The pronounced crystal structure anisotropy of BP further underpins the emergence of strikingly anisotropic macroscopic properties, as exemplified by its thermal [10][11][12] and electrical conductivities, 1,5,13,14 as well as its optical response. 5,[15][16][17] Since practical applications based on these properties invariably exploit non-equilibrium states of either the lattice or hot carriers, it is desirable to attain a detailed understanding of the ultrafast dynamics of electronic and vibrational degrees of freedom in BP. Following photo-excitation, hot carriers relax to the band edges by transferring their excess energy to the lattice via the emission of phonons, which triggers coupled carrier-lattice non-equilibrium dynamics. Optical and photoemission spectroscopies have been employed extensively to investigate carrier-phonon scattering channels and their influence on the carrier dynamics in BP. 14,[18][19][20][21][22][23][24][25] While these techniques provide direct information on the electrons, the nonequilibrium dynamics of the lattice can only be inferred indirectly through its effects on the electronic structure. Femtosecond electron diffuse scattering (FEDS), conversely, is an ideal technique to circumvent these limitations and complement optical and photoemission spectroscopies. FEDS provides direct access to lattice dynamics and electron-phonon scattering processes with time and momentum resolution. [26][27][28] Owing to its sensitivity to both electron-phonon and phonon-phonon scattering processes in reciprocal space, FEDS is thus well-suited to establish a microscopic picture of the energy flow between hot electrons and the BP lattice.
Here, we combine FEDS experiments with ab-initio calculations to investigate the coupled electron-phonon dynamics in BP following photo-excitation. Our study reveals that strongly anisotropic non-thermal phonon populations are established throughout the first picoseconds of the dynamics. A regime of thermal equilibrium is only re-established by the ensuing anharmonic decay pathways (phonon-phonon coupling) on timescales of the order of 50-100 ps. To unravel the origin of the non-equilibrium lattice dynamics and its signatures in FEDS experiments, we conduct first-principles calculations of the coupled electron-phonon dynamics based on the time-dependent Boltzmann formalism, whereby electron-phonon and phonon-phonon scattering processes are explicitly accounted for. Calculations of the structure factor further enable a direct comparison with the experimental data. We show that the relaxation of excited electrons and holes is governed by the emission of high-energy optical phonons within a restricted region of the Brillouin zone, and it is responsible for the anisotropic non-thermal phonon populations revealed by FEDS.

Results and discussion
The layered orthorhombic crystal structure of BP is illustrated in Figure 1 Bright high-intensity features arise for transferred momenta matching the reciprocal lattice vectors G, according to Bragg's law for the interference condition. These measurements are in good agreement with previous TEM experiments. 31 Besides the pronounced anisotropy of the BP crystal lattice, marked by different structural motifs along the armchair and zigzag

Femtosecond electron diffuse scattering measurements
To investigate the non-equilibrium lattice dynamics of BP with momentum and time resolutions, we perform FEDS measurements on a free-standing thin film of BP. The sample has an estimated thickness of 39 ± 5 nm and it has been obtained by mechanical exfoliation of a bulk crystal. In FEDS, a laser pulse is employed to drive the system into an excited electronic state. After a controllable time delay t, the sample is probed by an electron pulse which diffracts off the lattice. The distribution of the diffracted electrons generated by this procedure provides a direct probe of the non-equilibrium dynamics of the lattice in reciprocal space. 32 A schematic illustration of the experiment is reported in Figure 1 Figure 2 (b) for the A (circles) and X (pentagons) points in the BZ. As diffuse scattering occurs primarily through phonon-induced scattering processes, the signal measured at a given point q in the BZ reflects the phonon population with the same momentum. [26][27][28][34][35][36] The red curve in Figure 2

Theoretical modeling of non-equilibrium lattice dynamics
To gain further insight into the non-equilibrium dynamics of the lattice, we perform firstprinciples calculations of the coupled electron-phonon dynamics of BP based on the timedependent Boltzmann equations: 37 Here, ∂ t = ∂/∂t, f nk denotes the electron distribution function for electron band index n and electron momentum k, and n qν the phonon distribution function for wavevector q and branch index ν. Equations (1) and (2)  is accounted for by the collision integral Γ ep nk (Γ pe qν ), whereas the phonon-phonon scattering due to lattice anharmonicities is accounted for via Γ pp qν . In short, Eqs. (1) and (2)  Owing to the absence of local minima in conduction band along the armchair direction (i.e., Γ-X and Z-Q), the photo-excited electrons are constrained to occupy states with crystal momenta along the zigzag direction, i.e., where the available local minima are located (arrows in Figure 1 (d)). This scenario is illustrated by highly-anisotropic electronic occupations f 0 nk in the conduction band, reported in Figure 1 Figure 4(a)) throughout the first 50 ps of the dynamics. For momenta around Γ and A, the temperature reaches a maximum at 1.7 and 2.3 ps, respectively, whereas no maximum is observed around X. These timescales indicate the time required for the electrons to transfer energy to the lattice via electron-phonon scattering. The good agreement with the experimental time constant of 1.7 ps extracted from the rise of the FEDS intensity at A (Figure 2(b)), suggests that transient changes of the FEDS intensities for timescales smaller than 2 ps reflect primarily the energy transfer from the electrons to the lattice driven by the electron-phonon coupling.
In Figure 4 To inspect directly the influence of the non-equilibrium lattice dynamics on the scattering intensity probed in the FEDS experiments, we conduct first-principles calculations of the structure factor by explicitly accounting for the influence of electron-phonon interactions and anisotropic population of the vibrational modes in the unit cell. Specifically, we perform computations of the all-phonon structure factor I all (Q, T ). 39,40 Indeed we find that taking into account multi-phonon effects is essential for an accurate reproduction of the experimentally observed diffraction patterns of BP seen in Figures 3(a-c). The expression for I all (Q, T ) reads: Here N p is the number of q-points used to sample the first Brillouin Zone, f κ (Q) denotes the scattering amplitude of atom κ, W κκ (Q, T ) is the Debye-Waller factor, τ κ represents the atomic positions and R p defines the position vector of unit cell p contained in a Born-von Kármán supercell. The phononic factor, e P p,κκ (Q,T ) , includes all orders of phonon processes and its exponent is given by: where M κ and M 0 are the atomic and reference masses, and the phonons are described by the eigenmodes e κ,ν (q) and frequencies ω qν . A key quantity entering the equation of the structure factor is the mean-squared displacement of the atoms due to mode qν, defined as The time-dependence of the all-phonon structure factor is encoded in u 2 qν T , which is directly related to phonon populations n qν (T ). To account for the influence of the non-equilibrium lattice dynamics on the FEDS intensity, we evaluated Eq.
(3) at each time snapshot by populating phonons according to the vibrational temperatures obtained from the solution of the time-dependent Boltzmann equation (Figure 4).
The calculated (non-equilibrium) all-phonon structure factor is shown in Figure 3  diamond-shaped diffraction pattern that characterises the return to thermal equilibrium.
These findings enable us to establish the picture sketched in Figure 5 for the non-  Figure 4(b-d)). Distinctive fingerprints of this regime are visible in the FEDS intensity at t = 2 ps (Figure 3(b)). The ensuing hot-phonon population subsequently thermalizes with other lattice vibrations via phonon-phonon scattering, thereby driving the lattice towards thermal equilibrium (i.e., T qν = const.) within 50 ps, and leading to the thermalized FEDS intensity reported in Figure 3(c).

Conclusions
We

Sample preparation and thickness determination
The thin black phosphorus (BP) flake was obtained by standard mechanical exfoliation performed in air. The samples were then quickly imaged in the optical microscope and subsequently transferred to a load-lock chamber connected to our main experimental chamber in ultra-high vacuum. We estimate the total exposure to air to be less than one hour. We found that this method yielded diffraction patterns consistent with previous experimental works. 31 Given the multilayer nature of the samples (40 nm corresponds to roughly 80 layers), the observed scattering signals predominantly arise from the bulk as opposed to the oxidized surface layers of the flake. We note, however, the presence of forbidden reflections (h + l = 2n + 1) in the diffraction patterns. Such forbidden reflections were also observed in previous works, but their origin could not be attributed with certainty. 31 We postulate that they may be caused by stacking faults or structural deviations at the surface. These additional reflections do not alter the overall agreement between experiment and theory.
The flake thickness was estimated by transmission measurements in an optical microscope in combination with transfer matrix calculations and the optical constants of BP. 17

Computational details
First-principles calculations employed the primitive cell of bulk BP (point group D 2h and space group Cmce) that contains 4 atoms. 30 All calculations were performed using the PBE generalized gradient approximation 41 to density functional theory. We employed planewaves basis sets and Troullier-Martins norm-conserving pseudopotentials 42 as implemented in the Quantum ESPRESSO suite. 43 The planewaves kinetic energy cutoff was set to 90 Ry and the sampling of the Brillouin zone was performed using a uniform 12×10×10 k-point grid. We determine the interatomic force constants by means of density-functional perturbation theory calculations 44 using a 5×5×5 Brillouin-zone q-grid. The full set of phonon eigenmodes and eigenfrequencies was obtained by using standard Fourier interpolation of dynamical matrices on a 50×50×50 q-point grid. Such a dense grid guarantees a fine resolution of the calculated structure factor maps. The phonon band structure over a chosen high-symmetry path is shown in Figure S2.

Estimation of excited carrier density
The incident fluence on the BP flake is I inc = 98 J/m 2 . The pump was polarized along the armchair direction, determined by rotating a waveplate in the pump arm and maximizing the Debye-Waller effect at 50 ps (maximum absorption). Given the highly anisotropic optical absorption of BP, specifying the pump polarization is essential to employ the proper value of refractive index, which matters for the estimation of excited carried density. The complex refractive index of BP along the armchair direction at 800 nm is 3.19 + 0.29i, as estimated from a previous work. 1 We use transfer matrices to calculate the transmitted part, T, as well as the reflected part, R, of the incident fluence. This yields the absorbed fluence I abs = (1 − R − T ) · I inc . The carrier density per square centimeter is then: n = 10 −6 · 1 E ph · I abs d · c · 1 2 = (7.3 ± 0.9) · 10 13 electrons/cm 2 , where E ph is the energy of one pump photon, d is the flake thickness, which we estimate to 39 ± 5 nm, c = 10.46 · 10 −10 m is the unit cell length in the out-of-plane direction and the factor 1/2 accounts for the two layers per unit cell.
Diffuse scattering maps: raw data and influence of pump photon energy In supplementary Figure 1 we show a series of diffuse scattering maps at chosen delays and for different pump photon energies and/or data processing steps. Panels (a-c) and (j-l) correspond to the experimental and theoretical results shown in Figure 3 of the main text, respectively. They are reproduced here for convenience. Panels (d-f) display the diffuse scattering signals obtained when pumping the sample with 0.59 eV photons. The same initial excited electron density is assumed for both the 0.59 eV and 1.61 eV experiments, based on the Bragg peaks' reduction at 300 ps. While the diffuse scattering signatures at 2 ps seem slightly more intense for the 1.61 eV pump, we observe the same qualitative thermalization process in both cases. Since the picosecond dynamics does not seem to depend significantly on the initial carrier distribution after excitation, the assumption of a thermalized electron system for our purposes seems justified. In the main text and in panels (a-f), we show two-fold symmetrized data. This processed data is obtained using the open-source Python environment developed by De Cotret et al. 2 Prior to symmetrization, it was verified that the time-dependence of the elastic scattering for Friedel pairs showed the same dynamics within error margin. The raw data for the 1.61 eV pump are shown in panels (g-i). Comparing panels (a-c) and (g-i), we observe peak splitting effects at large scattering vectors for the symmetrized data. This artifact arises from field distortions of our magnetic lens, which contains aberrations at higher scattering vectors.

Supplementary Figures
Supplementary Figure 1: Diffuse scattering maps at 2 ps, 10 ps, and 50 ps for different pump photon energies and data processing steps.