Anisotropic Nonequilibrium Lattice Dynamics of Black Phosphorus

Black phosphorus has recently attracted significant attention for its highly anisotropic properties. A variety of ultrafast optical spectroscopies has been applied to probe the carrier response to photoexcitation, but the complementary lattice response has remained unaddressed. Here we employ femtosecond electron diffraction to explore how the structural anisotropy impacts the lattice dynamics after photoexcitation. We observe two time scales in the lattice response, which we attribute to electron–phonon and phonon–phonon thermalization. Pronounced differences between armchair and zigzag directions are observed, indicating a nonthermal state of the lattice lasting up to ∼60 ps. This nonthermal state is characterized by a modified anisotropy of the atomic vibrations compared to equilibrium. Our findings provide insights in both electron–phonon as well as phonon–phonon coupling and bear direct relevance for any application of black phosphorus in nonequilibrium conditions.

L ayered van der Waals (vdW) materials have attracted significant research interest in recent years due to their potential device applications. 1−4 The most prominent 2D material, graphene, exhibits high carrier mobility but lacks a band gap, which is required in many applications. In contrast, transition metal dichalcogenides possess a band gap in the visible range but a lower carrier mobility. With a thicknessdependent band gap extending from the infrared to the visible 5−7 and a high carrier mobility, 8−10 black phosphorus provides an important complementary building block for vdW heterostructure devices. A central aspect of black phosphorus is its in-plane anisotropic structure, shown in Figure 1a. The layers have two inequivalent high-symmetry directions, the socalled zigzag and armchair directions. This structural anisotropy is also reflected in many macroscopic material properties, such as optical absorption 9,11−14 and in-plane anisotropic thermal 15−18 and electrical 5,9,19,20 conductivities. These anisotropic properties offer additional tunability in device design.
Since any device operates in nonequilibrium conditions, a microscopic understanding of nonequilibrium states in vdW materials is of particular interest. For optoelectronic devices, knowledge of the evolution of the system after optical excitation is desired. Carrier dynamics in black phosphorus have been studied using a variety of time-resolved optical spectroscopies 20−26 as well as time-and angle-resolved photoemission spectroscopy 27,28 (trARPES). An important relaxation pathway for excited carriers is via coupling to the lattice. However, to date, no study has directly reported on the ultrafast lattice response of black phosphorus upon photoexcitation, which reflects the strength of electron−phonon as well as phonon−phonon interactions. In this work, we employ femtosecond electron diffraction 29 (FED) to directly probe the structural dynamics of photoexcited black phosphorus.
The measurement principle is sketched in Figure 1b. The sample is excited with an ultrashort laser pulse (pump), and the lattice response is probed using an ultrashort electron pulse (probe) with a kinetic energy of 70 keV. The electrons are diffracted by the sample and diffraction patterns are recorded in transmission for different time delays between the pump and probe pulses. To excite the sample, we use optical pulses with a wavelength of 770 nm (1.61 eV). The polarization of the pump pulse is set to the armchair direction of the crystal. All measurements are performed at a base temperature of 100 K.
Since diffraction patterns are measured in transmission, the samples need to be thin films. We prepared a thin film of black phosphorus by mechanical exfoliation from a bulk crystal (HQ Graphene) using water-soluble glue and transferred it on a standard copper transmission electron microscopy (TEM) grid using the floating technique. 30 Based on the transmission of the film and previously reported optical constants of black phosphorus, 14 we estimate the film thickness to be 39 ± 5 nm. The sample was transferred to vacuum directly after preparation to minimize degradation. Figure 2a shows a typical transmission diffraction pattern of black phosphorus. The orthorhombic crystal structure of black phosphorus (see Figure 1a) results in an anisotropic diffraction pattern. Bragg reflections along the zigzag and armchair directions, (h00) and (00l), are marked with blue and green boxes, respectively. Only reflections with even h and l are allowed due to the crystal symmetry (space group Cmce with phosphorus atoms at Wyckoff positions 8f). We observe additional "forbidden" reflections caused by stacking faults, multiple scattering, or structural deviations at the surfaces, which are not taken into account in the following analysis.
In this work, we focus on the Bragg reflections along the armchair and zigzag directions. Our primary observables are the intensities of the Bragg reflections, which decrease with increasing displacement of the atoms due to lattice vibrations (Debye−Waller effect). Since the probe electrons propagate in parallel to the van der Waals stacking direction through the crystal, our measurement is sensitive to the in-plane atomic vibrations.
To extract intensities from the diffraction patterns, we fit the observed peaks to the sum of a 2D pseudo-Voigt profile and a tilted background. The same Gaussian−Lorenzian mixing is assumed in all directions. The integrated intensity of the peak is then used to obtain the intensity change as a function of pump−probe delay. Figure 2b displays the resulting evolution of the Bragg reflection intensities after photoexcitation. We observe two time scales in the intensity decrease and pronounced differences between reflections along the armchair and zigzag directions. The amplitudes of the intensity decrease are larger for reflections with higher scattering vectors, as expected by the Debye−Waller theory.
To analyze the structural dynamics, we convert the Bragg reflection intensities into changes in atomic mean squared displacement (MSD). For anisotropic crystals, the temperature factor is 31 Here, h, k, and l are the Miller indices and a*, b*, and c* are the magnitudes of the reciprocal space lattice vectors, defined such that a·a* = 2π. For black phosphorus in the standard setting of Cmce, a, b, and c are the (real space) lattice vectors in zigzag (a = 3.31 Å), out-of-plane (b = 10.46 Å), and armchair (c = 4.37 Å) directions. 34 U 11 , U 22 , and U 33 correspond to the MSDs in the zigzag, out-of-plane, and armchair directions, respectively (see refs. 31 and 32 for more details about the U-matrix). Note that in black phosphorus, U 12 and U 13 are zero due to symmetry. 34,35 For Bragg reflections purely along the zigzag direction, that is (h00), the relative intensity change after laser excitation is given by Similarly, for reflections purely along the armchair direction, (00l), the relative intensity change after laser excitation reads Here, U ii 0 denotes the MSD before laser excitation and U ii (t) denotes the MSD as a function of pump−probe delay. I(t) denotes the intensity as a function of pump−probe delay and I 0 is the intensity before laser excitation.
The changes in MSD, shown in Figure 2c, are markedly different for the armchair and zigzag directions. From the different amplitudes we conclude that the interatomic potential is anisotropic. It is energetically less costly to displace atoms along the armchair direction compared to the zigzag direction. Hence, as the lattice temperature rises, the additional energy leads to a larger increase of the MSD in the armchair direction. These findings are in qualitative agreement with lattice dynamical calculations 36 and static X-ray diffraction measurements. 34 The MSD dynamics in the armchair and zigzag directions are fitted with a biexponential function, see solid lines in Figure  2c. The finite time resolution is taken into account by convolving the fit function with a Gaussian with a full width at half maximum (FWHM) of 150 fs. The fit results for the amplitudes A i and time constants τ i are summarized in Table 1.
The fast time constants τ 1 are very similar for the armchair and zigzag directions. We attribute this initial rise in MSD to energy transfer from the photoexcited electrons to the lattice. The observed ∼0.5 ps electron−lattice equilibration time constant is consistent with subpicosecond dynamics observed with time-resolved optical spectroscopy 23,24,26 and trARPES. 28 While the amplitudes of the initial MSD rise are the same for the two directions, the thermalized state at late delays exhibits Nano Letters pubs.acs.org/NanoLett Letter a higher MSD increase in the armchair direction compared to the zigzag direction. This indicates that electron−phonon equilibration leads to a nonthermal phonon distribution. Compared to the thermal case, this nonthermal phonon distribution is characterized by a higher ratio of the zigzag to armchair MSD. Hence, on average, phonons with a high displacement in zigzag direction couple more strongly to the electrons. The persistence of the nonthermal phonon distribution is illustrated in Figure 2d, in which the timedependent MSDs in the armchair and zigzag directions are normalized to the respective fit values at 100 ps. Within tens of picoseconds, the nonthermal phonon distribution relaxes into a thermal phonon distribution. We therefore attribute the slower time constants τ 2 of the MSD dynamics to these redistribution processes. The fact that the MSD rises further during phonon thermalization indicates that high-energy phonons decay into lower-energy phonons because one high-energy phonon decays into multiple low-energy phonons, and in addition, low-energy phonons produce a higher atomic displacement per phonon. 33 Since the conduction band minimum and the valence band maximum are located at the Z-point of the Brillouin zone, 5 thermalized carriers can mostly absorb and emit phonons with small momenta. Hence, we can conclude that lattice thermalization is achieved mostly by direct phonon−phonon coupling via anharmonicities. Finally, a thermal state is reached after ∼60 ps. This time scale of phonon thermalization is similar to time scales reported for the vdW materials WSe 2 37 and graphite. 38 In contrast to these materials, however, the nonthermal phonon population in black phosphorus manifests itself in an anisotropic evolution of the MSD due to the inplane anisotropy.
We performed similar experiments for different laser pump fluences, pump pulse polarizations, and for a different sample  The errors correspond to 68.3% confidence intervals of the fit.

