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.

Layered 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 thickness-dependent 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 a. The layers have two inequivalent high-symmetry directions, the so-called 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.

Figure 1

Anisotropic structure of black phosphorus and schematic illustration of the experiment. (a) Side and top views of the crystal structure of black phosphorus, showing anisotropy not only between the in-plane and out-of-plane directions but also between the in-plane directions (armchair and zigzag). (b) Schematic representation of a time-resolved electron diffraction setup (see text for details).

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 b. 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 a shows a typical transmission diffraction pattern of black phosphorus. The orthorhombic crystal structure of black phosphorus (see Figure a) 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.

Figure 2

Overview of anisotropic lattice dynamics in photoexcited black phosphorus. (a) Transmission diffraction pattern of thin-film black phosphorus. We focus on the high-symmetry Bragg reflections along the armchair and zigzag directions, indicated by colored boxes. (b) Relative changes in Bragg reflection intensities as functions of pump–probe delay. Here, we average over the Friedel pairs, e.g., (002) and (002̅), since they show the same dynamics. The measurement was conducted with an incident fluence of (9.8 ± 1.4) mJ/cm2. Based on the optical constants of black phosphorus[14] and the film thickness, we estimate the absorbed energy density to be (380 ± 70) J/cm3. The data presented are the average of several delay scans, and the error estimates represent the standard error of the mean. (c) Changes of atomic mean squared displacement (MSD) in the armchair (green circles) and zigzag (blue triangles) directions as functions of pump–probe delay. The anisotropy of the lattice is reflected in an anisotropic MSD change in the two directions. The higher MSD change in the armchair direction indicates that bonds are softer in this direction. The MSD values presented are the weighted average of MSDs calculated from each Friedel pair, and the error bars are calculated using error propagation. The inset is a close-up of the data at early time delays. (d) MSD changes normalized to the fit values at 100 ps. A two-step time scale as well as a transient nonequilibrium between the zigzag and armchair directions is observed.

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 b 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−33]

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]U11, U22, and U33 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, U12 and U13 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, U0 denotes the MSD before laser excitation and U(t) denotes the MSD as a function of pump–probe delay. I(t) denotes the intensity as a function of pump–probe delay and I0 is the intensity before laser excitation.

The changes in MSD, shown in Figure c, 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 c. 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 and time constants τ are summarized in Table .

Table 1

Fit Results of the MSD in the Armchair and Zigzag Directions with a Biexponential Function Convolved with a Gaussiana

	armchair	zigzag
A1 [10–3 Å2]	0.66 ± 0.03	0.68 ± 0.02
τ1 [ps]	0.48 ± 0.05	0.58 ± 0.04
A2 [10–3 Å2]	1.84 ± 0.03	0.85 ± 0.03
τ2 [ps]	22 ± 1	20 ± 2

The errors correspond to 68.3% confidence intervals of the fit.

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 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 d, in which the time-dependent 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 WSe2[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 in-plane anisotropy.

We performed similar experiments for different laser pump fluences, pump pulse polarizations, and for a different sample 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/cm2, 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 visualizes the evolution of the in-plane MSD, using the biexponential fit of the data shown in Figure c. 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]

Figure 3

Photoexcitation transiently modifies the anisotropy of the atomic mean squared displacement (MSD). (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.

In equilibrium, the shape of the in-plane MSD is already anisotropic (blue curves of Figure ) due to the anisotropic bond stiffness. After laser excitation, the MSD increases as shown in Figure c. The evolution of the MSD is displayed in Figure a. Note that not only the size but also the shape of the MSD evolves with time. This is visualized in Figure b 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 in-plane and out-of-plane MSD could be different, with implications for energy flow across heterostructure interfaces.