{'fraction_non_alphanumeric': 0.026035270422950337, 'fraction_numerical': 0.013214127818440832, 'mean_word_length': 4.3223529411764705, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective mass theory for the anisotropic exciton in 2D crystals: Application to phosphorene

We present a theoretical study of the exciton binding energy for anisotropic two-dimensional crystals. We obtain analytical expressions from variational wave functions in different limits of the screening length to exciton size ratio and compare them with numerical solutions, both variational and exact. As an example, we apply these results to phosphorene, a monolayer of black phosphorous. Aided by density functional theory calculations for the evaluation of the two-dimensional polarizability, our analytical solution for the exciton binding energy gives a result which is very close to the numerical ones and, in turn, is comparable to the experimental value, as recently reported.

Since the mechanical exfoliation of graphene 1 , research in understanding the properties of two-dimensional (2D) crystals has increased in many folds. Atomically thin single-layered materials obtained from transition metal dichalcogenides 2 , boron nitride 3 , bismuth 4 , etc. are being extensively studied for applications as electronic and photoelectronic devices. Few-layered black phosphorous (BP) is a recent addition to the list of graphene-inspired materials [5][6][7][8] . Apart from having a sizeable band gap which can be tuned by the manipulation of the number of layers, the atomic structure of BP is highly anisotropic which leads to high asymmetry of the electronic band structure even for few layers. In particular, a single layer of BP or phosphorene is attracting most of the attention. The peculiar anisotropic nature of the band gap distinguishes phosphorene from other 2D crystals, increasing its potential functionality.
Excitons are a bound state of an electron and an hole and play an important role in the optical properties of the material. Understanding the nature of excitons and their dependence on the electronic structure of the host material is critical and lends a deeper perspective into the many-body physics involved in 2D crystals. The 2D nature of the polarizability of these crystals introduces an important length scale (screening length) r 0 . For distances between charges in the crystal plane, r, greater than r 0 the electron-hole binding potential behaves like in a 3D system i.e., it goes as ∼ 1/r. However, for the case where r is less than r 0 the potential is 2D-like, i.e., logarithmic. This behavior makes excitons in 2D crystals different from their 3D counterparts 9,10 .
Since most common 2D crystals are isotropic, the effect of anisotropy on the optical properties of these materials has remained essentially unexplored. The appearance of phosphorene has, however, changed this view and recent works address this issue from an analytical 8 , numerical 11 , and first-principles 7 standpoints. Here we give a detailed account of a variational approach, introduced by us in Ref. 8, to the calculation of the exciton binding energy E X in anisotropic 2D crystals. Several analytical expressions are derived in certain limits of the 2D interaction potential. The accuracy of our analytical expressions for E X is tested against both variational and exact numerical solutions to the actual 2D potential, finding excellent agreement in a wide and experimentally relevant range of screening lengths. In particular, the value of E X for phosphorene, as obtained from our analytical expression, compares almost exactly to the numerical results. Furthermore, this value nicely agrees with the recently reported experimental result 12 .
The present work is divided as follows. In Sec. II we review the form and limiting behaviour of the Coulomb interaction potential for charged particles in 2D systems. In Sec. III we present our variational approach based on an anisotropic exciton wavefunction. We first present the analytical result for E X in the limiting case where the 2D potential reduces to the standard 3D Coulomb potential ∼ 1/r for isotropic 2D systems to later introduce the anisotropy and re-derive the binding energy for this case. In the same manner we derive analytical expressions for the isotropic and anisotropic binding energies in the opposite limit where the 2D interaction potential behaves logarithmically. We also compare our analytical expressions with the numerically solved variational problem as well as with the exact numerical solution. In Sec. IV we propose an alternative variational approach based on gaussian orbitals. In Sec. V, after computing the 2D polarizability with density functional theory (DFT), our analytical approach is applied to the case phosphorene. Finally we present our conclusions in Sec. VI.

II. BINDING PARTICLE-HOLE POTENTIALS IN 2D
As originally derived by Keldysh 13 , the Coulomb potential energy created by a point charge at the origin that electrons feel in 2D follows the expression: where r 0 ≡ dǫ/(ǫ 1 + ǫ 2 ) andǭ = (ǫ 1 + ǫ 2 )/2. Here d is the thickness of the 2D material, ǫ is its bulk dielectric constant, and ǫ 1 and ǫ 2 are the dielectric constants of the surrounding media, typically substrate and vacuum.
Here r 0 plays the role of a screening length and sets the boundary between two different behaviours of the potential. For r < r 0 the potential diverges logarithmically, as if created by line charges. In this limit, the potential takes the simplified form also given by Keldysh 13 : where γ is the Euler constant. For r > r 0 the potential becomes the standard Coulomb potential created by point charges which goes as 1/r: A very good approximation to the Keldysh potential, fairly accurate in both limits and simpler to use, was introduced by Cudazzo et al. 9 : (2))e −r/r0 .
(4) It is interesting to compare these four expressions as a function of the distance r in a range of several orders of magnitude both above and below r 0 . We present such a comparison in Fig. 1. There it can seen the range of validity of each approximation, the Cudazzo et al. expression being remarkably accurate for all distances.

