Ultrafast anisotropic exciton transport in phosphorene

Abstract

Phosphorene is a two-dimensional (2D) material exhibiting strong in-plane structural anisotropy. In this work, we investigate the influence of structural anisotropy on the optics, dynamics, and transport of excitons in phosphorene by combining microscopic many-body theory with first principles calculations. Our framework offers a complete and material specific description of the excitonic properties of phosphorene, including exciton states and exciton–phonon interactions, which allow us to quantitatively evaluate the optical absorption spectra, exciton relaxation, and exciton transport, revealing direction-dependent characteristics. Interestingly, we identify the critical role of long-range exchange interactions, which significantly enhance the anisotropy of exciton diffusion, particularly at low temperatures. Our work provides fundamental insights into exciton dynamics in an intrinsically anisotropic 2D material, offering guiding principles for the design of next-generation optoelectronic devices.

---

1. Introduction

Atomically thin two-dimensional (2D) materials have garnered significant attention over the past two decades due to their unique properties arising from dimensional confinement.^1–4 Compared to their bulk counterparts, 2D materials exhibit distinct electronic, optical, mechanical, and thermal characteristics, making them highly suitable for next-generation ultrathin devices such as high-speed transistors,^5,6 flexible displays,⁷ and sensitive photodetectors,⁸ while also holding promise for revolutionary energy technologies including high-performance batteries,^9,10 supercapacitors,^9,11 and efficient catalysts.^12,13

Coulomb-bound electron–hole pairs, known as excitons, play a fundamental role in the optical properties of semiconducting 2D systems.¹⁴ The reduced dimensionality and reduced dielectric screening in 2D materials compared to 3D materials enhance electron–hole interactions,¹⁵ resulting in large exciton binding energies and strong oscillator strengths.^16–19 These features allow excitons to remain stable at room temperature^20,21 and they can be tuned via external fields,^22,23 strain,^24–26 or environmental factors.^27,28 This rich landscape of excitonic phenomena in 2D materials is not only of scientific interest but also central to the development of optoelectronic devices, including light-emitting diodes (LEDs),²⁹ photodetectors,⁸ solar cells,³⁰ optical modulators,^31,32 and quantum light sources.³³

Anisotropy in materials offers a powerful means to control optical properties and energy transport,³⁴ opening new pathways for device design and functionality. Phosphorene, a two-dimensional material with inherently large in-plane anisotropy, has emerged as a promising platform for exploring such effects. The anisotropic optical response of phosphorene has been widely observed,^35–38 evident in linearly polarised optical absorption and the emergence of hyperbolic exciton-polaritons. However, despite these intriguing findings, the dynamics and transport behavior of excitons in phosphorene remain largely unexplored.

In this work, we propose that the intrinsic anisotropy of phosphorene leads to highly directional exciton energy transport. We use microscopic many-body particle theory³⁹ parametrised with first principles calculations to investigate the exciton properties of phosphorene. We compute: (i) the exciton states and corresponding optical absorption spectra to elucidate their directional dependence; (ii) the exciton relaxation dynamics following photoexcitation; and (iii) the anisotropic exciton diffusion constants along distinct crystallographic directions, which are found to be dramatically affected by long-range exchange interactions. Our findings demonstrate that the inherent exciton anisotropy in phosphorene is robust against phonon scattering, opening the door for intrinsic directional control of energy in a device. This will pave the way for new applications including exciton highways for light harvesting applications, linearly polarised light-emitting diodes, and future excitonic circuitry.

2. Methods

2.1. First principles calculations

We perform structural optimization, and calculate the electronic band structure and the macroscopic dielectric constant using density functional theory as implemented in the Vienna Ab initio Simulation Package (VASP)⁴⁰ with the projector augmented wave method.^41,42 We use a 16 × 12 × 1 Monkhorst–Pack k-grid to sample the Brillouin zone⁴³ and set the energy cut-off to 500 eV. To achieve an accurate treatment of the exchange–correlation energy, we adopt the hybrid HSE06 functional for structural optimization and electronic band structure calculations.⁴⁴ We maintain a vacuum spacing greater than 15 Å to prevent interlayer interactions, and optimise the structure to achieve residual forces below 0.01 eV Å⁻¹.

Using wannier90,⁴⁵ we construct a tight binding model parametrised by the HSE06 calculations. We use this tight binding model to finely sample the electronic bands around select k-points, we then fit the resulting two-dimensional band energy surfaces to an analytical functional form using linear regression, which we then use to construct the associated Hessian matrix, and finally evaluate the electron and hole effective masses.

For the macroscopic dielectric constant, we adopt the Perdew–Burke–Ernzerhof (PBE) generalised gradient approximation (GGA) to the exchange–correlation functional.⁴⁶ We use the finite difference method implemented in VASP to compute the dielectric tensor. The effective monolayer thickness and correction of the 2D dielectric constant follows the methodology proposed in previous works.⁴⁷

We use Quantum ESPRESSO⁴⁸ to calculate the phonon dispersion and electron–phonon matrix elements of phosphorene. We adopt norm-conserving pseudopotentials in conjunction with the PBE exchange–correlation functional.⁴⁶ For phonons, we construct the dynamical matrix from the calculated force constants to get phonon energies at arbitrary q-points. We also use QUANTUM ESPRESSO to calculate the interband dipole moment, defined as the probability of electron transitions between different energy bands, for light absorption calculations.

For electron–phonon interactions, the matrix elements have the following form:

where ω^j_q is the phonon frequency associated with wave vector q and branch j, α is the electronic band index, and ρ is the mass density of the material. In phosphorene, the optical phonon modes are relatively flat and excitons are highly localised in momentum space. As a result, the exciton–phonon interaction associated with optical modes can be effectively described by transitions involving a narrow range of phonon energies and momenta, and we employ a constant deformation potential with respect to the phonon wave vector q. For acoustic modes, we find that the deformation potential scales linearly with phonon wave vector for small q, in agreement with previous works,⁴⁹ that is, D^k,j,α_q = D^k,j,α_AC |q|. Overall, we interpolate the electron–phonon coupling coefficients using for acoustic linear modes and keep the value constant for optical modes. The acoustic quadratic mode characteristic of 2D materials corresponds to out-of-plane atomic vibrations which have vanishing electron–phonon coupling in centrosymmetric materials. Due to the computational cost, we account for the directional dependence by evaluating the first principles values along several high-symmetry directions, and then use linear interpolation for general directions. While the interpolation scheme may overlook subtle variations in the electron–phonon coupling, it still effectively captures the dominant anisotropic trends driven by the crystal structure, capturing the behavior along high symmetry lines.

2.2 Wannier equation

To model the exciton behavior, we use the Wannier equation to describe the electron–hole interaction:^50,51where m_r is the reduced mass arising from the electron m_e and hole m_h effective masses, and V_q is the screened Coulomb electron–hole interaction with momentum difference q between electron and hole. The excitonic binding enegy E^b_μ and wavefunction φ^μ(k) for state μ can be found by diagonalising the above equation, which gives a hydrogen-like series of excitonic energy levels. The parameter k of the exciton wavefunction is the relative electron–hole momentum, k = αk_h + βk_e, where α = m_e/(m_e + m_h) and β = m_h/(m_e + m_h).

The Coulomb interaction in the Wannier equation can be modelled using the two-dimensional Keldysh potential:^52,53

where ε₀ is the vacuum dielectric constant, ε_background is the background screening, d is the thickness of the material, and A is the surface area per unit cell. Throughout this study, we consider phosphorene on a SiO₂ substrate with a background ε_background of 3.9,⁵⁴ to simulate realistic experimental conditions.

We treat ε^∥_material as the average of the dielectric constant in the two in-plane directions, and we compare the excitonic calculation results with those using the interpolation method for the Keldysh potential.⁵⁵ The resulting excitonic wavefunction distributions are nearly identical. Our calculations show that the linewidth difference between the average and elliptical dielectric screening models is less than 0.001 eV, indicating that phonon-driven process, such as linewidths and diffusion, are not significantly influenced by the anisotropy of the dielectric constants. Following eqn (2), electron–hole interactions and the |q|-dependence of the Keldysh potential, eqn (3), it is clear that the excitonic binding energies and wavefunctions are driven mostly by small |q|. This means the Keldysh potential takes the form V ∝ 1/ε_background|q|, such that substrate screening, rather than material screening (anisotropic or otherwise) dominates the interaction.

The exciton–phonon scattering matrix elements:³⁹

describe scattering between exciton state η with centre-of-mass momentum Q to the exciton state η′ with centre-of-mass momentum Q^Prime. In this expression, G^ηη_Q−Q′ denotes the exciton–phonon scattering coefficient, which has the form , where c and v denotes the conduction and valence band, respectively; A^ν,±_Q−Q′ = (δ_±1,1 + n^ν_Q−Q′)describes emission (+) and absorption (−) of a phonon; and the δ function imposes the exciton–phonon scattering channel, which in numerical calculations is replaced by a Lorentzian. The use of the Lorentzian can be justified by noting that truncating the semiconductor Bloch equation at higher order (including two-phonon processes) softens the usual delta function, giving rise to a self-consistent exciton–phonon scattering rate.⁵⁶ This expression for the exciton–phonon scattering rate can be derived using the semiconductor Bloch equations.⁵⁰

The dephasing describes the total scattering rate of an exciton in a specific state:³⁹

arising from the summation of the exciton–phonon scatterings with all other excitons of momenta Q′ in stateη′.

3. Results and discussion

3.1 Crystal structure and physical parameters

Fig. 1(a) shows the crystal structure of phosphorene, a single layer of black phosphorous. Unlike the completely planar structure of graphene, phosphorene has a buckling structure. The primitive cell of phosphorene adopts an anisotropic rectangular shape and contains four phosphorous atoms. The calculated lattice constants of phosphorene are a₁ = 3.277 Å and a₂ = 4.576 Å, and the calculated effective monolayer thickness is 5.3 Å. By comparison, the experimental lattice constants of bulk black phosphorous are a₁ = 3.31 Å, a₂ = 4.38 Å, and a₃ = 10.5 Å, and the effective thickness of a single layer is 5.25 Å.⁵⁷

Fig. 1 (a) Top view and side view of phosphorene. The directions along a₁ and a₂ are referred to as zigzag and armchair directions, respectively. (b) HSE06 band structure of phosphorene. The valley along Γ–X is labeled V throughout this work. The Brillouin zone is shown in the inset. (c) Phonon spectrum of phosphorene.

In this work, we follow the convention to name a₁ and a₂ as armchair (AC) and zigzag (ZZ) directions, respectively.

The calculated macroscopic dielectric constants are ε_∥(ZZ) = 12.17, ε_∥(AC) = 16.13, and ε_⊥ = 6.73.

