Optical Properties of Anisotropic Excitons in Phosphorene
Matthew N. Brunetti Oleg L. Berman Roman Ya. Kezerashvili
cond-mat.mes-hall











^1Physics Department,  New York City College of Technology

The City University of New York,
  300 Jay Street,   Brooklyn NY, 11201, USA 

^2The Graduate School and University Center

The City University of New York,
New York, NY 10016, USA 






  We study the eigenenergies and optical properties of both direct excitons in a phosphorene monolayer in different dielectric environments, and indirect excitons in heterostructures of phosphorene with hexagonal boron nitride.
  For these systems, we solve the 2D Schrödinger equation using the Rytova-Keldysh (RK) potential for direct, and both the RK and Coulomb potentials for indirect excitons.
  The results show that excitons formed from charge carriers with anisotropic effective mass exhibit enhanced (suppressed) optical absorption, compared to their 2D isotropic counterparts, under linearly polarized excitations along the crystal axis with relatively smaller (larger) effective carrier masses.
  This anisotropy leads to dramatically different excited states than the isotropic exciton.
  The direct exciton binding energy depends strongly on the dielectric environment, and shows good agreement with previously published data.
  For indirect excitons, the oscillator strength and absorption coefficient increase as the interlayer separation increases.
  The choice of RK or Coulomb potential does not significantly change the indirect exciton optical properties, but leads to significant differences in the binding energy for small interlayer separation.
  



Optical Properties of Anisotropic Excitons in Phosphorene
    Matthew N. Brunetti^1,2, Oleg L. Berman^1,2, and Roman Ya. Kezerashvili^1,2
    December 30, 2023
===============================================================================





§ INTRODUCTION


  The experimental discovery of graphene in 2004 <cit.> was a fascinating first glimpse into the world of two-dimensional (2D) materials – its exceptional mechanical, thermal, and electrical properties suggested a new paradigm of flexible, durable, and highly efficient 2D electronic devices.
  In 2010, when monolayers of both insulating hexagonal boron nitride (h-BN) <cit.> and semiconducting transition metal dichalcogenides (TMDCs) <cit.> were first exfoliated, research efforts towards the development of next-generation 2D devices accelerated.
  The discovery of monolayer (ML) black phosphorus, referred to as phosphorene, came a full decade after the advent of graphene when a flurry of publications in 2014 heralded the arrival of a new addition to the 2D materials universe <cit.>.
  Within the first year of its discovery, the number of publications on phosphorene grew tenfold <cit.>, an unprecedented rate of growth even within the rapidly expanding field of 2D materials research.
  The sudden shift in intense research focus towards phosphorene is clearly justified, due to phosphorene's unique properties which make it a promising candidate for a variety of applications unsuited to its 2D relatives.

  Perhaps the primary distinguishing feature of phosphorene is its highly corrugated crystal structure, where a single monolayer appears to be composed of two distinct planes of phosphorus atoms.
  Each atom bonds to three neighbors, two of which are in the same plane and one which occupies the opposite plane <cit.>.
  The in-plane and out-of-plane bonds are characterized by very different bond lengths and bond angles, which in turn leads to extreme anisotropy in its mechanical, thermal, and electronic properties <cit.>.
  This intrinsic structural anisotropy is expressed in nearly every property of phosphorene, manifesting itself in charge carrier effective masses <cit.> and mobilities <cit.>, DC conductivity <cit.>, Raman spectra <cit.>, optical absorption <cit.> and photoluminescence <cit.> spectra, and in its response to mechanical strain <cit.>.
  Furthermore, the rippled structure of phosphorene leads to a larger surface area, which makes it ideal <cit.> for a variety of environmental <cit.> and biomedical <cit.> sensing applications.
  In addition to the already unique and intriguing structure of black phosphorus, it was recently shown that there are four more allotropes of ML phosphorus, each exhibiting unique crystal structures and distinct material properties <cit.>.

  Unlike the TMDCs, which are indirect gap semiconductors for all but their ML forms <cit.>, phosphorene remains a direct gap semiconductor from its bulk form (≈ 0.3 eV) <cit.> down to a single ML (≈ 2 eV) <cit.>.
  In addition to its dependence on layer number <cit.>, the band gap, as well as many other properties, is also sensitive to the magnitude and direction of an applied mechanical strain <cit.>, giving researchers a variety of ways to tailor the electronic properties of phosphorene to suit a particular task.
  Due to its broadly tunable band gap, phosphorene has also been identified as a promising material for converting solar energy to chemical energy <cit.>.

  
  
  

  Characterizing the optical properties of phosphorene is not only important in the context of phosphorene's potential applications to optoelectronic devices, but is an essential tool in understanding its fundamental properties, e.g. its electronic band structure.
  Indeed, some of the first experimental studies of phosphorene measured its photoluminescence (PL) <cit.>, optical absorption <cit.>, and Raman <cit.> spectra, as well as photocurrent generation <cit.>.
  In particular, one group observed an “extraordinary” PL peak in a phosphorene bilayer sample <cit.>, and another found a similarly strong PL signal in ML phosphorene, centered at 1.45 eV <cit.>.
  It was quickly realized that these remarkable optical properties were due to strongly bound and optically active excitons within the phosphorene ML.

  
  
  
  
  Exciton binding energies in bulk semiconductors tend to be on the order of a few tens of meV due to strong dielectric screening, adversely affecting their stability at room temperature and restricting the energy ranges in which they are optically active.
  By contrast, excitons in 2D semiconductors exhibit dramatically increased binding energies compared to their bulk counterparts, due on the one hand to quantum confinement effects reducing the degrees of freedom and therefore the average kinetic energy of the system <cit.>, while on the other hand experiencing stronger electrostatic screening.

  Around the same time that experimentalists first observed evidence of strongly bound and optically active excitons in mono- and few-layer phosphorene, a number of ab-initio studies of the electronic band structure and excitonic properties of phosphorene were published.
  Using a variety of theoretical approaches and numerical methods, the binding energy of the direct exciton in a phosphorene ML was calculated to be between 0.7-0.8 eV <cit.>, while the exciton binding energy of ML phosphorene on an Si/SiO_2 substrate was calcuated <cit.> and measured <cit.> to be around 0.3 eV.
  Another ab-initio study <cit.> predicted that the exciton-forming optical transition centered at 1.45 eV in a freestanding phosphorene monolayer could absorb a staggering 15% of incident light, but only if the excitation was linearly polarized along the armchair crystal direction – the extreme anisotropy of phosphorene left the transition almost completely dark for light polarized perpendicular to the armchair crystal axis.

  While research on excitons in phosphorene has largely focused on direct excitons in ML phosphorene, the field has recently expanded to consider spatially separated excitons formed in heterostructures (HS) consisting of two phosphorene monolayers separated by few-layer insulating h-BN, abbreviated as PHP HS.
  
  In this configuration, the electron and hole occupy different parallel phosphorene monolayers and their recombination is suppressed by the tunneling barrier created by the dielectric separating the phosphorene.
  As a result, indirect excitons exhibit much longer lifetimes than their direct counterparts.
  Indirect excitons (also called dipolar excitons due to their intrinsic dipole moment) exhibit many of the same properties as direct excitons, but importantly, this intrinsic dipole moment creates a weak, repulsive exciton-exciton interaction.
  As a result of their enhanced binding energy, the typically small effective mass of charge carriers in semiconductors, and weakly repulsive inter-particle interactions, indirect excitons in e.g. the TMDCs have been identified by theorists as promising candidates for high-temperature Bose-Einstein condensation (BEC) and superfluidity <cit.>.
  By extension, indirect excitons in phosphorene recently attracted interest when it was proposed that they could exhibit directionally-dependent BEC and superfluidity <cit.>.
  

  In this work, we calculate the eigenfunctions and eigenenergies of (i) direct excitons in ML phosphorene and (ii) indirect excitons in a PHP HS.
  
  First, the Schrödinger equation for an interacting electron and hole with anisotropic effective masses is solved numerically, yielding the eigenenergies and corresponding eigenfunctions of the excitonic ground and excited states.
  We then use well-established methods <cit.> for analyzing inter-excitonic optical transitions to study the anisotropic exciton eigensystem, obtaining the optical transition energies, oscillator strengths and absorption coefficients of the transitions.

  This paper is organized as follows.
  In Sec. <ref>, we summarize the theoretical approach for solving the 2D Schrödinger equation of the electron-hole system with anisotropic effective masses.
  In Sec. <ref>, we present the theoretical framework for calculating the optical properties of excitons.
  We describe our computational approach and discuss the choice of input parameters in Sec. <ref>.
  The results of our calculations for direct and indirect excitons follow in Sec. <ref>.
  We compare the calculated properties of excitons in phosphorene to the properties of excitons in other 2D materials in Sec. <ref>.
  Our conclusions follow in Sec. <ref>.