III. VARIATIONAL WAVEFUNCTION APPROACH
For generic 2D crystals with electrons and holes presenting anisotropic effective masses, m φ(x, y) = 2 where λ is the variational anisotropy scaling factor relating the exciton extension along the x direction, a x (which is also a variational parameter) and the one along the y direction (a y = λa x ). With this variational wavefunction we can evaluate the expectation value of the kinetic energy: where µ x and µ y are the reduced effective masses, m e m h /(m e + m h ), along x and y directions, respectively. The expectation value of the potential energy is given by and the variational exciton binding energy is obtained from the addition of these two quantities, Upon minimization with respect to a x and λ, one obtains the optimal parameters defining the extension and shape of the exciton and the actual binding energy E X . Results from three minimization procedures, one analytical and two numerical, are presented in next section.

A. Analytical Results
The integral for the potential energy in Eq. (6) turns out to be too difficult for an exact variational analytical solution. The main goal of this section is to make use of the asymptotic behaviour of the Keldysh potential to get analytical expressions for E X in the limits r ≫ r 0 and r ≪ r 0 , namely, valid for large and small excitons, respectively.

r ≫ r0 limit
We begin by evaluating E X in the isotropic case (λ = 1, a x = a y = a), considering only the long-range behaviour of the Keldysh potential (see Eq. 3). The contribution of the potential energy to E X is given in this limit by Now minimizing E X (a) with respect to the variational exciton radius one obtains where the minimal exciton radiusã x is given bỹ and m and a 0 are the free electron mass and the Bohr radius, respectively. For the anisotropic case (λ = 1) the exciton extension along the x-direction is now given bỹ In the previous expression we find a function of λ defined through the elliptic integral Defining now µ xy as the exciton extension along the x axis can now be written asã The exciton extension along the y direction is thus and the λ-dependent binding energy of the exciton now becomes Defining now we can see that the minimal λ,λ, satisfies in general the following equation: which has no analytical solution for λ. However, it can be shown that for λ 1 Finally, notice that the results obtained in this subsection will be valid as long as the exciton extension in both x and y directions is much larger than r 0 . The consistency of this approximation for given experimental parameters (r 0 , µ x , µ y , andǭ) has to be checked a posteriori.

r ≪ r0 limit
As r → 0 the logarithmic behaviour of the Keldysh potential dominates. The potential energy in 2D takes now the form given in Eq. 2. In the isotropic case the exciton radius is now given bỹ and the binding energy of the exciton is For an anisotropic system the λ-dependent exciton extension along the x direction is given bỹ Using now a different definition for µ xy the exciton x-extension now becomes Again, taking into account that a y = λa x , the λdependent minimal exciton extension along the y direction isã Finally we obtain the exciton energy for this case:  where the minimal λ is Again notice that this result will be valid as long as the x and y minimal extensions of the excitonic wave function are small compared to r 0 . Finally, for completeness, we have also obtained an analytical expression for the binding energy using the Cudazzo potential in the isotropic case:

B. Numerical Optimization and exact solution
To validate and test the accuracy of the limiting analytical expressions given in the previous section, we now use the wavefunction in Eq. (5) to numerically compute the potential energy given by Eq. 6 for the exact Keldysh potential. We also solve, numerically as well, the 2D Schrödinger equation for the same potential, which will give us the exact value of E X (down to the required numerical precision). Exciton binding energies for the isotropic case (λ = 1) are presented in Fig. 2 as a function of the screening length r 0 . For comparison's sake, we takeǭ = 1, i.e., the 2D crystal is suspended in vacuum, and µ x = µ y = m/2. Thus, according to Eq. 9, E x = −2 Ryd for r 0 = 0. The numerical variational result compares very well with the exact numerical value in the large range of explored screening lengths. For r 0 ≫ a 0 the analytical solution in Eq. 21 works fairly well. There the size of the exciton is smaller than r 0 and the 1/r contribution to the Keldysh potential is negligible. As expected, the analytical solution starts to fail as r 0 → a 0 since there the size of the exciton becomes comparable to r 0 and the long-range 1/r contribution to the Keldysh potential becomes dominant. (One should keep in mind that the limit of validity of the analytical result, as shown in Eq. 20, logically depends on the values ofǭ and µ.) We also compare with the result given by Eq. 28, obtained using the approximate expression to the potential in Eq. 4. This expression, although not as friendly as the previous one, extends the limit of validity of our analytical results down to r 0 ≈ a 0 . The results for the anisotropic case are presented in Fig. 3 for r 0 = 10a 0 and r 0 = 200a 0 as a function of the anisotropy ratio µx µy (notice a difference of one order of magnitude in the energy scales of each plot). Once again there is close agreement between the analytical solution, the numerical optimization, and the exact numerical solution for large r 0 , while for the smaller value, the analytical solution visibly deviates from the other two.
An important prediction of our analytical results is the relation between the anisotropy in the exciton extension and the effective masses:λ = ax ay ∼ µx µy 1/3 , which becomes exact in the limit of small excitons. To test this relation we fitted the optimal value of the variational parameterλ to the law for a large range of r 0 . The results of this fit are presented in Fig. 4. They confirm our analytical results and recover the exact 1/3 exponent in the limit of large r 0 . We finally provide a comparison between the exact and variational wave functions for several values of r 0 in the isotropic limit (see Fig. 5). Note that the distance is rescaled with the optimal radius and the amplitude of the wave function with the normalization constant .This representation illustrates to what extent the exact and variational wave functions satisfy similar scaling relations. Note that at r = 0 the exact wave functions do not show the prominent cusp of a 1s Slater-type orbital. This softened behavior at the origin suggests than a combination of gaussian functions may capture more accurately this feature of the wave function, as shown in next section.

IV. GAUSSIAN-BASIS VARIATIONAL METHOD
We have found that, in the limit of very small r 0 the binding energy is very sensitive to small changes in r 0 . Furthermore, the numerical solution of the Schrödinger equation requires a very fine mesh to reproduce the bound state in such a limit. On the other hand the analytical result found in the 1/r limit of the potential constitutes an isolated point an thus cannot be easily extended to small but finite values of r 0 . It is therefore interesting to find an alternative numerical method to study anisotropic excitons in the limit of small r 0 . Moreover, as we have presented in Fig. 5, the behavior of the exciton exact wave functions for different values of the screening length r 0 resembles more a 1s Gaussian than a Slater-type orbital. Gaussian-type orbitals (GTO) are very efficient basis sets used intensively in quantum chemistry and solid state calculations.
The gaussian basis functions {χ p } follow the expression: The index p is an integer that here has been chosen to run from 1 to 4 and the exponents α are coefficients that have to be optimized to minimize the ground state energy obtained by the variational method. Unlike the conventional gaussian approach, the anisotropy of the problem introduces two coefficients α p x , α p y per GTO. We limit variational freedom, assuming that the anisotropy is identical for the four basis functions.
In this equation κ is a constant that does not depend on p so we reduce the number of exponents that we have to optimize from eight to four. The variational wave function in the GTO basis is given by For fixed values of α p , the energy is computed by generalized diagnonalization of a 4x4 matrix. The matrix elements of the kinetic energy are: The matrix elements of the potential energy are computed by numerical integration and the overlap matrix is where M pq = (α p x + α q x )(α p y + α q y ). The binding energy and the optimal wave function are obtained by minimizing numerically the energy with respect to the five variational parameters. Efficient optimization of the energy requires a careful choice of the initial guess for the values of the exponents α p . In our case, we choose the optimal values for a 1s orbital of a "2D hydrogen atom", taking µ x = µ y = 1 and r 0 approaching zero. Once these optimal exponents are obtained, r 0 is changed slightly and the problem is solved again, using this time the exponents α obtained in the previous step. The procedure continues until the ground state energy for r 0 = 10a 0 is reached.
Fixing r 0 = 10a 0 and µ x = m, the effective mass along the y axis µ y is varied from 1 to 40m. In this case the initial guess for the exponents is the last set of coefficients α obtained when changing r 0 . A comparison between the values of the ground state energy as function of µ y variationally obtained with those computed by numerical solution of the Schrödinger equation with the Keldysh potential can be seen in Fig. 6 In Fig. 7, the length scales a p x = 1/ √ α p x of the whole set of optimal GTO's are plotted vs µ y . κ is also plotted in Fig. 8 against µ y . In Fig. 9 we show a log-log representation of κ vs µ y where a linear fitting has been made.