The calculated band structure of phosphorene is shown in Fig. 1(b). It exhibits strong anisotropy around the valence band maximum (VBM), with the effective masses along the Γ–X (ZZ) and Γ–Y (AC) directions being 19.670m₀ and 0.849m₀, respectively, where m₀ is the free electron mass. For the conduction band minimum (CBM), the anisotropy is less pronounced: at the Γ valley the effective masses along the ZZ and AC directions are 6.840m₀ and 0.963m₀, respectively; and around the secondary CBM valley (denoted as V), the effective masses are 1.544m₀ and 1.249m₀ for the ZZ and AC directions, respectively.

ReSe₂ is a notable member in TMDs materials showing anisotropy,⁵⁸ where the effective masses can be 2.81 times difference for the hole, and 6.16 times difference for electrons. Nevertheless, as mentioned above, the effective masses can be 23.2 times difference for the hole, and 6.8 times difference for electrons for phosphorene, indicating that phosphorene has much larger anisotropy than ReSe₂.

Fig. 1(c) shows the phonon dispersion of phosphorene. The lowest energy acoustic mode is the flexural ZA mode,⁵⁹ exhibiting the quadratic dispersion characteristic of 2D materials. Tables 1 and 2 list the calculated electron–phonon coupling coefficients. The modes exhibiting the largest electron–phonon coupling coefficients are selected for exciton–phonon coupling calculations. In general, low-energy acoustic modes drive intraband scattering, while the higher-energy optical modescan drive interband transitions or large Q transitions. As a result, optical phonon modes dominate the relaxation dynamics, while acoustic phonon modes dictate the linewidth and diffusion, especially at low temperature.

Table 1 Electron–phonon coupling coefficients g_c and the corresponding phonon energies E_ph for selected phonon modes coupling with CBM states. The electron–phonon coupling coefficients show anisotropy around q = 0, so we use the format for this case. For the optical and intervalley acoustic modes, the units are given by meV, and for intravalley acoustic modes, the units are given by meV nm^−1/2. The units of phonon energy are meV

Electron–phonon coupling coefficients/phonon energy (CBM)
Phonon mode	Γ → Γ		V → V		Γ → V
	gc	Eph	gc	Eph	gc	Eph
Acoustic 1	∼0		∼0		10.7	10.5
Acoustic 2	[∼0, 11.3, ∼0]	[5.1, 15.7, ∼0]	[∼0, 9.6, 6.6]	[5.1, 15.7, ∼0]	∼0
Acoustic 3	[44.1, 31.3, 31.2]	[10.4, 8.8, 5.1]	[39.4, 30.0, 18.1]	[10.4, 8.8, 5.1]	∼0
Optical 3	∼0		∼0		92.8	32.4
Optical 4	116.0	42.5	∼0		∼0
Optical 8	105.0	55.8	131.0	55.8	∼0

---

Table 2 Electron–phonon coupling coefficients g_v and the corresponding phonon energies E_ph for selected phonon modes coupling with VBM states. The format of data and the units are the same as in Table 1

Electron–phonon coupling coefficients/phonon energy (VBM)
Phonon mode	Γ → Γ
	gv	Eph
Acoustic 3	[14.9, 7.9, 3.1]	[10.4, 8.8, 5.1]
Optical 4	6.4	42.5
Optical 8	73.2	55.8

---

For thermal expansion, previous work⁶⁰ has shown lattice expansions in phosphorene of 0.1% and 0.08% for the AC and ZZ directions from 0 K to 300 K, respectively. These values indicate that even though phosphorene does show anisotropic behavior in its thermal lattice expansion, the value is too small to have a significant impact on the band dispersion and hence on the excitonic properties. Therefore, the role of thermal expansion is not considered throughout this work.

3.2. Exciton landscape and optical absorption

Fig. 2(a–d) shows the calculated spatial distribution |φ²_μ| of the direct ΓΓ excitonic wavefuntions obtained by solving the Wannier equation. The anisotropy of phosphorene leads to the 1s state exhibiting an elliptical shape in momentum space, which is consistent with the band structure shown in Fig. 1(b). Specifically, the electronic band structure has a broader band dispersion along Γ–X (the ZZ direction), leading to a more extended exciton wavefunction in this direction. Similarly, the spatial distribution of the 2p_x state is broader than that of the 2p_y state. Additionally, the anisotropy lifts the degeneracy of the phosphorene 2p_x and 2p_y excitons, a scenario that is markedly distinct to the corresponding excitonic landscape in isotropic systems such as transition metal dichalchogenides. Notably, our calculated exciton binding energy of E^b_1s = 0.29 eV for the 1s state is in a good agreement with the experimental value of 0.3 eV measured for phosphorene one SiO₂.⁶¹ Unlike the exciton energy, binding energies calculated by the Wannier equation are independent from the band gap value. Their relationship can be described as E_exciton = E_gap − E_binding.

Fig. 2 Spatial distributions of the (a) 1s, (b) 2p_x, (c) 2s, and (d) 2p_y states of bright ΓΓ excitons, with corresponding binding energies written in yellow. (e) Calculated light absorptions for the first three bright exciton states. The total absorptions along the AC and ZZ directions are depicted with solid lines, and contributions from individual states with dashed lines. Note that the absorption vanishes when the light polarization is along the ZZ direction.

The characteristic linear polarization of excitons in phosphorene is determined by the interband dipole matrix elements connecting the VBM and CBM. The numerical evaluation of these matrix elements, combined with the associated optical selection rules, show that the transition dipole is oriented along the AC direction, with a magnitude of 6.862 nm⁻¹.