§ EXCITONS WITH ANISOTROPIC EFFECTIVE MASS


  In order to analyze the optical properties of excitons in phosphorene, we must first calculate the eigenenergies and eigenfunctions of the exciton by solving the Schrödinger equation describing an interacting electron and hole with anisotropic effective masses.
  This in turn requires providing the material properties of phosphorene as input parameters to the Schrödinger equation, in particular the anisotropic effective carrier masses, the dielectric screening length of ML phosphorene, and the dielectric constant of the environment.
  Significant effort has already been dedicated to characterizing the electronic structure of phosphorene both experimentally <cit.> and using theoretical <cit.> and ab-initio techniques <cit.>.
  Importantly, these analyses have yielded, among other things, the anisotropic effective masses of both electrons and holes.
  Both the static dielectric constant <cit.> and ML thickness of phosphorene <cit.> are also known, which is important for characterizing the electrostatic interaction between the electron and hole.
  These parameters can be inserted directly into the Schrödinger equation describing the electron-hole system, enabling a straightforward solution of the anisotropic exciton eigensystem.
  We therefore present the quantum mechanical description of the electron and hole using the 2D Schrödinger equation in such a way that it can be applied to either direct or indirect excitons – further discussion on the formal differences between the two systems will be given as necessary.
  Finally, we note that the orthogonal crystal axes of phosphorene are referred to as the armchair and zigzag directions <cit.> – following the convention in the literature <cit.>, we associate the x- and y-axes with the armchair and zigzag directions, respectively.

  
  Within the effective mass approximation, the Hamiltonian for an interacting electron and hole with anisotropic effective mass, constrained to move in the plane of their respective monolayers, is given by:

  
    Ĥ_0 = - ħ^2/2( 1/m_e^x∂^2/∂ x_e^2 + 1/m_e^y∂^2/∂ y_e^2 + 1/m_h^x∂^2/∂ x_h^2 + 1/m_h^y∂^2/∂ y_h^2) + V( 𝐫_e - 𝐫_h ),

  where the m_i^j, j=x,y, i=e,h correspond to the effective mass of the electron or hole in the x or y direction, respectively, the positions of the electron and hole are given by 𝐫_i = (x_i,y_i,z_i), and V( 𝐫_e - 𝐫_𝐡) describes the electrostatic interaction between the electron and hole.
  Eq. (<ref>) can be used to treat both direct excitons (z_e - z_h = 0) and indirect excitons (z_e - z_h ≡ D), where the interlayer separation D = l_phos + N_BN l_BN is the distance between the middle of the phosphorene monolayers, l_phos and l_BN are the thicknesses of ML phosphorene and h-BN, respectively, and N_BN is the number of h-BN monolayers separating the phosphorene.
  For indirect excitons in a PHP HS, we consider the average z-position of the electron and hole to be in the middle of their respective phosphorene monolayers, since the narrow vertical confinement of each particle resembles a one-dimensional particle in a box.
  Therefore, the electron-hole separation D must account for (i) the thickness of one phosphorene ML, which can be pictured as being “split” between the upper half of the lower ML and the lower half of the upper ML, as well as (ii) the vertical distance between the phosphorene monolayers themselves due to the intervening h-BN.
  
  
  
  

  Applying the standard procedure for separation of variables in the two-body problem <cit.> to the anisotropic Hamiltonian (<ref>), we define the center-of-mass coordinate as 𝐑 = (X,Y), X = (m_e^x x_e + m_h^x x_h)/(m_e^x + m_h^x), Y = (m_e^y y_e + m_h^y y_h)/(m_e^y + m_h^y), and the relative separation between the electron and hole as 𝐫 = 𝐫_e - 𝐫_h = (x,y,D), x = x_e - x_h, y = y_e - y_h, D = z_e - z_h.
  After separation of variables in (<ref>), the Schrödinger equation for the relative motion of the electron-hole system is given by:

  
    [ -ħ^2/2 μ^x∂^2/∂ x^2 - ħ^2/2 μ^y∂^2/∂ y^2 + V( 𝐫) ] ψ( 𝐫) = E ψ( 𝐫),

  where μ^j = (m_e^j m_h^j)/(m_e^j + m_h^j), j=x,y, is the reduced mass of the exciton in the x and y directions and E and ψ( 𝐫) are the eigenenergies and eigenfunctions of the exciton, respectively.

  For the direct exciton, the electron-hole interaction V( 𝐫) is described by the Rytova-Keldysh (RK) potential <cit.>,

  
    V_RK( 𝐫) ≡ V_RK( r ) = -π k e^2/2 κρ_0[ H_0 ( r/ρ_0) - Y_0 ( r/ρ_0) ].

  In Eq. (<ref>), k = 9 × 10^9 N·m^2/C^2, r ≡|𝐫| = √(x^2 + y^2) (for the direct exciton in a phosphorene ML) or r = √(x^2 + y^2 + D^2) (for the indirect exciton in a PHP HS) is the magnitude of the relative electron-hole separation, κ = (ϵ_1 + ϵ_2)/2 describes the surrounding dielectric environment, where ϵ_1 and ϵ_2 correspond to the dielectric constants of the materials a) above and below the ML for the direct exciton, or b) between and surrounding the phosphorene monolayers for the indirect exciton in a PHP HS, H_0 and Y_0 are the Struve and Bessel functions of the second kind, respectively, and ρ_0 is the screening length, given by <cit.>:

  
    ρ_0 = 2 πχ_2D/κ,

  where χ_2D is the 2D polarizability, which can be calculated via ab-initio methods.

  The asymptotic behavior of the RK potential with respect to the interparticle separation r is given by:

  
    V_RK (r) = k e^2/κρ_0[ ln( r/2 ρ_0) + γ]	   	r ≪ρ_0 
    
    
    	  - k e^2/κ r   	r ≫ρ_0
    	.

  where γ is Euler's constant.

  
  
  
  

  

  It was determined theoretically <cit.> that the static dielectric constant of phosphorene is anisotropic, ϵ^x = 12.5, ϵ^y = 10.2, and a more recent ab-initio study <cit.> likewise found that the 2D polarizability χ_2D was anisotropic, χ_2D^x=0.42 nm, χ_2D^y = 0.397 nm.
  However, the authors of Ref. Rodin2014 found that if χ_2D^x ≈χ_2D^y, one can approximate the 2D polarizability as isotropic by taking the average χ_2D = (χ_2D^x + χ_2D^y)/2 without changing the results significantly, and we will employ the same approach here.

  While the RK potential has been applied to the study of indirect excitons in 2D heterostructures <cit.>, it is still common to use the Coulomb potential to study these systems.
  Therefore, for indirect excitons, we will solve the 2D Schrödinger equation using both the RK and Coulomb potentials and compare the results.
  The Coulomb potential describing the interaction between spatially separated electrons and holes in a PHP HS can be written as:

  
    V_C( 𝐫) ≡ V_C( r ) = -k e^2/ϵ' √(x^2 + y^2 + D^2),

  where the dielectric constant ϵ' takes the value of the environment, i.e. ϵ' = κ.

  Now, the eigenenergies E and corresponding eigenfunctions ψ(𝐫) of the exciton are obtained by solving the Schrödinger equation (<ref>) with the RK potential (<ref>) and D=0 for direct excitons, or for indirect excitons with either the RK (<ref>) or Coulomb (<ref>) potentials and D = l_phos + N_BN l_BN.



§ EXCITON OPTICAL ABSORPTION


  Calculations of the optical properties of excitons in phosphorene employ well-established methods <cit.> for modeling the response of atomic-like systems to an incident electromagnetic (EM) wave of frequency ω and polarization ê.