V. APPLICATION TO PHOSPHORENE
A paramount example of an anisotropic 2D crystal is phosphorene 5,6,8 , where the effective masses along x and y directions can differ by even an order of magnitude. From the start we have chosen to express the 2D potential constant in terms of the bulk dielectric constant ǫ and the effective thickness of the 2D crystal d (see Eq. 1). Equivalent expressions for the 2D potential, which rely on the evaluation of the actual 2D polarizability of the 2D crystal, χ, have recently been proposed [9][10][11] . These are probably more appropriate for actual crystals, although it has also been shown that Eq. 1 works well as long as ǫ is taken as the in-plane component of the bulk dielectric tensor 10 of the 3D crystal. Here we will compare both possibilities.
As shown in Ref. 9, the screening length r 0 depends on the polarizability χ as r 0 ≡ 2πχ. The polarizability for 2D materials can be computed using the expression where L is the distance between layers in a 3D layered structure. As can be seen, the dielectric function ǫ tends to unity as the inter layer distance L tends to infinity. We have computed the dielectric function at different inter layer distances within the density functional theory framework using the PBE functional 15 and norm conserving TM pseudopotentials as available in the SIESTA package 16 . The atomic and electronic structures have been duly converged on all parameters. SIESTA calculates the imaginary part of dielectric function from which the real part of it is obtained using the Kramers-Kronig relations. In order to account for the under-estimated band gap, the scissor approximation, as implemented in SIESTA has been utilized. The scissor shift of 1.2505 eV was made to match our previously reported band gap value of 2.15 eV 8 .
Using the plane-averaged static dielectric function calculated with SIESTA (see Fig. 10), a value for χ of 3.9 A is obtained across different values of L (in particular as L → ∞ as shown in Fig. 10), which provides a compelling proof of the expected behaviour given by Eq. 36. This value for χ yields a screening length of r 0 = 24.5 A. From Eqs. 23, and 24, 25, and 27, and using the values for the reduced effective masses previously calculated by us 8 , we obtainλ = 0.52 and an exciton extension of a x = 14.4Å andã y = 7.5Å along x and y directions, respectively. Since these values are smaller than r 0 , the use of the analytical expressions obtained in the logarithmic limit of the potential is justified. Also this was expected from Fig. 2 and, in particular, from the comparison shown in Fig. 3 for anisotropic cases. There it can be seen that already for r 0 = 10a 0 the deviation between the analytical result and the numerical ones is less than 10% for a ratio µ y /µ x ≈ 7 (which corresponds to phosphorene). Using now Eq. 26 we obtain an exciton binding energy for phosphorene in vacuum of E X = 0.78 eV. This result is remarkably close to a recently reported value of ≈ 0.9 eV, measured for phosphorene on a SiO substrate 12 . Similarly, we can use the value of r 0 obtained from the real part of the bulk ǫ (at zero frequency) and the thickness d of the monolayer. Since the thickness is somewhat undetermined, we have computed the binding energy for varying d, as shown in Fig. 11. It can be observed that the binding energy of the monolayer computed using the microscopically derived χ matches the binding energy obtained using r 0 ≡ dǫ/(ǫ 1 +ǫ 2 ) for d ≈ 7Å, which certainly can be considered the thickness of phosphorene.

VI. CONCLUSIONS
We have shown that a variational approach to the exciton binding energy in anisotropic 2D crystals can give excellent results when compared to numerical approaches. Furthermore, we have studied the range of validity of analytical solutions to the variational approach and found that these can give highly satisfactory results in a range of values of screening lengths which is relevant for actual 2D crystals such as phosphorene. We have computed the exciton binding energy in this case and found a very good agreement with a recently reported experimental result 12 . As long as the screening length to exciton size is large, our analytical results can be trivially used to predict the exciton binding energy of any 2D crystal.