Light absorption can be described by the Elliot formula,⁶² which is obtained from the semiconductor luminescence equations³⁹ and derived from the Heisenberg equation of motion:⁵⁰

where on the left-hand-side, I denotes the absorption intensity, and σ represents the polarization of the absorbed light. On the right-hand-side, in the numerator, M̃^σ_η is the optical matrix element, given by the expression , where d_η is the interband dipole of the valley for the given stateη, and ê_σ is the unit vector of incident light polarization. In the denominator, the first term E^b_η is the exciton binding energy for the exciton stateη; the second term ħω corresponds to the energy of the incident light; and the final term accounts for two broadening mechanisms: the optical broadening γ^σ_η arising from radiative decay is defined by:where n is the refractive index and M^cv_σ is the interband dipole between conduction and valence bands. Th phonon-induced broadening Γ⁰_η represents the dephasing of state η at zero center-of-mass momentum.

Indirect Γ–V excitons are optically dark due to large momentum transfer and hence do not couple to light.

The calculated interband dipole indicates that the dipole is aligned along the AC direction, and light absorption will be greatly suppressed when the polarization of the incident light is along the ZZ direction which is perpendicular to the AC direction. Fig. 2(e) shows the absorption spectrum of phosphorene, with a nearly vanishing absorption when light polarization is along the ZZ direction. There are two visible peaks when the light polarization is parallel to the interband diople (the AC direction), which stem from the 1s and 2s exciton states. In contrast, since the parity of p_x is odd, the absorption of state 2p_x is forbidden by symmetry. The sum in M̃^σ_η is zero due to the odd parity. These anisotropic results are in good agreement with previous works.⁶³

3.3. Phonon-mediated exciton dynamics

Exciton dynamics, the time evolution of excitons scattering between states mediated by phonons, is a key factor in understanding exciton linewidths, relaxation, and transport.

3.3.1. Exciton dispersion and linewidth. In studying exciton dynamics, an accurate description of the exciton energy is needed in order to reliably track the evolution of exciton states. The excitonic dispersion can be described as a parabolic dispersion, stemming from the electronic and hole valleys, and a Q-dependent long-range exchange interaction K_Q:where M = m_e + m_h is the total mass of the exciton. The long-range exchange interaction resembles a dipole–dipole interaction and has the form:⁶⁴where A is the area per unit cell, and d^valley_k denotes the interband dipole matrix element within a given valley. We include the exchange interaction only for the ΓΓ excitons, since the large momentum mismatch in ΓV excitons suppresses the long-range interaction, and we neglect the short-range exchange interaction as it is typically weak in two-dimensional materials.⁶⁵ The short-range exchange interaction has a negligible effect on the exciton dispersion in 2D materials.^65,66 While the direct and long-range exchange interactions depend on small momentum transfers q ≪ G, such that the Coulomb interaction decays slowly as 1/q, the short-range depends on large momentum transfers q > G and decays very quickly 1/q². Therefore the short range interaction is typically small in these systems. The long-range exchange interaction can be expressed as a dipole–dipole interaction, and therefore it depends on the transition dipole which is highly anisotropic. In phosphorene, the optical transition dipole is largely aligned along the AC direction. Since exctions can be regarded as dipoles within the material, their energies are strongly affected by the long-range interaction, leading to the change of energy dispersion. In other words, the long-range interaction couples to the exciton's center-of-mass momentum. The importance of long-range exchange interactions has been observed in many other works.^67–70 For example, first-principles calculations of phosphorene have demonstrated the importance of these long-range exchange interactions in the resulting exciton band dispersion.⁶⁶ The Keldysh potential takes the form V(Q) ≈ c/(aQ² + bQ), implying that for small Q, the linear term bQ dominates and K_Q scales linearly with Q. For large Q, K_Q tends to a constant, and the energy dispersion adopts to parabolic form driven by the shape of the electron and hole valleys.

Fig. 3(a) shows the energy dispersion of the 1s state of the ΓΓ excitons. The intrinsic anisotropy of the effective masses leads to more dispersive excitonic energy bands along the AC direction compared to the ZZ direction. Moreover, since the dipole moment d is oriented along the AC direction, the scalar product Q·d vanishes in the ZZ direction, further enhancing the anisotropy in the energy dispersion. As discussed above, the pink dashed line of Fig. 3(a) highlights how the original parabolic energy band is distorted by the exchange-induced linear term for small Q, but exhibits the expected parabolic profile at larger Q. The exchange contribution for small Q significantly affects exciton diffusion, as will be discussed later.

Fig. 3 (a) Total exciton energy of the 1s state. The figure is divided into two parts: the left shows the ZZ direction, and the right shows the AC direction. The band distortion indicated by the pink dashed line shows the effect of the exchange energy. (b) Linewidth Γ_Q=0 with and without exchange interaction. The red and black curves correspond to results without and with exchange interaction, respectively. (c) Schematic illustration of the exciton relaxation cascade. (d) Time evolution of the incoherent exciton population. The light pulse is shown with orange shadow area, and the total population is shown with solid line. Individual population of each states is demonstrated with different colours, as indicated in the legend.

We show the linewidth Γ^1s₀ in Fig. 3(b), which is significantly different when the long-range exchange interactions are included or excluded. Including the exchange interaction significantly reduces the linewidth because the linear dispersion at low Q makes scattering from Q = 0 less likely, owing to the reduced excitonic density of states and more restrictive energy and momentum conservation.