This approach was successfully used to study optical transitions in excitons in semiconductor quantum wells <cit.>, and has recently been applied to excitons in 2D materials <cit.>.
  The optical transition energy E_tr corresponds to the difference in energy between the initial state, ψ_i, and the final state, ψ_f, and must coincide with the energy of the incident EM wave, i.e. E_tr = E_f - E_i = ħω.

  
  The oscillator strength, f_0, is a dimensionless quantity which gives the relative strength of a particular optical transition.
  For the isotropic 2D exciton, f_0 is proportional to the exciton reduced mass μ and does not depend on the in-plane orientation of the linearly polarized EM wave <cit.>.
  Modifying the standard expression for f_0 to account for anisotropy, we consider and calculate two distinct oscillator strengths, f_0^j, which correspond to the oscillator strengths of optical transitions induced by linearly polarized light oriented along the x and y axes, respectively.
  The polarization-dependent oscillator strength is thus given by:

  
    f_0^j = 2 μ^j (E_f - E_i) |⟨ψ_f |ê|ψ_i ⟩|^2/ħ^2.


  In Eq. (<ref>), |ψ_i ⟩ and |ψ_f ⟩ are the wavefunctions of the initial and final states, respectively.
  The allowed and forbidden optical transitions for a particular polarization ê can be determined by calculating f_0^j, or more specifically by computing the dipole transition matrix element, |⟨ψ_f |ê|ψ_i ⟩|^2, which represents the overlap integral between the initial and final wavefunctions when the initial state interacts with an external electric dipole moment.
  The dipole transition matrix element is zero for forbidden transitions and non-zero if the transition is allowed.
  For allowed transitions, the oscillator strength is positive under photon absorption (E_f - E_i > 0) and negative under photon emission (E_f - E_i < 0).

  A theoretical study <cit.> of the eigenstates of excitons in phosphorene with the RK potential using both Gaussian and sinudoidal basis functions provides crucial insight into the allowed optical transitions of the anisotropic exciton.
  It is well known that for 2D-hydrogen-like systems with either the Coulomb (isotropic dielectric environment) or RK (thin semiconducting film in an inhomogeneous dielectric environment) potentials, the allowed optical transitions of linearly polarized light are strictly limited to those in which the angular momentum quantum number differs by 1 between the initial and final eigenstates <cit.>.
  For the anisotropic exciton, on the other hand, the authors of Ref. Rodin2014 found that linearly polarized light can induce a transition between any two states in which the symmetry of the eigenfunction along the polarization axis changes from even to odd, or vice versa.
  For example, the ground state eigenfunction is even along both the x and y axes, n_x = n_y = 0.
  Therefore, linearly polarized light along the y direction can induce a transition to any state which is odd with respect to y, that is, n_y' = n_y + 1,3,5,…; n_x' = n_x, while x-polarized light can likewise induce a transition from the ground state to the eigenstates characterized by n_x' = n_x + 1,3,5,…; n_y' = n_y.

  While the oscillator strength gives us insight into the relationship between the eigenstates of the system and their response to an external EM force, there is a related quantity, the absorption coefficient α(ω), which describes how strongly a particular material absorbs light of a given frequency due to an optical transition specified by f_0^j.
  Let us consider the attenuation of an EM wave propagating through a homogeneous material.
  The intensity I of an EM wave of frequency ω is a function of the propagation distance z, given by:

  
    I(z,ω) = I_0 e^- α(ω) z,

  where I_0 is the initial intensity of the wave.
  Eq. (<ref>) illustrates the physical meaning of the absorption coefficient α(ω), i.e. it is the reciprocal of the propagation distance z over which the intensity of the EM wave of frequency ω decreases by a factor e.
  The absorption coefficient of optical transitions in isotropic atomic systems is given by <cit.>

  
    α(ω) = ( ω/ω_0 cπ e^2/2 ϵ_0 √(κ)μn_X/l_eff f_0 ) ( (Γ/2)/(ω_0^2 - ω^2)^2 + (Γ/2)^2),

  where ω_0 = (E_f - E_i)/ħ is the Bohr angular frequency of the optical transition, c is the speed of light, n_X is the 2D concentration of excitons in the system, l_eff is the effective vertical spatial extent of the exciton wavefunction, and Γ is the full width at half maximum of the optical transition, often referred to as the broadening, line broadening, or damping.

  The fraction n_X/l_eff physically represents the 3D exciton denisty, which for direct excitons in ML phosphorene can be straightforwardly written as the 2D concentration divided by the ML thickness, l_eff = l_phos.
  Considering indirect excitons in a PHP HS with the same 2D concentration n_X, it follows that the 3D exciton density must be reduced compared to direct excitons, i.e., l_eff > l_phos.
  In this case we consider the electron and hole to be bound to their respective phosphorene monolayers, and that the wavefunctions of the electron and hole do not penetrate far into the surrounding h-BN, so that the indirect excitons are effective contained within the two phosphorene monolayers.
  As the EM wave passes through the PHP HS, it only interacts with the exciton in regions where the exciton wavefunction is appreciably non-zero, i.e. the interaction only occurs within the phosphorene monolayers themselves.
  Therefore, we use l_eff = 2 l_phos for indirect excitons.

  Summing Eq. (<ref>) over all possible optically induced transitions in a given material (that is, not restricted to excitonic transitions) yields the absorption spectrum, a thorough description of how strongly the material absorbs light of frequency ω.
  Since we only consider a very limited subset of all possible optical transitions in phosphorene, let us focus on the scenario where the energy of the incident excitation is equal to the energy of the transition given by f_0^j, ħω = E_f - E_i, which corresponds to a local maximum in the absorption spectrum α(ω):

  
    α^j(ω=ω_0) ≡α^j = ( π e^2/2 c ϵ_0 √(κ)μ^jn_X/l_eff f_0^j ) ( 2/Γ).

  Eq. (<ref>) can be used to characterize how strongly a particular optical transition absorbs the incident excitation.
  Additionally, the value of α^j can be used to conveniently compare the relative strengths of different excitonic transitions.

  As previously stated, the anisotropic absorption coefficient, α^j, describes the attenuation of an EM wave of frequency ω and polarization ê as a propagates an arbitrary distance z through a dielectric.
  2D materials, however, do not have arbitrary thickness – indeed, 2D materials are noteworthy precisely because each monolayer has a well-defined thickness.
  Therefore, it would be instructive to consider how the intensity of the incoming EM wave is reduced due to the wave propagating a distance which corresponds exactly to the thickness of the phosphorene ML(s) occupied by the excitons.
  Recalling Eq. (<ref>), we now define the polarization-dependent absorption factor as 𝒜^j≡ 1 - ( I(z = l_eff, ω = ω_0)/I_0 ), or:

  
    𝒜^j = 1 - exp[ - α^j l_eff] = 1 - exp[ - ( π e^2 n_X/2 c ϵ_0 √(κ)μ^j f_0^j ) ( 2/Γ) ].

  Eq. (<ref>) therefore gives the fractional decrease in the intensity of the EM wave as it propagates through one exciton layer (that is, one ML for direct excitons, or one PHP HS for indirect excitons), i.e. 𝒜 = 0.01 means that each exciton layer absorbs 1% of the incident EM wave.
  The absorption factor is particularly convenient when comparing absorption between direct and indirect excitons, or between excitons in 2D materials with different thicknesses.

  Simplified forms of Eqs. (<ref>) and (<ref>) are given by Eqs. (<ref>)- (<ref>) in Appendix <ref>.



