Inverse Funnel Effect of Excitons in Strained Black Phosphorus

We study the effects of strain on the properties and dynamics of Wannier excitons in monolayer (phosphorene) and few-layer black phosphorus (BP), a promising two-dimensional material for optoelectronic applications due to its high mobility, mechanical strength and strain-tuneable direct band gap. We compare the results to the case of molybdenum disulphide (MoS$_2$) monolayers. We find that the so-called funnel effect, i.e. the possibility of controlling exciton motion by means of inhomogeneous strains, is much stronger in few-layer BP than in MoS$_2$ monolayers and, crucially, is of opposite sign. Instead of excitons accumulating isotropically around regions of high tensile strain like in MoS$_2$, excitons in BP are pushed away from said regions. This \emph{inverse} funnel effect is moreover highly anisotropic, with much larger funnel distances along the armchair crystallographic direction, leading to a directional focusing of exciton flow. A strong inverse funnel effect could enable simpler designs of funnel solar cells, and offer new possibilities for the manipulation and harvesting of light.


I. INTRODUCTION
Two-dimensional (2D) crystals, such as graphene, transition metal dichalcogenides (TMDs) and, more recently, few-layer black phosphorus (BP), have revealed great technological potential thanks in particular to their unique combination of mechanical and optoelectronic properties. On the one hand, many of these atomically thin membranes can withstand unprecedented strains of up to 10−25% without plastically deforming or rupturing [1]. This is in stark contrast to most bulk semiconductors that fail mechanically at strains of about 0.1 − 0.4%. On the other hand, these materials cover a wide range of optically active (i.e. direct) band gaps. A particularly important 2D crystal in this regard is few-layer BP [2][3][4][5], as it is the only member of the family with a direct gap that covers the range between 0.3 eV and 1.8 eV as the number of layers is decreased. This is a crucial range of energies for many semiconductor technologies [6,7], including infrared photodetectors [8], telecommunications [9], and even photovoltaics [10][11][12], which could furthermore benefit from BP's high mobilities [13][14][15]. Finally, several of these 2D crystals exhibit an extraordinarily strong coupling between these two aspects, strain and optical activity [16]. The gap of TMDs monolayers, for example, decreases by up to ∼ 1.5% under 1% of uniaxial tension [17]. Once more, few-layer BP is remarkable in this regard, as it shows one of the strongest modulation of its band gap, ranging from a ∼ 6% (monolayer) [18] to a ∼ 23% (bulk) increase under 1% of uniaxial tension [19]. Not only is the strain sensitivity of the BP gap stronger and of opposite sign to that of TMDs, but it is also expected to be anisotropic with the direction of applied uniaxial strain [20]. This bestows few-layer BP with rather unique opportunities in the field of optoelectronics that are only recently beginning to be explored.
FIG. 1. Color online: An exciton funnel for MoS2 is represented in (a), where an indenter creates an inhomogeneous strain profile that modulates the gap (bottom), and pushes photogenerated excitons (in green) isotropically towards the center of indentation. In BP (b), the same strain profile creates a stronger, highly anisotropic inverse funnel effect that pushes excitons away from the indentation along the armchair direction.
A fundamental optoelectronic phenomenon in semiconductors is the generation and recombination of excitons, i.e. particle-hole pairs that become bound by Coulomb interaction and form a state with energy E ex inside the semiconductor gap E gap . When illuminated by a light source of frequency ω > E gap / , electrons are excited from the valence to the conduction band. The electron and the hole lose energy through several mechanisms [21], eventually relaxing to the exciton state in the gap. After a finite lifetime τ , typically longer than the preceeding relaxation, the exciton recombines through the emission of a photon of energy E ex < E gap , producing photoluminescence. A number of 2D crystals have remarkably strong excitonic photoluminescence [22][23][24][25][26], with binding energies E b = E gap − E ex typically exceeding those of semiconductors.
It has been proposed that strain-engineering of the gaps of 2D crystals could be used to efficiently manipulate excitons. By creating a strain gradient, e.g. by localized elastic indentation of the crystal, the exciton energy E ex is expected to vary spatially in a similar way as the gap itself. Feng et al. [1] predicted that strain gradients in MoS 2 monolayers create a force on (neutral) excitons that pushes them towards the regions of maximum tension (least gap) [27], in what was dubbed an "exciton funnel", see Fig. 1a. They argued that in a photovoltaic solar cell, the modulation of the optical gap through strain and the efficient funnelling of excitons to specific locations could lead to significant performance gains when compared to the standard photocarrier diffusion in conventional, fixed-gap cells, even beating the Shockley-Quessier limit [28]. Various aspects of exciton funnelling in MoS 2 have been explored experimentally [29,30]. More generally, a range of promising applications for strain-engineered optical properties in 2D crystals have been proposed [31][32][33].
In this work we study the properties of excitons in strained few-layer BP and their dynamics under strain gradients, and compare them to the case of MoS 2 monolayers. We find that BP exhibits a strongly anisotropic inverse funnel effect, whereby excitons are efficiently driven away from regions with tensile strains in a specific direction relative to the crystal axis. This behaviour is rare amongst known 2D crystals, and it could prove preferable to the original funnel effect of TMDs, by separating the source of strain and the location of exciton accumulation. An example is the inverse funnel solar cell of Fig. 1b, wherein the optically active regions are strained, and separated from the unstrained regions under the electrodes. We furthermore show that the absolute funnelling efficiency in few-layer BP is potentially far better than in MoS 2 monolayers, particularly as the number of layers increases.