3.3.2 Exciton relaxation. Fig. 3(c) provides a schematic illustration of the exciton relaxation cascade. The final thermalised exciton population is independent of the initial exciton state that is excited. In our calculations, we choose to excite coherent excitons in 2s states by a resonant light pulse in order to illustrate how excitons relax from a higher to lower energy states. While 2p states are not optically active, we do include 2p states in the relaxation cascade because the mechanism driving the relaxation cascade is exciton–phonon coupling, so the excitonic population redistribution does not rely on the optical selection rules. During the relaxation process, coherent excitons in 2s states are scattered by phonons to other incoherent exciton states.³⁹ This polarization-to-population transfer is modeled by the following two microscopic differential equations:

The first equation describes the process of creating a coherent excitonic polarization, where is the excitonic polarization driven by a Lorentzian-like vector field A_σ(t) representing the laser pulse. We assume that the light polarization is parallel to the interband dipole. The second equation describes the exciton dynamics of population N^η_Q across exciton state η at momentum Q. The first term on the right-hand side determines the population transfer among exciton states, which will lead to the thermal equilibration among exciton states at long times. The final term accounts for the incoherent population generated from coherent excitons via exciton–phonon scattering.

As shown in Fig. 3(d), a 60 fs Lorentzian pulse excites the 2s state of the ΓΓ exciton, which also increases the population of all other exciton states via phonon scattering. We do not consider the recombination of the electron–hole pair forming the exciton, so the total exciton population saturates after the pulse. At the beginning, the populations of 2p_x, 2s and 1s states of ΓΓ excitons grow. In this regime, the scattering from the 2s-ΓΓ state into the 2p_x-ΓΓ exciton is much faster than that into the 1s-ΓΓ exciton because the energy difference between the final 2p_x-ΓΓ state and the initial 2s-ΓΓ state is approximately 40 meV at Q = 0, which closely matches the phonon energy of the third optical phonon mode, causing stronger intravalley scattering matrix elements for this transition compared to the 2s to 1s transition. From about 40 fs, the subsequent dynamics between all states leads to the populations of the 2p_x and 2s states starting to fall, while that of the 1s state continues to rise, and its population surpasses those of the 2p_x and 2s-ΓΓ excitons around 80 and 170 fs, respectively. As the system relaxes towards thermal equilibrium, the long-time population becomes largest for the lowest energy 1s-ΓΓ state, while the populations of the 2p_x ΓΓ -exciton and the 1s ΓV-exciton converge to similar values beyond 300 fs because the energies of these states are very similar. The 2s ΓΓ-exciton has the lowest population since its exciton energy is higher than that of other states. Finally, we omit the 2p_x and 2s states of the ΓV exciton in the figure due to their negligible contributions. We refer the reader to Fig. 2 and Tables 1, 2 for the excitonic binding energies, electron–phonon coupling coefficients, and phonon mode energies of the system that drive the above relaxation processes.

3.3.3 Exciton transport. Exciton transport can be described by tracking the spatiotemporal evolution for an initial exciton distribution. In an anisotropic system, this is well-described by Fick's second law of diffusion,^34,71 and is given by ∂_tN(r) = (D_x∂_x² + D_y∂_y²)N(r). In the formula, t denotes time, N is the exciton population, and D_i represents component i of the diffusion constant:where v^η,i_Q is the group velocity associated with the exciton in state η and centre-of-mass momentum Q along the ith direction, Γ^η_Q is the dephasing, E_Q,η is the exciton energy, and is the partition function.

Here, we assume the system is under thermal equilibrium, representing steady state diffusion. While time-dependent exciton transport is beyond the scope of this work, we predict that the impact of intermediate scattering channels, such as the ΓV and 2p excitons, will not significantly increase the diffusion. In fact, the lack of long-range exchange in these optically dark excitons will reduce their group velocity making them slower than the 1s ΓΓ exciton states.

In Fig. 4(a and b), we compare the diffusion constants with and without the inclusion of the exchange interaction. In both cases, D_AC (D_x) is much larger than D_ZZ (D_y), which can be rationalized through the associated band dispersions. Fig. 4(b) shows that even without long-exchange exchange interactions there is a large diffusion anisotropy in phosphorene. Specifically, the diffusion coefficient is proportional to the square of the group velocity, which in turn can be expressed as , where E_Q,η has the form in eqn (8). This implies that the diffusion constant is inversely proportional to the squared value of the effective mass, explaining why the diffusion coefficient D_AC is much larger than D_ZZ. The directional dependence of diffusion coefficients is significantly enhanced by the exchange interaction. Since the exchange energy is proportional to the scalar product Q·d, the energy influence is maximized along the dipole direction (the AC direction).

Fig. 4 Calculated diffusion coefficients of phosphorene. (a) With exchange interaction (b) Without exchange interaction. The red and blue curves correspond to AC and ZZ directions, respectively. The 2D dephasing strength at (c) 50 K and (d) 300 K. Thermal expectation values of at (e) 50 K and (f) 300 K. Arrows indicate the direction of .

Another important feature of the diffusion coefficients is their temperature dependence. As shown in Fig. 4(a and b), the diffusion coefficients (with or without exchange) decrease in the high temperature regime. This is driven by an increase in dephasing with rising temperature, leading to a reduction in the diffusion constants at high temperature. Fig. 4(c) and (d) show the dephasing landscapes at temperatures of 50 K and 300 K, respectively, illustrating that the values at 300 K are significantly greater than those at 50 K. In the opposite low-temperature limit, the value of the diffusion constant in the AC direction soars. This is driven by a combination of the fact that the exciton population accumulates around Q ≈ 0 at low temperature, and the steep gradient of the energy dispersion for small Q due to exchange along the AC direction.