§ COMPUTATIONAL APPROACH




 §.§ Discussion of input parameters used in numerical calculations


  Calculating the optical properties of the exciton using Eqs. (<ref>), (<ref>), and (<ref>) requires the excitonic eigenenergies E_n and eigenfunctions |ψ_n ⟩, which are obtained by solving the Schrödinger equation (<ref>) with either the RK potential (<ref>) (for both direct and indirect excitons) or the Coulomb potential (<ref>) (for indirect excitons only).
  The Schrödinger equation takes as input parameters the anisotropic exciton reduced masses μ^x and μ^y and either the uniform dielectric constant ϵ' for the indirect exciton with the Coulomb potential, or the average environmental dielectric constant κ and 2D polarizability χ_2D for either direct or indirect excitons with the RK potential.

  Numerical solution of the Schrödinger equation using the aforementioned interaction potentials and input parameters is performed using the finite element method (FEM), which yields N pairs of eigenenergies and eigenfunctions which are solutions to the Schrödinger equation, corresponding to the N most-strongly-bound states.
  The eigenenergies E_n and eigenfunctions ψ_n, along with the appropriate anisotropic reduced mass μ^j, are then used as input parameters to calculate the oscillator strength, f_0^j according to Eq. (<ref>).
  The corresponding polarization-dependent absorption coefficients α^j and absorption factors 𝒜^j can then be calculated using as inputs the oscillator strength, f_0^j, the anisotropic exciton reduced mass μ^j, the 2D exciton concentration n_X, the phosphorene ML thickness l, environmental dielectric constant κ (or ϵ' for the Coulomb potential), and the excitonic optical broadening Γ.



  The standard values of these input parameters are given in Table <ref> – unless otherwise noted, all subsequent results were obtained using these values.
  Whereas the sets of carrier masses μ_i, i=a,b,c,d, were straightforwardly taken from the corresponding references, some additional discussion of the other input parameters is necessary.

  The phosphorene ML thickness, l_phos, obtained via ab-initio calculations in Ref. Kumar2016, agrees well with theoretical results from other works, namely 0.53 nm <cit.> and 0.6 nm <cit.> – we note that these last two references also measured the ML thickness using atomic force microscopy (AFM), obtaining values of 0.85 and 0.7 nm, respectively, but the authors themselves note that AFM measurements tend to over-estimate ML thickness.
  The 2D polarizability, χ_2D, was calculated from first-principles in Ref. Rodin2014 and agrees well with the value of 0.38 nm, also obtained from first-principles in Ref. Prada2015.

  The 2D exciton concentration, n_X, differs from the previous quantities in that it is not a material property that can be definitively measured or calculated – instead, n_X depends mainly on the excitation intensity, that is, a high-intensity laser will excite a higher concentration of excitons than a low-intensity laser.
  Therefore, it is reasonable to expect that n_X can and will vary significantly between experimental configurations, and even from one trial to the next.
  Instead of exhaustively considering a wide range of possible values of n_X, we instead choose one value of n_X which is representative of a typical experiment to use throughout our calculations.
  One recent study <cit.> of exciton-exciton annihilation rates in a phosphorene ML found that exciton-exciton annhiliation becomes the dominant recombination mechanism (as opposed to e.g. thermal decomposition or radiative recombination) at an exciton concentration of about 6.1 × 10^16 m^-2.
  While this value is not representative of a typical experiment, it may still be helpful to consider as an upper bound.
  Lacking an appropriate result from experimental studies in phosphorene, we turn instead to excitons in TMDCs, where we find a reasonable value of n_X in Ref. You2015, which studied excitons in a WSe_2 ML.

  Based on previous optical studies of excitons in 2D materials, we will use two values of Γ.
  It is important to note that the Γ obtained from experimental measurements is, like n_X, dependent on several external factors, including but not limited to: the sample temperature; the presence of structural defects within the sample; and/or surface contaminants at either the substrate/monolayer interface or the monolayer/air interface.
  These confounding variables can significantly alter the observed optical properties of the material, especially the presence of defects and contaminants which may be difficult to identify, characterize, isolate, and prevent.

  In choosing a value for Γ, we therefore adopt a similar approach to our choice of a value for n_X – we will choose a value for Γ which is generally appropriate for the system in question, but need not correspond exactly to one particular observed value.
  Many experimental PL/absorption studies of excitons in phosphorene are conducted with the phosphorene ML placed on a substrate (typically SiO_2), while the opposite side of the ML is left exposed to the atmosphere.
  These studies all observed significant broadening of the excitonic emission/absorption peak, with reported values of 70 meV <cit.>, 100 meV <cit.>, and 150 meV <cit.>.
  Therefore, when calculating α and 𝒜 for FS or uncapped phosphorene (that is, on an SiO_2 or h-BN substrate), we will use the value Γ = (70 meV)/ħ≈ 10^14 s^-1.
  For phosphorene encapsulated by h-BN, we again turn to similar studies on the TMDCs, where large excitonic broadening was observed in uncapped TMDC samples at room temperature <cit.>, but encapsulating the TMDC with h-BN was found to drastically reduce the excitonic linewidths to their cryogenic limit <cit.>, Γ = (11 meV)/ħ≈ 10^13 s^-1.
  Results for α and 𝒜 for direct excitons will therefore be presented using two different values of Γ, depending on the dielectric environment, while for indirect excitons only Γ = 10^13 s^-1 will be used since only h-BN encapsulation of the PHP HS is considered.



§ RESULTS OF CALCULATIONS


  

 §.§ Direct Excitons


  In this Section we present the results of our calculations of the eigenenergies and optical properties of the direct exciton using the input parameters listed in Table <ref>, focusing in particular on how our results change depending on the four sets of masses μ_i and the four dielectric environments denoted by FS, SS, HS, and HE.
  The notation X[i,k] will be used as a shorthand for “the value of the quantity X calculated using the set μ_i in the dielectric environment k ∈[FS, SS, HS, HE]”, i.e. f_0^y[a,FS] means “the value of f_0^y in FS phosphorene calculated using μ_a.”
  For the input parameters that produce the minimum or maximum value of a particular quantity, a min or max subscript will be added to the corresponding parameter, i.e. E_b[d_max,k] means that μ_d yields the maximum value of E_b for the given k.
  The percent difference between the maximum and minimum values, with respect to e.g. the μ_i, of a particular quantity will be denoted with a % subscript on the parameter, i.e. E_b[i_%,k].
  Averaging a quantity over a set of parameters will be denoted by the subscript avg., as in E_b[i_avg.,k].
  If i or k have been previously established in context, either index may be omitted from the notation X[i,k].

  
  
  
  
  
  
  
  
  

  

  In Table <ref>, the eigenenergies of the direct exciton in a phosphorene ML are calculated for four different dielectric environments and for each of the four sets of μ_i from Table <ref>.
  The binding energies follow the relation E_b[d,k] > E_b[c,k] > E_b[b,k] > E_b[a,k] for all dielectric environments k.
  The choice of the μ_i does not change E_b significantly, though the difference between the minimum and maximum E_b increases with κ, e.g. we find E_b[i_%,FS] ≈ 5%, while E_b[i_%,HE] ≈ 8.5%.
  This percent difference also increases for higher excited states.
  Of course, the parameter which most significantly changes E_b is κ, where E_b[i,FS] ≈ 2 E_b[i,(SS,HS)] and E_b[i,SS] ≈ 2 E_b[i,HE].

  In addition, our results for the binding energy of the direct exciton shown in Table <ref> agree very well with previously reported results in a variety of different dielectric environments.
  The binding energy of the direct exciton in a FS phosphorene ML has been calculated via ab-initio methods on several occasions – prior calculations vary between 700-850 meV <cit.>, which agrees quite well with our average value of about 740 meV, considering the previous results were obtained using a variety of methods and, therefore, a variety of input parameters.
  Regarding phosphorene on an SiO_2 substrate, the direct exciton binding energy was theoretically calculated to be about 400 meV <cit.>, which is within the range shown in Table <ref>.
  Another experimental study of phosphorene on an SiO_2 substrate <cit.> determined the binding energy to be about 300 meV, while a separate experimental investigation performed around the same time <cit.> obtained a surprisingly high value of 900 meV.
  Finally, in Ref. Rodin2014, where the electron-hole interaction was also modeled using the RK potential, the direct exciton binding energy was calculated to be about 200 meV for κ = 5, which falls within our calculated range for HE (κ = 4.89).

  

  The oscillator strengths of the first allowed optical transitions for x- and y-polarized light are shown for all four dielectric environments in Fig. <ref>.
  In particular, f_0^y (shown with the open markers) refers to the ( 0,0 ) → (0,1) transition and f_0^x (shown by the solid markers) refers to the (0,0)→(1,0) transition.
  Since both E_tr^j and f_0^j depend on μ_i, both quantities are averaged across the four μ_i and the average value is denoted by the plot marker.
  The major and minor axes of the ellipses encircling each data point mark the minimum and maximum values of E_tr^j and f_0^j.

  From Fig. <ref>, we see that both f_0^j and E_tr^j are decreasing functions of κ.
  The effect of anisotropy is also evident in the relative magnitudes of f_0^x and f_0^y, where f_0^y[i,k] ≫ f_0^x[i,k] for all k, and furthermore, f_0^y[i,HE] > f_0^x[i,FS]

  

  Table <ref> shows the ratios f_0^j[i_avg.,k] along with the corresponding C_D[k].
  Following the procedure outlined in Appendix <ref>, the average absorption coefficient with respect to the μ_i can be easily calculated as α^j[i_avg.,k] = C_D[i_avg.,k] f_0^j[i_avg.,k].

  Interestingly, whereas f_0^y[i,k] > f_0^x[i,k] for all i and k, we find that the opposite is true for α^x and α^y.
  The reason for this can be seen from Tables <ref> and <ref> – although f_0^y can exceed f_0^x by anywhere between about 30% (in FS) and over 100% (in HE), μ_i^x can be more than an order of magnitude larger than μ_i^y, so that the ratio f_0^x is always greater than f_0^y, and hence α^x will always be larger than α^y.

  A prior study <cit.> of the optical absorption and PL properties of ML phosphorene found that exciton-forming excitations were much more strongly absorbed if the excitation was polarized along x than along y.
  While the underlying theory of exciton-forming transitions differs substantially from the treatment of intra-excitonic transitions, it is plausible that in both cases, the fact that μ^x is much smaller than μ^y leads to enhanced absorption of x-polarized light.
  In other words, the amplitude of the oscillatory response of the exciton to an x-polarized driving force is much larger than the amplitude of oscillations induced by a y-polarized excitation, due to the fact that μ^y≫μ^x.

  From Eqs. (<ref>) and (<ref>), we see that α^j and 𝒜^j are inversely proportional to √(κ).
  At the same time, we find that f_0^y[i,FS] ≈ 1.25 f_0^y[i,HE], while f_0^x[i, FS] > 2 f_0^x[i, HE].
  Combining these two trends and assuming for the moment that Γ remains the same in all dielectric environments, we would expect a significant change in both α^x and α^y, approximately α^x[i,FS] ≈ 5 α^x[i,HE] and α^y[i,FS] ≈ 3 α^y[i,HE].
  Assuming instead that h-BN encapsulation significantly reduces Γ, we find that α^x[HE] is about twice as large as α^x[FS], while α^y[HE] is greater than α^y[FS] by nearly a factor of four.
  The aborption factor 𝒜^j reveals the significant difference in excitonic optical activity between the two polarization directions and the four dielectric environments, where we obtain 𝒜^x[i_avg.,k] as 1.36%, 0.68%, 0.55%, 3.0% and 𝒜^y[i_avg.,k] as 0.017%, 0.011%, 0.009%, 0.066%, for k = FS, SS, HS, HE, respectively.

  Comparing the optical quantities E_tr^j, f_0^j, α^j, and 𝒜^j for all μ_i and across all dielectric environments, we observe some general trends.
  Let us first address the quantities related to the y-transitions, followed by the x-related quantities.

  The optical transition energies follow the relation E_tr^y[c,k] > E_tr^y[d,k] > E_tr^y[b,k] > E_tr^y[a,k] in all k.
  Curiously, the oscillator strengths in FS phosphorene are reversed compared to the transition energies, i.e. f_0^y[a,FS] > f_0^y[b,FS] > f_0^y[d,FS] > f_0^y[c,FS], though we instead obtain f_0^y[a,k'] > f_0^y[d,k'] > f_0^y[b,k'] > f_0^y[c,k'] for k' = (SS, HS, HE).
  The ordering with respect to the μ_i is reversed again for the f_0^y, and therefore for α^y and 𝒜^y as well, i.e. f_0^j[c,k] > f_0^j[b,k] > f_0^j[d,k] > f_0^j[a,k], for all k.
  Recalling from Table <ref> that μ_a^y > μ_d^y > μ_b^y > μ_c^y, it appears that f_0^j is an increasing function of μ^j, while E_tr^j is a decreasing function of μ^j.
  
  

  In contrast to the y-polarized quantities, whose relative magnitudes were constant across the four dielectric environments for E_tr^y and f_0^y but were inconsistent in f_0^y, we find that the ordering of the E_tr^x[i,k] is different for each k, while the relative magnitudes of both f_0^x and f_0^x are consistent for all k.
  The transition energies E_tr^x show significant variation between different environments, with the only constant being that E_tr^x[d] is always the largest value.
  Whereas E_tr^x[d,FS] > E_tr^x[a,FS] > E_tr^x[b,FS] > E_tr^x[c,FS], we find that the relative magnitude of E_tr^x[c] increases as the dielectric screening increases, while E_tr^x[a] decreases relative to the other values, such that E_tr^x[d,HE] > E_tr^x[c,HE] > E_tr^x[b,HE] > E_tr^x[a,HE].
  The oscillator strengths, on the other hand, follow the same order as the μ_i^x themselves, i.e. μ_c^x > μ_d^x > μ_b^x > μ_a^x, for all dielectric environments k.
  Additionally, for all k, the ordering of the f_0^x, α^x, and 𝒜^x is reversed with respect to the μ_i^x.

  These observations suggest that while the optical properties corresponding to a particular excitation polarization ê are primarily determined by the corresponding μ^j, these quantities also exhibit some dependence on the opposite μ^j' ≠ j, stemming from the dependence of the optical properties on the excitonic ground state, whose properties must represent both μ^x and μ^y.
  Considering the uniquely strong response of the material properties of phosphorene, e.g. the anisotropic effective charge carrier masses, to external stimuli such as mechanical strain <cit.>, the preceeding analysis should prove useful in guiding future efforts to engineer phosphorene MLs with specific optical properties.



 §.§ Indirect Excitons


  In this Section we present and analyze the same calculated quantities as for the direct exciton, now for the indirect exciton in a PHP HS.
  All quantities were calculated by solving the Schrödinger equation with both the RK and Coulomb potentials, for different interlayer separations corresponding to N_BN = 1,2,3,…,8.
  We will adapt the notation X[i,k] used in the previous section to accomodate the different input parameters for the indirect exciton.
  Here, the notation X[i,p,N_BN] will denote “the quantity X calculated using μ_i, the potential p = (RK, C), and interlayer separation N_BN.”

  

  In Fig. <ref>, we plot the dependence of the indirect exciton eigenenergies E_n, n=1,2,3,6, on the interlayer separation, N_BN, where all E_n were calculated for both the RK (solid markers) and Coulomb (open markers) interaction potentials.
  Calculations were performed for all four μ_i.
  The plot marker denotes the average value of the μ_i, while the boundaries of the shaded regions denote the minimum and maximum values.

  The difference between the RK and Coulomb potentials is significant only for the first couple eigenstates at small interlayer separations.
  We find that the percent difference between the RK and Coulomb potentials decreases as N_BN increases, that is, E_b[i, p_%, N_BN] ≈ 11%, 7.7%, 5.5%,…, 2% for N_BN = 1, 2, 3,…, 8.
  The excited state energies follow a similar trend, where we find E_2[i, p_%, (1,8)] ≈ (7.5%, 1.5%).
  As shown in Eq. (<ref>) when the relative separation |𝐫| exceeds the screening length, ρ_0 = (2 πχ_2D)/κ, the RK potential converges to the Coulomb potential.
  For a PHP HS with κ = 4.89 and χ_2D = 0.41 nm, we calculate ρ_0 = 0.526 nm.
  Therefore, one would expect the RK and Coulomb potentials to converge as the total electron-hole separation exceeds 0.526 nm.
  Considering that l_BN = 0.333 nm, it is unsurprising that the indirect exciton binding energies for the RK and Coulomb potentials start to overlap as N_BN>2.
  The convergence of the excited state eigenenergies is also the result of increasing electron-hole separation, since the average separation of a two-particle bound state increases as progressively higher excited states are accessed.

  As with the direct exciton, the choice of μ_i does not significantly change the indirect exciton binding energy – for example, we calculate E_b[i_%, p, 1] ≈ 6%, decreasing to about E_b[i_%, p, 8] ≈ 4%.
  Although the value of κ is the same, the indirect exciton binding energy is reduced by about 40% compared to the direct exciton in HE due to the increased electron-hole separation in the PHP HS, from E_b[i,HE] ≈ 200 meV to E_b[i, p, 1] ≈ 120 meV.

  

  
  In Fig. <ref>, we present our calculations of f_0^j for both the RK and Coulomb potentials for N_BN = 1 - 6.
  The f_0^j are shown in separate plots for each value of N_BN.
  Our calculations show that increasing N_BN leads to an increase in f_0^j and a decrease in E_tr^j.
  Similar numerical studies of the indirect exciton in Xene <cit.> and TMDC <cit.> heterostructures with interlayer h-BN also indicated that f_0^j is an increasing function of N_BN.
  

  In general, f_0^y does not change much as N_BN increases because f_0^y was already quite large for the direct exciton.
  On the other hand, since f_0^x was small in the case of the direct exciton, we observe a significant increase in f_0^x as N_BN is incrementally increased.
  We also find that E_tr[i,C, N_BN] > E_tr[i,RK,N_BN], while f_0^j[i,RK,N_BN] > f_0^j[i,C,N_BN], for any i and N_BN.

  
  Let us also mention an unusual trend in the range of calculated f_0^x[i,p,N_BN] with respect to increasing N_BN – whereas the f_0^x[i] become more tightly clustered as N_BN increases from 1 to 5, the values become more spread out as N_BN continues to increase from 5 to 8.
  The relative magnitudes of the f_0^x[i] do not change with N_BN nor with p – they are always related by f_0^x[c] > f_0^x[d] > f_0^x[b] > f_0^x[a], as is the case with the direct exciton for k = (SS, HS, HE).
  Considering instead the incremental increase in f_0^x[i] with N_BN provides insight into the observed behavior.
  For example, the calculated values of f_0^x[a] increase nearly linearly at small N_BN before their growth is suddenly and strongly suppressed around N_BN = 6, i.e. f_0^x[a,RK,(1,2,3,4,5,6,7,8)] = (0.523,0.603,0.671,0.728,0.773,0.804,0.820,0.821).
  By contrast, f_0^x[c] increases nearly linearly for all N_BN, i.e. f_0^x[c,RK,2] - f_0^x[c,RK,1] = 0.067, while f_0^x[c,RK,8] - f_0^x[c,RK,7] = 0.057.
  For comparison, the change in f_0^x[b], which in general is only slightly larger than f_0^x[a], starts to taper off between N_BN = 7 and N_BN = 8, suggesting that it is also approaching some kind of asymptotic limit, while f_0^x[d], only slightly smaller than f_0^x[c], also shows nearly linear growth throughout the range of N_BN calculated here.

  We are therefore led to the conclusion that the incremental increase of f_0^x[c] must be strongly suppressed for N_BN≥ 8, such that f_0^x[c] approaches some constant value less than 1.
  On the other hand, it appears to be the case that the f_0^x[a] has already converged towards its asymptotic value, which must be close to the calculated value of 0.821 at N_BN = 8 – meanwhile, f_0^x[b] has an asymptotic maximum which is probably not much greater than f_0^x[b,RK,8] = 0.875.
  Recalling μ_c^x > μ_d^x > μ_b^x > μ_a^x, it appears that the magnitude of μ^j is directly related to this asymptotic value of f_0^j at large separations N_BN, and furthermore that the N_BN at which this asymptotic value is reached also increases with μ^j.

  Heterostructures of 2D materials exhibit a variety of interesting excitonic and optical behavior, including but not limited to the ability to tune the excitonic optical absorption strength and the corresponding transition energies.
  A more comprehensive study, one which, for example, systematically varies each input parameter individually over a broad yet physically plausible range, is necessary.
  By examining in detail how the excitonic properties change with respect to variations in the individual input parameters, we can deepen our understanding of which input parameters determine the maximal asymptotic value of f_0^j and the interlayer separation N_BN at which the asymptotic value is reached, and whether or not the asymptotic properties of f_0^j can be further tuned by strain, dielectric environment, external electromagnetic fields, etc., and if so, by how much these quantities may change when these external tuning mechanisms are applied.

  
  Let us now analyze in depth the effect that the choice of interaction potential has on the optical properties of the indirect exciton.

  Our calculations show that the percent differences for E_tr^j and f_0^j between the RK and C potentials generally decrease as N_BN increases.
  First, we find that E_tr^j[i,p_%,N_BN] > E_n[i,p_%,N_BN], for all i, j, n, and N_BN, i.e. the choice of interaction potential leads to a larger difference in the optical transition energies than in the corresponding individual eigenenergies, though in general the E_tr^j[i,p,N_BN] follow the same trends with respect to increasing N_BN as the eigenenergies themselves.
  These differences between the RK and Coulomb potentials decrease quickly with N_BN, from E_tr^j[i,p_%,1] ≈ (14.5%, 16.5%) for j = (x,y), respectively, to E_tr^j[i,p_%,5] ≈ (4.2%, 4.9%).

  Turning now to the oscillator strengths, we observe some unusual deviations from the consistent patterns observed for E_tr^j.
  First, let us discuss the general relationship between f_0^j, N_BN, and μ_i, returning later to the exceptions mentioned earlier.

  As mentioned earlier, since f_0^y is already quite large for the direct exciton, it does not change significantly as N_BN increases.
  By the same logic, the percent difference in f_0^y[p_%] is similarly small and decreases sharply as N_BN increases.
  In particular, f_0^y[p_%,(1,8)] ≈ (1%, 0.06%), and furthermore, f_0^y[c,p_%,N_BN] > f_0^y[b,p_%,N_BN] > f_0^y[a,p_%,N_BN] > f_0^y[d,p_%,N_BN], for all N_BN.
  By contrast, the relationship between f_0^x[i] and the interaction potential is less straightforward.
  Whereas f_0^x[(c,a),p_%,1] ≈ (6.4%, 8.87%), corresponding to the minimum and maximum values, we find unexpectedly that the relationship is reversed at large interlayer separations, i.e. f_0^x[(c,a),p_%,7] = (1.28%,0.19%).
  Furthermore, the percent difference f_0^x[a,p_%] actually increases from N_BN=7 to N_BN=8, from 0.19% to 0.25%.
  This is the only time that we observe an increase in the RK/C percent difference of any quantity with increasing N_BN.

  While these quantities may not be noteworthy on their own, they are analyzed in-depth here because of their unusual deviation from the trends which until now have consistently held true.
  It is unclear why only f_0^x shows this abnormal progression, even as E_tr^x and f_0^y do not.

  

  In Table <ref>, we present the calculated ratios f_0^j, averaged over the four μ_i, for all N_BN.
  We find that 𝒜^x[RK] increases from about 5.5% to about 8.8% as N_BN increases from 1 to 8, while 𝒜^y[RK] increases from about 0.81% to 0.87% across the same range of N_BN.
  The results for the Coulomb potential are very similar, especially for 𝒜^y and at larger N_BN for both x and y, though we find that 𝒜^x[C,1] ≈ 5.1%, nearly a 10% decrease compared to 𝒜^x[RK,1].


  



§ ANALYSIS AND DISCUSSION


  
  Now let us compare the properties of excitons in phosphorene to the properties of excitons in the TMDCs <cit.> and the buckled 2D allotropes of silicon (Si), germanium (Ge), and tin (Sn), known as silicene, germanene, and stanene, and collectively as the Xenes <cit.>.
  In Ref. Brunetti2018a, the properties of indirect excitons in a TMDC/h-BN heterostructure (THT HS) were calculated using a similar method to the one used in this work.
  The study focused on four of the most common TMDCs, namely MoS_2, MoSe_2, WS_2, and WSe_2.
  Calculations were performed for a THT HS with N_BN∈ [1,9] using only the RK potential.
  The relevant material parameters were the exciton reduced mass, μ, the 2D polarizability, χ_2D, and the TMDC ML thickness, l_T.
  Many ab-initio studies had previously calculated the material properties of the TMDCs, so for each material, two calculations were performed with two different values of the input parameters μ and χ_2D.
  Out of the multitude of possible choices, the input parameters were chosen based on the combination of values which produced the largest and smallest exciton binding energy, corresponding to upper and lower bounds on all calculated quantities.
  In particular, the smallest reported value of μ and the largest reported value of χ_2D were used together to provide the lower bound, while the largest μ and smallest χ_2D found in the literature were used to provide the upper bound.

  
  For indirect excitons in a THT HS, the binding energies were calculated to be between E_b[RK,N_BN=1] = (90-110, 100-105, 90-105, 90-110) meV in MoS_2, MoSe_2, WS_2, and WSe_2, respectively.
  Increasing the separation to N_BN = 5, the binding energies were reduced to between 50-70 meV for all materials, decreasing to about 40-55 meV at N_BN = 8.
  Comparing these values to the results shown in Fig. <ref>, we find that the binding energy of indirect excitons in a THT HS is smaller than in a PHP HS by about 5-10%.
  The 1s → 2p optical transition energy of the indirect exciton in a THT HS was calculated to be about E_tr[RK,(1,5,8)] = (50-60, 30, 20) meV.
  By comparison, Fig. <ref> demonstrates that the anisotropic exciton reduced mass causes E_tr^x (E_tr^y) to be significantly larger (smaller) than the analogous optical transition energy of the isotropic exciton.

  
  To facilitate the comparison of the optical properties of excitons in different materials, we use the absorption factor 𝒜^j to control for the factor of l, which is different for each 2D material, in the denominator of Eq. (<ref>).
  For a THT HS, the indirect exciton absorption factor was calculated to be 𝒜[RK,1] = 2-3.7%, while in a PHP HS, we calculate 𝒜^x[p,1] ≈ 5.1-5.5% and 𝒜^y[p,1] ≈ 0.81%, with the Coulomb potential yielding slightly smaller values of 𝒜 than the RK potential.
  The calculated values of 𝒜 in the TMDCs do not change significantly with increasing N_BN, reaching a maximum of about 2.5-4.4% at N_BN = 8, because f_0 is already quite large at N_BN = 1, similar to the observed behavior of f_0^y in a PHP HS.
  Again, we see here that the anisotropy of excitons in phosphorene leads to strongly enhanced (suppressed) optical activity under x- (y)-polarized excitations.

  
  Turning now to the properties of excitons in Xenes, we note that a direct comparison is complicated by the uniquely tunable nature of excitons in the Xenes.
  Briefly, the buckled crystal structure of the Xenes allows the band gap, and therefore, the effective mass of charge carriers, to be tuned by an external electric field, E_⊥, oriented perpendicular to the plane of the Xene ML.
  As a result, the excitonic properties can be dramatically altered by changing the magnitude of the applied electric field.
  We will restrict the discussion here to a range of E_⊥ which lead to binding energies that are comparable to excitons in phosphorene.
  Also, a previous ab-initio study predicted that the crystal structure of silicene became unstable around E_⊥≈ 2.6 V/Å, so all calculations in Ref. Brunetti2018b were performed for E_⊥≤ 2.7 V/Å.

  
  In Ref. Brunetti2018b, the properties of both direct and indirect excitons were calculated.
  For direct excitons, results were obtained for freestanding (FS) Xene monolayers and for Si monolayers encapsulated by h-BN.
  The properties of indirect excitons were calculated using both the RK and Coulomb potentials in Xene/h-BN heterostructures, primarily focusing on silicene (SHS HS).

  
  For the direct exciton in ML Xenes, it was calculated that E_b[Si,FS] ≈ 740 meV for E_⊥≈ 1.5 V/Å, E_b[Ge,FS] ≈ 740 meV at E_⊥ = 2.7 V/Å, while the binding energy in FS Sn reached a maximum of about 550 meV.
  Compared to E_b[HE] ≈ 200 meV in phosphorene, the direct exciton binding energy in HE Si reached a maximum of about 350 meV at the maximum electric field of E_⊥ = 2.7 V/Å, while the binding energy was about 200 meV for E_⊥≈ 0.8-1.2 V/Å.

  Originally, calculations of α and 𝒜 of excitons in the FS Xenes were performed for Γ = 10^13 s^-1, but for consistency we will instead assume Γ = 10^14 s^-1 as used here.
  Since the tuning mechanism of excitons in Xenes involves changing the charge carrier effective mass, the absorption coefficient and absorption factor are strongly suppressed at moderate to high electric fields, while the oscillator strength increases with increasing electric field.
  In general, when the electric field is large enough that the exciton binding energy is comparable to that of phosphorene, the value of 𝒜 in the FS Xenes is only about 1%, much weaker than 𝒜^x[FS] but comparable to 𝒜^y[FS].
  On the other hand, 𝒜[Si,HE] ≈ 2% at E_⊥ = 1.0 V/Å, while 𝒜^x[HE] ≈ 3% and 𝒜^y[HE] ≈ 0.6%, comparable to 𝒜 in the FS Xenes.

  For indirect excitons in an SHS HS, the maximum binding energy at E_⊥ = 2.7 V/Å was calculated to be about E_b[Si,p,1] ≈ 150 meV, not much bigger than the value of E_b[i,p,1] ≈ 120 meV shown in Fig. <ref>.
  The data for 𝒜 in an SHS HS again shows that indirect excitons are more optically active than direct excitons, where 𝒜[Si,p,1] ≈ 3-4% for E_⊥≈ 1 V/Å.
  Also, 𝒜[Si,p,N_BN] was found to depend only weakly on the choice of interaction potential p, while the change in 𝒜 with respect to N_BN is again quite small in the SHS HS, comparable to 𝒜^y.

  By comparing the properties of the anisotropic exciton in phosphorene to isotropic excitons in the TMDCs and Xenes, the effects of anisotropy are clearly emphasized.
  Whereas binding energies were mostly comparable in all three types of materials, the polarization-dependent optical properties of anisotropic excitons in phosphorene are drastically different from the optical properties of isotropic excitons in the TMDCs and Xenes.
  In particular, the small value of μ^x in phosphorene leads to a larger optical transition energy and significantly enhanced optical absorption, while the corresponding optical quantities for y-polarized excitations are much smaller than in isotropic excitons.



§ CONCLUSIONS



  We study the optical properties of direct excitons in ML phosphorene, and of indirect excitons in a PHP HS, by calculating the x- or y-linear-polarization-dependent optical transition energies, oscillator strengths, absorption coefficients, and absorption factors.
  To calculate these properties, the eigenenergies and eigenfunctions of the exciton were calculated by solving the Schrödinger equation using four different sets of anisotropic exciton reduced masses found in the literature.
  Additionally, we considered four different dielectric environments for the direct exciton corresponding to four common experimental (or theoretical, in the case of FS phosphorene) configurations.
  For the indirect exciton, the Schrödinger equation was solved using both the Rytova-Keldysh and Coulomb interaction potentials, and at different interlayer separations D corresponding to an integer number N_BN=1-8 of h-BN monolayers separating the ML phosphorene.
  Further analysis of our results for direct and indirect excitons was performed by examining how the results changed with respect to the change in exciton reduced mass, dielectric environment, choice of interaction potential, and interlayer separation.

  The intrinsic anisotropy of phosphorene manifests itself most noticeably in the optical properties of both direct and indirect excitons, where for direct excitons we predict that α^x > α^y by as much as a factor of 8, with this difference decreasing to about a factor of four for ML phosphorene encapsulated by h-BN.
  By combining the calculated absorption coefficient with the known thickness of ML phosphorene, we predict that direct excitons in a single phosphorene ML may absorb as much as 3% of incident x-polarized light, though this figure depends strongly on the 2D exciton concentration in the ML as well as on the line broadening of the excitonic transition.
  Analysis of the relationship between the absorption coefficient and the input parameters, and subsequent comparison to the optical properties of isotropic excitons in TMDCs and Xenes, suggests that the anisotropic mass is directly responsible for enhancing (suppressing) optical activity along the crystal axis with relatively light (heavy) exciton reduced mass.

  While exciton binding energies were comparable between the TMDCs, Xenes, and phosphorene, the excited states of the anisotropic exciton exhibit significant deviations from those of the isotropic exciton, where we find for the direct exciton that E_tr^x[FS] > E_tr^y[FS] by nearly a factor of two, with this difference decreasing as dielectric screening increases.
  The exciton binding energy also strongly depends on the dielectric environment, where we calculate direct exciton binding energies of about 800 meV, 350 meV, and 200 meV, corresponding to FS phosphorene, uncapped phosphorene on an SiO_2 or h-BN substrate, and ML phosphorene encapsulated by h-BN.
  Furthermore, we find excellent agreement between our calculated binding energies and previous theoretical and experimental results.

  The increased spatial separation of the electron and hole in a PHP HS leads to a significant reduction in the indirect exciton binding energy compared to the direct exciton in the same dielectric environment.
  Specifically, we obtain an indirect exciton binding energy of about 120 meV in an PHP HS separated by only one ML of h-BN, compared to a direct exciton binding energy of 200 meV in HE phosphorene.
  Whereas the binding energy of the indirect exciton is reduced due to the increased interparticle separation, we find that the optical activity of the indirect exciton is enhanced compared to the direct exciton, and furthermore, that the oscillator strength is an increasing function of interlayer distance.
  As a result, we predict that indirect excitons in a PHP HS can absorb up to 5% of an incident x-polarized excitation when separated by one ML of h-BN, increasing to more than 8% absorption for 8 layers of h-BN, though we again note that the specific values of these quantities depend heavily on external factors such as exciton concentration and exciton broadening.

  In general, analysis of our results shows that increased dielectric screening leads to a decrease in all calculated quantities, i.e. the eigenenergies E_n, oscillator strength f_0^j, absorption coefficient α^j, and absorption factor 𝒜^j.
  The calculated binding energies are not particularly sensitive to the choice of μ_i, but the optical properties can vary significantly depending on the relative magnitudes of μ^x and μ^y.
  In particular, our results indicate that the optical transition energies E_tr^j, absorption coefficients α^j, and absorption factors 𝒜^j are decreasing functions of the corresponding reduced mass μ^j, while the oscillator strength f_0^j is an increasing function of μ^j.
  While the dependence of the optical properties on the μ^j is not completely straightforward, it is clear that any mechanism which affects the anisotropic charge carrier masses in phosphorene will in turn affect the optical properties of excitons in phosphorene.
  Considering that phosphorene is interesting to researchers precisely because of the external tunability of its properties via e.g. mechanical strain, an exhaustive study of the dependence of the excitonic and optical properties on parameters such as the anisotropic reduced mass, Rytova-Keldysh screening length, and environmental dielectric constant would be a welcome contribution to the literature.

  Our results represent the first comprehensive numerical calculations of the eigenenergies and optical properties of indirect excitons in a PHP HS with up to 8 layers of h-BN.
  Furthermore, our calculations support experimental observations and theoretical studies of the direct exciton binding energy in ML phosphorene.
  We then expand upon these results by analyzing the dependence of the optical properties of excitons in phosphorene on a variety of common input parameters.
  Our analysis indicates that the excitonic optical properties are highly sensitive to the anisotropic effective carrier masses, which can be tuned experimentally.
  Finally, our results demonstrate that an exhaustive study of the eigenstates and optical properties of the anisotropic exciton, in particular the dependence of these quantities on the input parameters shown in Table <ref>, is warranted.

  

§ ACKNOWLEDGEMENTS


  The authors are grateful to acknowledge that this work is supported by the U.S. Department of Defense under Gran No. W911NF1810433.



  

§ CONVENIENT SIMPLIFICATIONS TO THE ANALYTICAL EXPRESSIONS FOR THE EXCITONIC OPTICAL QUANTITIES


  In this Appendix we simplify the analytical expressions for the optical properties α^j and 𝒜^j presented in Sec. <ref>.
  Examining Eqs. (<ref>), (<ref>), and (<ref>), we see that f_0^j depends directly on the numerically calculated eigenenergies and eigenfunctions, while α^j and 𝒜^j are given by purely analytical expressions, provided f_0^j is known.
  In other words, f_0^j is the only optical quantity that depends directly on the numerical results – on the other hand, l and κ are specified for a particular scenario, while n_X and Γ do not have specific values.
  As a result, α^j and 𝒜^j will exhibit the same qualitative behavior as the corresponding f_0^j and μ_i^j.
  By distinguishing f_0^j and the associated μ^j from the constants and input parameters l, κ, n_X, and Γ, we aim to provide the reader with a simple way to calculate α^j and 𝒜^j using different parameters than those given in Table <ref>.
  These quantities, which we call the scale factors and denote by C for the absorption coefficient α^j, and 𝒞 for the absorption factor 𝒜^j, act as a sort of conversion factor between the cumbersome but straightforward analytical expressions for α^j and 𝒜^j and the values of e.g. μ_i and f_0^j which are unique to our numerical results.

  Let us begin with Eq. (<ref>), and as a first step separate the physical constants from the input parameters:

  
    α^j     = C ( n_X f_0^j/√(κ)μ^j l_effΓ),
    	  
    
    	C     = π e^2/c m_0 ϵ_0 = 3.335 × 10^-5 m^2/s,

  where m_0 is the rest mass of the electron.

  Now, the fraction within brackets in Eq. (<ref>) contains all possible input parameters used in calculating α^j, but we can further refine our expression for C by recognizing that not every quantity shown in brackets in Eq. (<ref>) is a free parameter.
  In particular, we consider only four values for κ for the direct exciton and only one value for the indirect exciton as shown in Table <ref>.
  Similarly, we use only l_eff = l_phos and l_eff = 2 l_phos for the direct and indirect exciton, respectively.

  There are now five possible values of C which are applicable to our results:

  
    C_D    = C/√(κ) l_phos = 
    	  6.164 × 10^4 m/s,	   FS
    
    	  3.980 × 10^4 m/s,	   SS
    
    	  3.592 × 10^4 m/s,	   HS
    
    	  2.788 × 10^4 m/s,	   HE
    
    
    	  C_I    = C/√(κ) (2 l_phos)  = 1.394 × 10^4 m/s,

  where the subscripts D and I denote direct and indirect excitons, respectively.

  Now Eq. (<ref>) can be further simplified to:

  
    α_j = C_D/I( f_0^j/μ_i^j) ( n_X/Γ).


  Now, the foundation of our results, which consist of the numerically calculated eigenvalues and eigenfunctions, are effectively contained within the fraction f_0^j/μ_i^j, for which we will use the notational shorthand f_0^j≡ f_0^j/μ_i^j.
  On the other hand, n_X and Γ are essentially free parameters, insofar as the values given in Table <ref> are rough estimates meant to represent typical values of these quantities.
  Using the default values of n_X and Γ given in Table <ref>, we define the absorption coefficient scale factor as C_D/I = C_D/I( n_X/Γ) and obtain:

  
    C_D    = 
    	  3.082 × 10^6 m^-1,	   FS
    
    	  1.990 × 10^6 m^-1,	   SS
    
    	  1.796 × 10^6 m^-1,	   HS
    
    	  1.393 × 10^7 m^-1,	   HE
    
    C_I    = 6.969 × 10^6 m^-1.


  We also note that C_D/I is independent of the x- or y-polarization of the excitation.
  Using Eqs. (<ref>), (<ref>), or (<ref>), one can easily modify parameters such as n_X or Γ to match a particular scenario while still facilitating direct comparisons with the results presented in Sec. <ref>.

  Ultimately, the values of α^j presented in Sec. <ref> can therefore be calculated using the following expression:

  
    α^j = C_D/I( f_0^j).


  The calculation of the absorption factor 𝒜^j can likewise be simplified:

  
    𝒜^j = 1 - exp[ - 𝒞^D/I( f_0^j) ],

  where the dimensionless quantity 𝒞^D/I≡C_D/I l_eff is given by,
  
    𝒞^D    = 
    	  1.668 × 10^-3,		   FS
    
    	  1.076 × 10^-3,		   SS
    
    	  9.717 × 10^-4,		   HS
    
    	  7.541 × 10^-3,		   HE
    
    𝒞^I    = 7.541 × 10^-3.


  Finally we note that 𝒞^D/I≪f_0^j, so that the exponent in Eq. (<ref>) is always much smaller than unity.
  Applying the well-known expansion of e^x for small x, e^x≈ 1 + x + (x^2)/2 + …, the absorption factor can be approximated by:

  
    𝒜^j≈ 1 - (1 + (- 𝒞^D/I( f_0^j) ) + …) = 𝒞^D/I( f_0^j).

  This convenient approximation may prove useful for quickly estimating 𝒜^j from the values of 𝒞^D/I given above, along with the values of f_0^j and μ_i^j presented in Sec. <ref>.

  

  



§ ANALYSIS OF CALCULATED EIGENFUNCTIONS


  Since the Schrödinger equation (<ref>) features anisotropy along the x and y axes, the calculated eigenfunctions are most conveniently characterized by the quantum numbers n_x and n_y, in contrast to the isotropic case, where 2D polar coordinates are used and the eigenfunctions are described in terms of the principal and angular momentum quantum numbers, n and l, analogous to the 2D hydrogen atom <cit.>.
  The quantum numbers n_x and n_y corresponding to a particular eigenstate can be deduced by inspecting the eigenfunction and counting the number of times the eigenfunction changes sign along each axis.
  This is because eigenfunctions obey the empirical rule that the number of times the eigenfunction crosses ψ(𝐫) = 0 increases as the quantum number increases.

  When referring to the excitonic eigenstates in terms of the quantum numbers n_x and n_y, we use the notation (n_x,n_y), and denote the excitonic ground state by (0,0).
  Similarly, the pairs (1,0) and (0,1) refer to the eigenstates where the exciton has absorbed one quantum of energy in the x- or y-directions, respectively – we colloquially refer to these states as the “first excited state in  (x or y)”.

  Discussion of the eigenstates of the anisotropic exciton is further complicated by our computational method, which does not explicitly characterize the excitonic eigenstates in terms of the quantum numbers n_x and n_y, or indeed, in terms of any set of quantum numbers.
  Instead, our calculations yield only the eigenvalues and eigenfunctions, sorted by decreasing eigenenergy.
  In the process of analyzing and discussing our results, it may be instructive to refer to a particular eigenstate not in terms of the quantum numbers n_x and n_y, but by denoting it by its “rank” amongst all eigenstates produced by a particlar calculation.
  In this case we use the standard ket notation | n ⟩, that is, the ground state (which of course has the largest eigenenergy) is | 1 ⟩ and its eigenenergy is E_1, the first excited state (i.e. the state with the second-largest eigenenergy) is | 2 ⟩, with corresponding eigenenergy E_2, the second excited state (corresponding to the state with the third-largest eigenenergy, E_3) is | 3 ⟩, and so on.

  

  

  

  Curiously, we find that the ordering of the eigenstates with respect to the quantum numbers n_x and n_y is sensitive to the choice of μ_i.
  In other words, the anisotropic reduced masses μ^x and μ^y change the eigenenergy of the first x excited state (1,0) relative to the eigenenergies of the higher excited states in y, in particular the states (0,3) and (0,4) as shown in Table <ref> and in Fig. <ref>.
  Further analysis of the eigenfunctions corresponding to | n > 6 ⟩ confirms that the ordering of the eigenstates in terms of (n_x,n_y) differs depending on the relative magnitudes of μ^x and μ^y.

  Finally, let us mention that while solutions to the Schrödinger equation were obtained up to the | 12 ⟩ eigenstate, we have restricted our discussion and presentation of the results in the text to the first 6 eigenstates.
  Our reasons for this are fourfold: (i) to reduce visual clutter in the figures and emphasize the lower eigenstates which are more experimentally relevant; (ii) the optically active state (1,0) can appear as late as | 6 ⟩, so we do not truncate our results before this state; (iii) the eigenenergies of states | n > 6 ⟩ change by a very small amount, and their optical activity is very strongly suppressed due to the presence of other allowed optical transitions to states with | n ≤ 6 ⟩; (iv) eigenstates beyond approximately | 8 ⟩ strain our numerical methods and can occasionaly yield physically ambiguous or even non-sensical results.

  In Fig. <ref> we plot slices of the direct exciton eigenfunctions along the x- and y-axes.
  Also shown on the plots are the corresponding μ_i, the eigenstate | n ⟩, and the eigenenergy of the state, E_n.
  We find that μ_a and μ_b share the same internal structure, and so plots for μ_b are not shown.
  The associations between | n ⟩ and (n_x,n_y) shown in Table <ref> were determined by first examining the eigenfunctions (some of which are shown in Fig. <ref>) and counting the number of times the function changes sign along each axis, then cross-referencing those associations with the allowed and forbidden optical transitions of the anisotropic exciton, as determined theoretically in Ref. Rodin2014 and supported by our numerical results in Fig. <ref>.

  However, we note that the plots of the (0,2) and (0,4) states, corresponding to the states | 3 ⟩ (for all μ_i) and either | 5 ⟩ (for i=a, b) or | 6 ⟩ (for i=c, d), respectively, clearly show that the eigenfunction changes sign twice with respect to the x coordinate, suggesting that the states should have quantum number n_x = 2.
  Considering that these anomalous eigenstates appear before the (1,0) eigenstate for all μ_i, we conclude that this behavior is an aberration, and not to be interpreted as an appearance of a symmetric excited state in x (e.g. a state characterized by n_x = 2,4,6,…).
  These eigenstates are optically dark, so it is difficult to assess how the abnormal behavior of the eigenfunction would affect calculations of the optical properties related to these states, if at all.

  In Fig. <ref>, the allowed and forbidden optical transitions between the first six eigenstates are shown in blue for y-polarized excitations and in red for x-polarized excitations.
  Counting from the top-left of the plot, the row numbers denote the initial eigenstate | n_i ⟩, while the column numbers correspond to the final eigenstate, | n_f ⟩.
  Boxes lying to the right (left) of the diagonal thus correspond to optical absorption (emission) transitions, where the location of each box in the array, specified by the ordered pair of (row,column) numbers, corresponds to the initial and final eigenstates of the transition.
  Each box corresponds to a possible optical transition, and the color of the box is based on the result of calculating f_0^x and f_0^y.
  The box was colored red (blue) if f_0^x (f_0^y) was calculated to be non-zero, and was colored white if neither calculation returned a non-zero result.

  Due to intrinsic error both in the numerical eigenfunctions themselves and resulting from numerical integration of the dipole transition matrix element, the oscillator strength was "non-zero" if it was greater than 10^-4.
  The cutoff value of 10^-4 was chosen after computing f_0^j for all 12 calculated eigenstates and observing that the oscillator strengths of allowed transitions decreased by roughly an order of magnitude for each successive allowed transition from a given initial state.
  On the other hand, the numerical error in the calculated oscillator strengths for "dark" transitions was of the order of 10^-10 or smaller for small n_i and n_f, but reached as high as 10^-7 for transitions involving eigenstates n>8.

  Qualitatively, the allowed and forbidden optical absorption transitions shown in Fig. <ref> agree exactly with the theoretically predicted optical selection rules of Ref. Rodin2014.
  Apparently, the aforementioned anomalous eigenfunctions had no effect on the calculation of the optical selection rules using Eq. (<ref>).