Nano Letters
pubs.acs.org/NanoLett Letter base temperature. In all cases, we obtain values of ∼0.5 ps for the fast time constants τ 1 , which reflect the energy transfer from the electrons to the lattice. In contrast, the slow time constants τ 2 decrease with increasing sample base temperature. We attribute this to enhanced phonon−phonon scattering due to a larger phonon population at higher temperature. Correspondingly, we also observe an acceleration of the phonon thermalization with increasing excitation density. For a sample base temperature of 295 K and a laser pump fluence of ∼2 mJ/cm 2 , we obtain τ 2 values of (11.5 ± 0.7) ps and (8 ± 1) ps for the armchair and zigzag directions, respectively. These data, in addition to the low-temperature data presented here, are available on a data repository. 39 So far, we have only considered the MSD along armchair and zigzag directions. However, we can calculate the MSD along all in-plane directions based on the MSD along these two directions. 31,35 Figure 3 visualizes the evolution of the in-plane MSD, using the biexponential fit of the data shown in Figure  2c. To estimate the MSD before laser excitation, we make the assumption that both before and after laser excitation, the MSD is proportional to the temperature (high-temperature limit). We estimate that the temperature of the system rises by (270 ± 50) K after laser excitation, based on the calculated absorbed energy density and the heat capacity of black phosphorus. The total heat capacity is approximated by its main contribution, the lattice heat capacity, calculated from the vibrational DOS. 40 In equilibrium, the shape of the in-plane MSD is already anisotropic (blue curves of Figure 3) due to the anisotropic bond stiffness. After laser excitation, the MSD increases as shown in Figure 2c. The evolution of the MSD is displayed in Figure 3a. Note that not only the size but also the shape of the MSD evolves with time. This is visualized in Figure 3b by showing the MSD curves normalized to their area. The nonthermal phonon distribution leads to a transient reduction of the MSD anisotropy, with a larger MSD in the zigzag direction compared to equilibrium. On longer time scales, as the phonons thermalize, the MSD relaxes back to its equilibrium shape (dashed red curve).
In summary, our measurements have explored how structural anisotropy impacts lattice thermalization in photoexcited black phosphorus. We have shown that the lattice response is well captured by biexponential dynamics: a subpicosecond component, common to both zigzag and armchair Bragg reflections, is assigned to electron−phonon equilibration, while a ∼20 ps component is found to be highly anisotropic and indicative of a nonthermal phonon population persisting for ∼60 ps. Our analysis reveals that the nonthermal phonon population results in a transient shape change of the MSD, with a higher displacement in zigzag direction compared to equilibrium.
We expect this nonthermal state of the lattice to have effects on the macroscopic properties of black phosphorus, such as thermal and electrical conductivities. In particular, for any application in which hot carriers are excited or injected in black phosphorus, transient changes of material properties will influence device performance. Beyond black phosphorus, we expect these results to be relevant to other vdW materials as well. For example, in any layered material, the evolution of inplane and out-of-plane MSD could be different, with implications for energy flow across heterostructure interfaces.  (a) Evolution of the MSD from before excitation (blue) to the thermalized state at 100 ps (red). The shape of the in-plane MSD before laser excitation is already anisotropic due to the in-plane structural anisotropy. Note that the MSD increases due to lattice heating but also the shape of the MSD changes transiently. (b) MSD curves from panel a, normalized to their area. This highlights the transient shape change of the MSD due to a nonthermal phonon distribution.