Fig. 4(e) and (f) show the landscapes of the thermal expectation value of at 50 K and 300 K, respectively. Exciton velocity vectors in these two figures evidently align more with the AC direction, and the magnitude is bigger with small Q, especially at low temperature. This anisotropic enhancement in mobility has direct implications for temperature-tunable excitonic device design.

Finally, we note that the calculated values of the diffusion coefficient for the AC and ZZ directions at 300 K are 7.40 and 0.48 cm² s⁻¹, with their average value similar to the experimental value for a 2 nm amorphous phosphorene sample at 5.0 cm² s⁻¹.⁷²

4. Conclusions

We have investigated the excitonic behavior of phosphorene using a fully microscopic framework, with all essential parameters calculated from first principles. Our results illustrate the large impact that the intrinsic structural anisotropy of phosphorene has on its excitonic properties. The anisotropy in optical absorption originates from the directionality of the interband dipole moment, which aligns along the armchair (AC) direction, in good agreement with experimental observations. For exciton transport, our calculations reveal that both exciton effective masses and long-range exchange interactions significantly contribute to highly different diffusion constants in the two in-plane directions, with exchange making a significant contribution in enhancing diffusion in the AC direction. Our findings unveil the intrinsic anisotropy in the excitonic behavior of phosphorene.

Conflicts of interest

There are no conflicts to declare.

Data availability