II. FORMALISM
The starting point to analyse exciton formation in 2D crystals is an accurate description of the non-interacting, strain-free bandstructure around the gap. To this end we employ a tight-binding description (see Appendix A) carefully fitted to ab-initio calculations, both for fewlayer BP (Ref. [34]) and MoS 2 (Ref. [35]). These tightbinding models include p z phosphorus orbitals for BP, and p x,y,z sulfur orbitals plus the five d molybdenum orbitals for MoS 2 , see lattices in Figs. 2(a,b). The models can be extended to incorporate arbitrary strain profiles, an important advantage over ab-initio approaches. The (e) Schematic representation of the single-particle energy bands, with an exciton state at energy Eex marked in green. (f) Schematic representation of the particle-hole excitation spectrum. Excitons within the light cone E < c|k| can decay radiatively, and are thus 'bright'.
resulting bands (Fig. 2c,d), gap values and carrier effective masses are consistent with experimental and theoretical results available in the literature. We ignore spinorbit coupling, which is responsible in MoS 2 for a range of interesting spin-dependent phenomena, that are however not essential for the present work.
As expected, in the case of few-layer BP the gap from the tight-binding model is direct and sits at the Γ point. Its value ranges from E gap = 1.84 eV for the monolayer (also known as 'phosphorene', Fig. 2d) to E gap = 0.41 eV for bulk BP [7]. The MoS 2 monolayer gap in our model is around E gap = 1.82 eV [36], but sits at the K point (Fig.  2c), and becomes indirect for multilayer MoS 2 samples. As an indirect gap is much less active optically, many optoelectronic applications of MoS 2 are mostly restricted to the monolayer. Another important difference between the bandstructure of the two materials is BP's strong anisotropy of carrier effective masses. While the effective mass of MoS 2 for the conduction and valence band is isotropic due to lattice symmetry, m c,v x = m c,v y ≈ 0.5m e , in monolayer BP we have a high anisotropy both in the conduction (m c x ≈ 0.2m e , m c y ≈ 1.1m e ) and in the valence bands (m v x ≈ 0.2m e , m v y ≈ 3.9m e ). (Throughout this work, x refers to the armchair orientation, and y to zigzag, see Fig. 1.) Significant anisotropies persist as the number of layers increases, a consequence of BP's puckered lattice structure (Fig. 2b). It is crucial to take into account the mass anisotropy when discussing exciton formation, as effective masses directly control their binding energies and spatial dimensions.
The problem of anisotropic excitons has only recently been analysed [25,[37][38][39][40]. In the limit of low exciton density, it is possible to treat the Coulomb interaction between a single electron and a single hole as a two-body problem [41]. Here we use the theory of Ref. [39], where approximate analytical solutions, valid for anisotropic electron and hole masses, were derived. The wave function of an exciton of total momentum Q can be expressed as Ψ ex ( R, r) = e i Q· R φ( r), where R is the center-of-mass coordinate, and r is the relative coordinate. In the effective mass approximation, the total exciton energy disperses with wavevector Q as Fig. 2f, where the total masses are M x,y = m c x,y + m v x,y . The function φ( r) satisfies the Schrödinger equation, with reduced masses µ −1 x,y = 1/m c x,y + 1/m v x,y and Coulomb interaction V ( r) between the electron and the hole. Its solution yields the binding energy E b from which E ex = E gap − E b (see Fig. 2e). As shown in Ref. [42], the Coulomb interaction V ( r) in a thin slab of thickness d and with dielectric constant ε, embedded between two dielectric media with constants ε 1 and ε 2 , depends on the screening length r 0 = dε/(ε 1 + ε 2 ), which marks the crossover between a logarithmic divergence for r < r 0 and the usual 1/r behavior for r > r 0 . For anisotropic materials, one can approximate ε = (ε x ε y ε z ) 1/3 [43]. In a suspended BP monolayer (d = 5.24Å, ε = 10.3), the screening length is around r 0 ≈ 25Å. Since the estimated excitonic radii a x,y in BP monolayer are smaller than r 0 (see Table I), the Coulomb interaction will be dominated by the logarithmic part. This is also the case for a suspended MoS 2 monolayer (d = 6.14Å, ε = 18.8) [44], with r 0 ≈ 58Å. We assume this configuration, ε 1 = ε 2 = 1, throughout the rest of this work. From the variational approach in Ref. [39] the expressions for the exciton radii in x and y directions read and a y = λa x , where a 0 is the Bohr radius, and λ = (µ x /µ y ) 1/3 measures the mass anisotropy. The binding energy E b in the same ap- Due to the presence of the electromagnetic environment, an exciton is merely a quasibound state of finite lifetime. Its main decay channel is through the emission of a photon of energy equal to that of the exciton E( Q) and of wavevector k = (Q x , Q y , k z ) for some out of plane k z . Since the photon energy is E = c| k|, this constraint can only be satisfied if c| Q| E ex , i.e. for small momentum excitons within the narrow light cone depicted in Fig. 2f. These 'bright' excitons decay with a finite rate Γ Q (see Appendix B for a derivation and general expressions). Around Q = 0, the decay rate reads while Γ Q = 0 for 'dark' excitons outside the cone, within this particular decay channel. Here α ≈ 1/137 is the fine structure constant, φ(0) = 2/(πa x a y ) is the exciton wavefunction [39] at r = 0, and v cv i are the valence-conduction dipole matrix elements v cv x,y = ψ c (0)|∂ kx,y H( k)|ψ v (0) , where H( k) denotes the tightbinding Bloch Hamiltonian and ψ c,v (0) are its singleparticle eigenstates at either side of the (direct) gap.
Typical intrinsic lifetimes τ 0 around Q = 0 are very short, at around τ 0 = 30 fs for monolayer BP and τ 0 = 100 fs for monolayer MoS 2 . It is known from experiments [45] that the exciton lifetime dramatically increases with temperature, likely due to fast phonon-scattering of excitons into non-decaying 'dark' states, such as those depicted in Fig. 2f. A simple argument based on instantaneous thermalization has been proposed [46] that, generalized to anisotropic exciton masses, yields the following lifetime for temperatures higher than ∼ 0.1 K, This simple estimate predicts greatly enhanced τ ≈ 249 ps and τ ≈ 525 ps room-temperature lifetimes for BP and MoS 2 excitons, respectively, both within order-ofmagnitude range of experimental results in pristine samples [26,45,47,48]. (It should be noted that currently available experimental results for time-resolved exciton decay remain notoriously sample dependent, probably due to the effect of sample preparation, disorder, environmental screening, and the intrinsic complexity of outof-equilibrium exciton dynamics.) Table I summarises the above exciton properties for unstrained BP and MoS 2 monolayers. A generic strain field (x, y) = ij (x, y) (i, j = x, y, z) can be efficiently incorporated into the hopping amplitudes of our tight-binding model. We denote by t 0 α,α the = 0 hopping amplitudes between any two Wannier orbitals α, α sitting at positions r 0 α , r 0 α , and connected by vector αα | are the dimensionless local electron-phonon couplings [49], and r αα = r 0 αα + · r 0 αα are the inter-site vectors modified by the strain tensor at the bond. For simplicity we assume that β αα depend solely on the L 2 angular momentum of the α and α orbitals, not on their L z projections. Thus, BP has a single parameter which we take as β pp ≈ 4.5, while MoS 2 has three, β pp = 3, β pd = 4, β dd = 5 (the latter are consistent with the Wills-Harrison rule [50]). Due to a lack of  accurate estimates of the above parameters in the literature, these values have once more been chosen here on the basis of ab-initio calculations, specifically by matching the direct-to-indirect gap transitions under strain in monolayers (at -4% and 6.7% uniaxial in BP [18], and at 2-3% biaxial in MoS 2 [1,51]). Appendix C shows a comparison between the above theory and state-of-theart ab-initio calculations for the exciton binding energy, both as a function of biaxial strain and number of layers.
The strain-dependence of exciton radii a x,y , reduced masses µ x,y , binding energy E b , band gap E gap and exciton energy E ex are presented in Fig. 3, both for MoS 2 (left column) and BP monolayers (right column). In the former, although the gap remains direct, it is shifted slightly away from the K point as a result of the strain. The most notable difference between the two materials is the strong anisotropy, apparent in the exciton shape and masses (panels a-d) and the opposite trend of the band gap with uniaxial strain: decreasing for MoS 2 (panel g) and increasing for monolayer BP (panel h, also true in multilayers). Due to the almost strain-independent binding energy E b in both cases (panels e,f), the exciton energy E ex , in green, also behaves this way under increasing uniaxial strain. In the case of biaxial strain the effect is even more pronounced (see Appendix C). Thus, an exciton generated on a sample with a finite strain gradient will be accelerated towards regions with higher tensile strain in monolayer MoS 2 (funnel effect), or away from said regions in few-layer BP (inverse funnel effect), as depicted in Figs. 1.
The unusual sign of gap modulation with strain in BP (∂E gap /∂ ii > 0) as compared to transition metal dichalcogenides in general (∂E gap /∂ ii < 0), see Fig. 3h, has been demonstrated in optical absorption experiments [19] and ab-initio calculations [18,57]. It is ultimately a consequence of the puckered crystal structure of BP. The gap in this material, E gap ≈ 2t 2 + 4t 1 > 0, is controlled directly by the partial cancellation between out-of-plane t 2 and in-plane t 1 hoppings, which have opposite sign (t 1 < 0 and t 2 > 0, as defined in Fig. 5). Due to the lattice puckering, tensile strains in the plane suppress t 1 , but increase t 2 due to the positive out-of-plane Poisson ratio, leading to a rapid gap increase.
The exciton lifetime τ for the two monolayers is shown in Figs. 4(a-d) versus temperature and strain. Its strain dependence is visibly stronger in BP than in MoS 2 , even diverging at the strain-induced direct-to-indirect transitions ( xx = 6.7% and yy = −4%), at which the BP monolayer valence band mass vanishes. Consider next the maximum distance an exciton may be funnelled across before it decays. We assume that the exciton does not dissociate under the acceleration (type-I funnel [1]), which is the relevant regime for realistic strains in both these systems given their large binding energies. Take a perfectly ballistic sample with a linear spatial dependence of E ex ( r) = F · r produced by a strain gradient, [58]. If the sample is disordered or temperature is high, the Drude scattering or phase coherence time τ D due to defects or phonons may become shorter than the exciton's lifetime τ . Its propagation then becomes diffusive before decaying, and the travelled distance is reduced 59,60]. Figs. 4(e,f) show the ballistic funnel distances B i at T = 5K (τ ∼ 4 ps) traveled by an exciton generated at initial point r = (x 0 , y 0 ) under a linear uniaxial strain profile xx = 0 xx + g(x − x 0 ) or yy = 0 yy + g(y − y 0 ). We consider a small strain gradient g = 1% per µm, and plot B i as a function of initial strain 0 ii . In a MoS 2 monolayer B is isotropic and of the order of ∼ 70 nm at zero initial strain, always towards increasingly strained regions. In BP monolayer, the ballistic (inverse) funnel distance is instead highly anisotropic, reaching ∼ 440 nm along armchair and ∼ 20 nm along zigzag directions, always away from strained regions. As a result, exciton flow becomes focused along the armchair direction in sufficiently ballistic samples, a phenomenon that may once more be exploited in various optoelectronic applications, as it will boost the device performance for specific orientations of electrodes or terminals, such as in the solar cell of Fig. 1b. A key aspect for photocurrent generation [61] in a funnel solar cell is the efficiency of exciton dissociation at the harvesting regions. This will critically depend on the contact properties, in particular the band alignment between BP and the semiconducting electrodes and the quality of the contact. The electrode materials and configuration should be chosen so as to form a p-n junction at the contact that may tear the exciton apart, converting its energy to electrical power with optimal efficiency. A number of recent works have been devoted to the properties of contacts to BP [62], with a focus on photovoltaic efficiency [11,[63][64][65]. It has been predicted in particular that MoS 2 [63,64] or compressed BP itself [11] could be ideal electrode materials for BP-based solar cells.
The remarkable performance of the inverse funnel effect in BP monolayers is largely due to the small exciton mass M x along the armchair direction, see Table I, but also to the strong sensitivity of E gap and E ex with strain, Fig. 3h. The strain modulation of the binding energy E b gives a relatively minor correction, so that more complex bound states that are formed at high excitation regimes, such as biexcitons [66], are expected to respond to strain gradients in a similar way as excitons, albeit possibly with reduced lifetimes at high densities [47]. The efficient modulation of optoelectronic properties with strain, a hallmark property of this material, was recently showcased by optical absorption experiments in elastically rippled few-layer BP [19]. Increasing the number of BP layers, moreover, the inverse funnel performance is expected to improve even further. As the gap is reduced, E ex shifts down to energies with a far smaller photon density (the photon density of states is ρ(E) = 8πE 2 /(hc) 3 ), and the range of bright excitons shrinks. This produces a sharp increase of exciton lifetimes, see Table I. Moreover, while the strain sensitivity of the exciton energy ∂ E ex remains mostly unchanged, the averaged exciton mass decreases by up to ∼ 40%, which conspires to increase the funnel distance even further as the number of layers increases (more details on multilayer funnelling can be found in Appendix D). As an example, a ballistic three-layer BP sample is expected to reach values of B x in the tens of micrometers at T = 5 K. A real BP trilayer would obviously be in the diffusive funnel regime in this case, and additional decay channels may also have to be considered [67], but even with a τ D ∼ 1ps, one would expect an D x of several micrometers. This renders few-layer BP a far more promising platform for exciton funnelling than MoS 2 .
To conclude, we have characterised the properties of Wannier excitons in few-layer BP and MoS 2 monolayers under strain. We have shown that the former presents strongly anisotropic exciton properties and a high sensitivity to strain. As a result we have demonstrated that few-layer BP should exhibit a remarkably strong anisotropic inverse funnel effect, which could be exploited for a number of optoelectronic technologies, such as high efficiency funnel solar cells. While we have focused on the case of BP, our proposal of inverse exciton funnelling could potentially be realised in other, structurally similar compounds, such as group-IV monochalcogenides (e.g. GeSe) [68,69]. Some of these new materials (e.g. GeS or SnS), have multiple valleys and indirect gaps without strain, however. If they prove to be as strain-tuneable as BP, this feature could perhaps be turned into an advantage for exciton control. Assuming the indirect gap can be made direct under strain, one can envision strained, optically-active regions funnelling excitons towards dark, unstrained regions, which would result in enhanced lifetimes of accumulated excitons. Further opportunities to exploit the remarkable interplay between strain and exciton dynamics are also expected in twisted multilayers [70], and heterostructures combining several of these materials. Moiré patterns due to a lattice mismatch between layers are expected to give rise to gap modulations and spontaneous strain superlattices, phenomena already familiar from twisted graphene bilayers and graphene/boron nitride heterostructures [71][72][73][74][75]. Moiré patterns and strain superlattices could thus open the door to two-dimensional crystalline materials with built-in, spontaneous funnelling, without the need of externally induced strains. The computation of the exciton properties in the main text rely on the ability to obtain the gap, effective masses and dipole matrix elements v cv x,y of the different materials under study, with and without strain. To this end, we employ a tight-binding model fitted to ab-initio calculations, and incorporate strains as discussed in the main text. For BP we used the model by Rudenko et al. of Ref. [34], which considers 14 hopping parameters between the p z orbital at each phosphorus atom, see  Table II. The model includes hoppings between atoms with relative distance up to 5.49Å, and requires no crystal fields. It has been shown to accurately describe the bandstructure of BP from the monolayer to the bulk [19,34].
For MoS 2 monolayers we have used the model in Ref.  Table II for their values.   [35], which includes all p orbitals in sulfur atoms, and all d orbitals in molybdenum atoms. We have obtained the tight-binding parameters and crystal fields for this model by fitting to LDA bandstructure results (see footnote [27]). The resulting values are shown in Table III.

Appendix B: Radiative Lifetime of excitons in Black Phosphorus
Radiative recombination occurs when an exciton in its ground state |Ψ ex ( Q) , of energy E ex ( Q), recombines with the consequent emission of a single photon in the state |γ k,ν = a † ν ( k)|0 em with energy ω k = c| k|. The operator a † ν ( k) creates a photon with momentum k polarized along the vector e ν . In this section we present a detailed derivation of general expressions for the decay rate of an exciton within this single-photon channel using a generic two-dimensional tight-binding description for the system, and the general description of excitons of Ref. [39].

Model
Assume a 2D system with a generic Bloch Hamiltonian H( k) obtained e.g. from a tight-binding model. For simplicity, we further assume the system has a direct gap at the Γ point (the final expressions will still be valid for expansions around a different point, as long as we measure momenta from that point), so that at small | k| we may expand with H i 1 = ∂ ki H and H ij 2 = ∂ ki ∂ kj H, evaluated at k = 0. In the presence of an electromagnetic environment, the minimal coupling enters as H( k − e A), where A is the electromagnetic field. The expanded Hamiltonian becomes, to first order in A with the electromagnetic vertex defined as The second-quantization version of H( k) is obtained as usual, where s, s are indices in a basis in the unit cell, and are implicitly summed over. The electromagnetic vertex is similarly expressed as where the k j term in V i above gives rise to a −i∂ xj ψ s ( r).
The Hamiltonian of the electromagnetic environment, derived from quantizing the electromagnetic action S = 1 4 d 4 xF µν F µν , can be written as where ω k = c| k|, and 2ω k Ω 0 e ν a ν ( k)e i k r + h.c. .
In the above equations, Ω is the total volume of the system, and ν = ±1 are the two possible polarizations of the photon field, so that k · e ν = 0. If k = kẑ, for example, then e ± = (∓x − iŷ)/ √ 2 for a basis with circular polarization. This form of A corresponds to the Coulomb (or transverse) gauge ∇ · A = 0, for which A oscillates in the plane perpendicular to the propagation direction of the photon.
The complete vertex in second quantization then reads The notation k above stands for the photon wavevector within the sample plane. Note the implicit summation over i (dot product of gauge field and fermionic current).

Exciton ground state
We consider now an exciton, i.e. an electron-hole pair bound by Coulomb interaction and with energy within the gap. As discussed in Ref. [39], the wavefunction of the pair may be obtained from a Schrödinger equation in the electron-hole relative coordinate r, which has a reduced mass tensor (µ ij ) −1 = (m c ij ) −1 + (m v ij ) −1 , in terms of the effective mass tensors of the electron and the hole at the conduction and valence bands. The solution for the wavefunction is φ( r). An exciton with a total momentum Q has energy Fig. 2f in the main text, where E gap ( Q) is the gap of H( Q) and E b ( Q) is the exciton binding energy. The particle-hole state is written as where S is the surface of the system, R is the centerof-mass coordinate and |0 is the electronic system's ground state. We may Fourier transform the above using c † ( k) = 1 √ S d 2 r e i k r ψ † ( r) and its converse ψ † ( r) = 1 √ S 2D k e −i k r c † ( k). This gives, for the exciton state at momentum Q,

Exciton decay
We wish to find the relaxation rate of the |Ψ ex ( Q) exciton due to its coupling W em to the electromagnetic environment. According to the Fermi golden rule, this rate is Note that |γ k,ν = a † ν ( k)|0 em is a single photon state, and |0 em is the electromagnetic vacuum. We insert the form of W em , Eq. (B2), to get, where we define V k = (V x k , V y k , 0) and The matrix element in Eq. (B3) produces the Kronecker constraint δ k − Q δ s v δ sc , and finally gives where the star denotes complex conjugation and the parenthesis in the second line is denoted by F i for brevity. From Eq. (B1), the dipole matrix elements v cv i and w cv ij above are defined as first and second order derivatives of H( k) around the gap at k = 0, respectively, We have also used 1 S 2D k e −i k r = δ( r). We finally eliminate the sum over the photon's k with the above δ k − Q and obtain We use the energy constraint to perform the sum kz = Lz 2π dk z (where L z S = Ω), taking into account the appropriate Jacobian Recall that e ν (with its dual e * ν ) form a twodimensional and Q-dependent orthonormal basis of the plane orthogonal to the photon's k, namely k = (Q x , Q y , k (0) z ) in this case. Note also that Q = | Q| < E ex / c above, otherwise Γ Q = 0. If Q = 0 (exciton ground state), the photon has only z momentum, and e ν can be chosen asx,ŷ, i.e e i ν = δ ν,i . We then recover the expression given in the main text, Note that the w cv φ (0) term can be neglected if φ( r) is an even (differentiable) function of position. Likewise, the v cv φ(0) term vanishes for an odd φ( r), which is expected of an excited state for the exciton (not considered in this work). Note also that the above derivation relies heavily on momentum conservation laws, which apply only if the sample size is much larger than the photon wavelength λ = 2π/k = hc/E ex . In the opposite limit, the decay rate is expected to be linear in sample area, and the photon is emitted isotropically.

Appendix C: Theory benchmarking
In this section we assess the accuracy of our theory for the exciton binding energy E b , by comparing to published ab-initio results based on the Bethe-Salpeter equation, see Refs. [11,38,57,76]. Figure 6a shows the comparison of E b in unstrained BP multilayers for increasing number of layers n. We find good agreement up to around n = 3 layers, with our theory well within the dispersion of published predictions. Beyond n = 3, the two-dimensional Keldysh potential employed in this work leads to an overestimation of the binding energy. The reason is that when the multilayer thickness exceeds the vertical exciton radius, a proper three-dimensional calculation of E b becomes necessary, which in turn leads to smaller binding energies, properly captured by the abinitio results. This is a generic effect, whereby an increased dimensionality leads to reduced binding energies from confining potentials [77]. Figure 6b shows a second comparison of E b in BP monolayers as a function of biaxial strain. We once more find reasonable agreement between our theory and abinitio calculation. Note in particular that the trend of increasing E b with tensile strain is correctly captured. We thus conclude that our analytical description of exciton properties is quite accurate for BP multilayers up to n = 3 within a wide range of realistic deformations.

Appendix D: Inverse funnel effect in BP multilayers
The theory and simulations presented in the main text focus on the case of BP and MoS 2 monolayers. While for MoS 2 , multilayers are less useful for optoelectronics, given their indirect gap, this limitation does not apply to BP multilayers, whose gap is direct irrespective of the number of layers n. Indeed, as anticipated in the main text, we expect the inverse funnel effect in multilayer BP to strongly outperform that of BP monolayers. The main reason is a strong enhancement of exciton lifetime as n increases, which is a result of the decreasing bandgap (see Fig. 7a), and the correspondingly suppressed density of photon densities at smaller energies, ρ(E) ∼ E 2 . The sensitivity of the gap as a function of strain ∂E gap /∂ ii is not dramatically affected by n, see Fig. 7b. Our theory thus predicts a strong increase for the BP exciton lifetime as n grows, see Fig. 7c. Similarly, the corresponding ballistic drift length, Fig. 7d, is enhanced into the tens of micrometers for ten layers.
These multilayer predictions should be taken with cau- tion, however. The theory, as presented, is strictly valid for excitons confined to two dimensions. When the multilayer thickness exceeds the vertical exciton diameter, the theory ceases to apply, strictly speaking. This happens already for a BP trilayer. The inclusion of the vertical dimension into the theory, however, is expected to only modify the results for the exciton binding energy, see Appendix C, which is not essential to the funnel effect, as discussed in the main text. Indeed, the main driver for the efficient inverse funnel mechanism is the modulation of the single-particle bandgap with strain, Fig. 7b, which is correctly captured by our theory. Therefore, we expect the results in Fig. 7(c,d) to remain qualitatively correct in a more general three-dimensional theory.
A second limitation, however, is not so straightforward to generalise. The only exciton decay channel considered here is radiative decay. Our theory does not include other non-radiative decay channels, such as Auger scattering, that while relatively unimportant for monolayers, should dominate exciton decay for thicker multilayers. Therefore, the results of Fig. 7 should rather be interpreted as an upper bound, absent a more complete description of all exciton decay channels.
One last consideration that should be kept in mind is the problem of inducing strain gradients in multilayers. As n grows, the material becomes stiffer, which may challenge efforts to create a given strain gradient. Moreover, the appearance of interlayer shear becomes a possibility (not considered here) which would modulate the interlayer coupling spatially, adding considerable complexity to the problem. A complete study of the multilayer elasticity problem is far beyond the scope of the present paper, but should be carefully considered in multilayer experiments.
Despite the above considerations, we anticipate that the optimal number of layers for the inverse funnel effect is probably greater than one, but not much greater than eight, for which the bandgap begins to saturate to its bulk value. The performance for the optimal BP multilayer is expected, in any case, to greatly exceed the already remarkable results predicted here for the BP monolayer.