First principles calculations were performed using Quantum Espresso (https://www.quantum-espresso.org/) and VASP (https://www.vasp.at/). Crystal structure data for phosphorene was obtained from https://www.crystallography.net/cod/. All other calculations can be reproduced using the equations in the main text.

Acknowledgements

K.-W. C. and B. M. acknowledge support from a UKRI Future Leaders Fellowship [MR/V023926/1] and J. J. P. T. and B. M. acknowledge support from a EPSRC Programme Grant [EP/W017091/1]. B. M. also acknowledges support from the Gianna Angelopoulos Programme for Science, Technology, and Innovation. The computational resources were provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service and funded by EPSRC [EP/P020259/1], and by the UK National Supercomputing Service, ARCHER2, for which access was obtained via the UKCP consortium and funded by EPSRC Grant [EP/P022561/1].

References

1. A. K. Geim and K. S. Novoselov, The rise of graphene, Nat. Mater., 2007, 6, 183–191.
2. S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev and A. Kis, 2D transition metal dichalcogenides, Nat. Rev. Mater., 2017, 2, 1–15.
3. D. Akinwande, C. Huyghebaert, C.-H. Wang, M. I. Serna, S. Goossens, L.-J. Li, H.-S. P. Wong and F. H. Koppens, Graphene and two-dimensional materials for silicon technology, Nature, 2019, 573, 507–518.
4. S. Das, A. Sebastian, E. Pop, C. J. McClellan, A. D. Franklin, T. Grasser, T. Knobloch, Y. Illarionov, A. V. Penumatcha and J. Appenzeller, et al., Transistors based on two-dimensional materials for future integrated circuits, Nat. Electron., 2021, 4, 786–799.
5. D. Lembke and A. Kis, Breakdown of high-performance monolayer MoS2 transistors, ACS Nano, 2012, 6, 10070–10075.
6. Y. Liu, X. Duan, H.-J. Shin, S. Park, Y. Huang and X. Duan, Promises and prospects of two-dimensional transistors, Nature, 2021, 591, 43–53.
7. A. K. Katiyar, A. T. Hoang, D. Xu, J. Hong, B. J. Kim, S. Ji and J.-H. Ahn, 2D materials in flexible electronics: recent advances and future prospectives, Chem. Rev., 2023, 124, 318–419.
8. M. Long, P. Wang, H. Fang and W. Hu, Progress, challenges, and opportunities for 2D material based photodetectors, Adv. Funct. Mater., 2019, 29, 1803807.
9. M. F. El-Kady, Y. Shao and R. B. Kaner, Graphene for batteries, supercapacitors and beyond, Nat. Rev. Mater., 2016, 1, 1–14.
10. Y. He, X. Shen and Y. Zhang, Layered 2D Materials in Batteries, ACS Appl. Nano Mater., 2024, 7, 27907–27939.
11. M. Tomy, A. Ambika Rajappan, V. M. Vimuna and X. Thankappan Suryabai, Emergence of novel 2D materials for high-performance supercapacitor electrode applications: a brief review, Energy Fuels, 2021, 35, 19881–19900.
12. T. A. Shifa, F. Wang, Y. Liu and J. He, Heterostructures based on 2D materials: a versatile platform for efficient catalysis, Adv. Mater., 2019, 31, 1804828.
13. D. Deng, K. Novoselov, Q. Fu, N. Zheng, Z. Tian and X. Bao, Catalysis with two-dimensional materials and their heterostructures, Nat. Nanotechnol., 2016, 11, 218–230.
14. R. Perea-Causin, D. Erkensten, J. M. Fitzgerald, J. J. Thompson, R. Rosati, S. Brem and E. Malic, Exciton optics, dynamics, and transport in atomically thin semiconductors, APL Mater., 2022, 10, 100701.
15. A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, Y. Li, B. Aslan, D. R. Reichman, M. S. Hybertsen and T. F. Heinz, Exciton binding energy and nonhydrogenic Rydberg series in monolayer WS 2, Phys. Rev. Lett., 2014, 113, 076802.
16. G. Zhang, A. Chaves, S. Huang, F. Wang, Q. Xing, T. Low and H. Yan, Determination of layer-dependent exciton binding energies in few-layer black phosphorus, Sci. Adv., 2018, 4, eaap9977.
17. E. Jung, J. C. Park, Y.-S. Seo, J.-H. Kim, J. Hwang and Y. H. Lee, Unusually large exciton binding energy in multilayered 2H-MoTe2, Sci. Rep., 2022, 12, 4543.
18. X. Zheng and X. Zhang, Excitons in two-dimensional materials, in Advances in Condensed-Matter and Materials Physics-Rudimentary Research to Topical Technology, 2019.
19. M. Dvorak, S.-H. Wei and Z. Wu, Origin of the variation of exciton binding energy in semiconductors, Phys. Rev. Lett., 2013, 110, 016402.
20. G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand and B. Urbaszek, Colloquium: Excitons in atomically thin transition metal dichalcogenides, Rev. Mod. Phys., 2018, 90, 021001.
21. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Atomically thin MoS 2: a new direct-gap semiconductor, Phys. Rev. Lett., 2010, 105, 136805.
22. J. Klein, J. Wierzbowski, A. Regler, J. Becker, F. Heimbach, K. Muller, M. Kaniber and J. J. Finley, Stark effect spectroscopy of mono-and few-layer MoS2, Nano Lett., 2016, 16, 1554–1559.
23. J. S. Ross, S. Wu, H. Yu, N. J. Ghimire, A. M. Jones, G. Aivazian, J. Yan, D. G. Mandrus, D. Xiao and W. Yao, et al., Electrical control of neutral and charged excitons in a monolayer semiconductor, Nat. Commun., 2013, 4, 1474.
24. S. Yang, Y. Chen and C. Jiang, Strain engineering of two-dimensional materials: Methods, properties, and applications, InfoMat, 2021, 3, 397–420.
25. K. He, C. Poole, K. F. Mak and J. Shan, Experimental demonstration of continuous electronic structure tuning via strain in atomically thin MoS2, Nano Lett., 2013, 13, 2931–2936.
26. A. Castellanos-Gomez, R. Roldán, E. Cappelluti, M. Buscema, F. Guinea, H. S. Van Der Zant and G. A. Steele, Local strain engineering in atomically thin MoS2, Nano Lett., 2013, 13, 5361–5366.
27. A. Raja, A. Chaves, J. Yu, G. Arefe, H. M. Hill, A. F. Rigosi, T. C. Berkelbach, P. Nagler, C. Schüller and T. Korn, et al., Coulomb engineering of the bandgap and excitons in two-dimensional materials, Nat. Commun., 2017, 8, 15251.
28. Y. Lin, X. Ling, L. Yu, S. Huang, A. L. Hsu, Y.-H. Lee, J. Kong, M. S. Dresselhaus and T. Palacios, Dielectric screening of excitons and trions in single-layer MoS2, Nano Lett., 2014, 14, 5569–5576.
29. J. S. Ross, P. Klement, A. M. Jones, N. J. Ghimire, J. Yan, D. Mandrus, T. Taniguchi, K. Watanabe, K. Kitamura and W. Yao, et al., Electrically tunable excitonic light-emitting diodes based on monolayer WSe2 p–n junctions, Nat. Nanotechnol., 2014, 9, 268–272.
30. S. Das, D. Pandey, J. Thomas and T. Roy, The role of graphene and other 2D materials in solar photovoltaics, Adv. Mater., 2019, 31, 1802722.
31. S. Yu, X. Wu, Y. Wang, X. Guo and L. Tong, 2D materials for optical modulation: challenges and opportunities, Adv. Mater., 2017, 29, 1606128.
32. Z. Sun, A. Martinez and F. Wang, Optical modulators with 2D layered materials, Nat. Photonics, 2016, 10, 227–238.
33. C. Chakraborty, N. Vamivakas and D. Englund, Advances in quantum light emission from 2D materials, Nanophotonics, 2019, 8, 2017–2032.
34. J. J. Thompson, S. Brem, M. Verjans, R. Schmidt, S. M. de Vasconcellos, R. Bratschitsch and E. Malic, Anisotropic exciton diffusion in atomically-thin semiconductors, 2D Mater., 2022, 9, 025008.
35. M. N. Brunetti, O. L. Berman and R. Y. Kezerashvili, Optical properties of anisotropic excitons in phosphorene, Phys. Rev. B, 2019, 100, 155433.
36. B. Ghosh, P. Kumar, A. Thakur, Y. S. Chauhan, S. Bhowmick and A. Agarwal, Anisotropic plasmons, excitons, and electron energy loss spectroscopy of phosphorene, Phys. Rev. B, 2017, 96, 035422.
37. X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia, H. Zhao, H. Wang, L. Yang, X. Xu and F. Xia, Highly anisotropic and robust excitons in monolayer black phosphorus, Nat. Nanotechnol., 2015, 10, 517–521.
38. J. Lu, J. Yang, A. Carvalho, H. Liu, Y. Lu and C. H. Sow, Light–matter interactions in phosphorene, Acc. Chem. Res., 2016, 49, 1806–1815.
39. S. Brem, M. Selig, G. Berghaeuser and E. Malic, Exciton relaxation cascade in two-dimensional transition metal dichalcogenides, Sci. Rep., 2018, 8, 1–8.
40. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
41. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
42. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
43. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
44. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, Influence of the exchange screening parameter on the performance of screened hybrid functionals, J. Chem. Phys., 2006, 125, 224106.
45. A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt and N. Marzari, An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions, Comput. Phys. Commun., 2014, 185, 2309–2310.
46. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
47. A. Laturia, M. L. Van de Put and W. G. Vandenberghe, Dielectric properties of hexagonal boron nitride and transition metal dichalcogenides: from monolayer to bulk, npj 2D Mater. Appl., 2018, 2, 6.
48. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli and M. Cococcioni, et al., Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter, 2017, 29, 465901.
49. M. Selig, G. Berghäuser, A. Raja, P. Nagler, C. Schüller, T. F. Heinz, T. Korn, A. Chernikov, E. Malic and A. Knorr, Excitonic linewidth and coherence lifetime in monolayer transition metal dichalcogenides, Nat. Commun., 2016, 7, 13279.
50. M. Kira and S. W. Koch, Many-body correlations and excitonic effects in semiconductor spectroscopy, Prog. Quantum Electron., 2006, 30, 155–296.
51. H. Haug and S. W. Koch, Quantum theory of the optical and electronic properties of semiconductors, world scientific, 2009.
52. N. S. Rytova, Screened potential of a point charge in a thin film, arXiv, 2018 preprint, arXiv:1806.00976,  DOI:10.48550/arXiv.1806.00976.
53. L. Keldysh, Selected Papers of Leonid V Keldysh, World Scientific, 2024, pp. 155–158.
54. J. W. McPherson, J. Kim, A. Shanware, H. Mogul and J. Rodriguez, Trends in the ultimate breakdown strength of high dielectric-constant materials, IEEE Trans. Electron Devices, 2003, 50, 1771–1778.
55. A. Galiautdinov, Anisotropic Keldysh interaction, Phys. Lett. A, 2019, 383, 3167–3174.
56. G. Meneghini, S. Brem and E. Malic, Excitonic thermalization bottleneck in twisted TMD heterostructures, Nano Lett., 2024, 24, 4505–4511.
57. A. Brown and S. Rundqvist, Refinement of the crystal structure of black phosphorus, Acta Crystallogr., 1965, 19, 684–685.
58. A. Usman, M. A. Aly, H. Masenda, J. J. Thompson, S. M. Gunasekera, M. Mucha-Kruczyński, S. Brem, E. Malic and M. Koch, Enhanced excitonic features in an anisotropic ReS 2/WSe 2 heterostructure, Nanoscale, 2022, 14, 10851–10861.
59. J.-W. Jiang, B.-S. Wang, J.-S. Wang and H. S. Park, A review on the flexural mode of graphene: lattice dynamics, thermal conduction, thermal expansion, elasticity and nanomechanical resonance, J. Phys.: Condens. Matter, 2015, 27, 083001.
60. Y. Aierken, D. Çakır, C. Sevik and F. M. Peeters, Thermal properties of black and blue phosphorenes from a first-principles quasiharmonic approach, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 081408.
61. J. Yang, R. Xu, J. Pei, Y. W. Myint, F. Wang, Z. Wang, S. Zhang, Z. Yu and Y. Lu, Optical tuning of exciton and trion emissions in monolayer phosphorene, Light: Sci. Appl., 2015, 4, e312.
62. S. Koch, M. Kira, G. Khitrova and H. Gibbs, Semiconductor excitons in new light, Nat. Mater., 2006, 5, 523–531.
63. J. Xiao, M. Zhao, Y. Wang and X. Zhang, Excitons in atomically thin 2D semiconductors and their applications, Nanophotonics, 2017, 6, 1309–1328.
64. J. J. Thompson, S. Brem, H. Fang, C. Antón-Solanas, B. Han, H. Shan, S. P. Dash, W. Wieczorek, C. Schneider and E. Malic, Valley-exchange coupling probed by angle-resolved photoluminescence, Nanoscale Horiz., 2022, 7, 77–84.
65. K. S. Thygesen, Calculating excitons, plasmons, and quasiparticles in 2D materials and van der Waals heterostructures, 2D Mater., 2017, 4, 022004.
66. D. Y. Qiu, G. Cohen, D. Novichkova and S. Refaely-Abramson, Signatures of dimensionality and symmetry in exciton band structure: Consequences for exciton dynamics and transport, Nano Lett., 2021, 21, 7644–7650.
67. B. Han, S. Stephan, J. J. Thompson, M. Esmann, C. Antón-Solanas, H. Shan, N. Kunte, S. Brem, S. Tongay and C. Lienau, et al., Angle-and polarization-resolved luminescence from suspended and hexagonal boron nitride encapsulated MoSe2 monolayers, Optica, 2022, 9, 1190–1196.
68. S. Prodhan, S. Giannini, L. Wang and D. Beljonne, Long-range interactions boost singlet exciton diffusion in nanofibers of π-extended polymer chains, J. Phys. Chem. Lett., 2021, 12, 8188–8193.
69. T. Yu and M. Wu, Valley depolarization due to intervalley and intravalley electron-hole exchange interactions in monolayer MoS 2, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 205303.
70. M. M. Glazov, E. L. Ivchenko, G. Wang, T. Amand, X. Marie, B. Urbaszek and B. Liu, Spin and valley dynamics of excitons in transition metal dichalcogenide monolayers, Phys. Status Solidi B, 2015, 252, 2349–2362.
71. A. Fick, Ueber diffusion, Annalen der Physik, 1855, 170, 59–86.
72. M. Z. Bellus, Z. Yang, J. Hao, S. P. Lau and H. Zhao, Amorphous two-dimensional black phosphorus with exceptional photocarrier transport properties, 2D Mater., 2017, 4, 025063.

---