Effective mass theory for the anisotropic exciton in 2D crystals: Application to phosphorene
Elsa Prada, J. V. Alvarez, K. L. Narasimha-Acharya, F. J. Bailen, J. J. Palacios
EEffective mass theory for the anisotropic exciton in 2D crystals: Application tophosphorene
Elsa Prada,
1, 2, 3
J. V. Alvarez,
1, 2, 3
K. L. Narasimha-Acharya, F. J. Bailen, and J. J. Palacios
1, 2, 3 Departamento de F´ısica de la Materia Condensada,Universidad Aut´onoma de Madrid, Cantoblanco, 28049 Madrid, Spain Instituto de Ciencia de Materiales Nicol´as Cabrera,Universidad Aut´onoma de Madrid, Cantoblanco, 28049 Madrid, Spain Condensed Matter Physics Center (IFIMAC),Universidad Aut´onoma de Madrid, Cantoblanco, 28049 Madrid, Spain (Dated: March 15, 2016)We present a theoretical study of the exciton binding energy for anisotropic two-dimensionalcrystals. We obtain analytical expressions from variational wave functions in different limits of thescreening length to exciton size ratio and compare them with numerical solutions, both variationaland exact. As an example, we apply these results to phosphorene, a monolayer of black phospho-rous. Aided by density functional theory calculations for the evaluation of the two-dimensionalpolarizability, our analytical solution for the exciton binding energy gives a result which compareswell with numerical ones and, in turn, with experimental values, as recently reported.
I. INTRODUCTION
Since the mechanical exfoliation of graphene , researchin understanding the properties of two-dimensional (2D)crystals has increased in many folds. Atomically thinsingle-layered materials obtained from transition metaldichalcogenides , boron nitride , bismuth , etc. are be-ing extensively studied for applications as electronic andphotoelectronic devices. Few-layered black phosphorous(BP) is a recent addition to the list of graphene-inspiredmaterials . Apart from having a sizeable band gapwhich can be tuned by the manipulation of the number oflayers, the atomic structure of BP is highly anisotropicwhich leads to high asymmetry of the electronic bandstructure even for few layers. In particular, a single layerof BP or phosphorene is attracting most of the attention.The peculiar anisotropic nature of the band gap distin-guishes phosphorene from other 2D crystals, increasingits potential functionality.Excitons are a bound state of an electron and an holeand play an important role in the optical properties ofthe material. Understanding the nature of excitons andtheir dependence on the electronic structure of the hostmaterial is critical and lends a deeper perspective intothe many-body physics involved in 2D crystals. The 2Dnature of the polarizability of these crystals introducesan important length scale (screening length) r . For dis-tances between charges in the crystal plane, r , greaterthan r the electron-hole binding potential behaves likein a 3D system i.e., it goes as ∼ /r . However, for thecase where r is less than r the potential is 2D-like, i.e.,logarithmic. This behavior makes excitons in 2D crystalsdifferent from their 3D counterparts .Since most common 2D crystals are isotropic, the effectof anisotropy on the optical properties of these materialshas remained essentially unexplored. The appearance ofphosphorene has, however, changed this view and recentworks address this issue from an analytical , numerical , and first-principles standpoints. Here we give a detailedaccount of a variational approach, introduced by us inRef. 8, to the calculation of the exciton binding energy E X in anisotropic 2D crystals. Several analytical expres-sions are derived in certain limits of the 2D interactionpotential. The accuracy of our analytical expressions for E X is tested against both variational and exact numeri-cal solutions to the actual 2D potential, finding excellentagreement in a wide and experimentally relevant rangeof screening lengths. In particular, the value of E X forphosphorene, as obtained from our analytical expression,compares almost exactly to the numerical results. Fur-thermore, this value nicely agrees with the recently re-ported experimental result .The present work is divided as follows. In Sec. II wereview the form and limiting behaviour of the Coulombinteraction potential for charged particles in 2D systems.In Sec. III we present our variational approach based onan anisotropic exciton wavefunction. We first present theanalytical result for E X in the limiting case where the 2Dpotential reduces to the standard 3D Coulomb potential ∼ /r for isotropic 2D systems to later introduce theanisotropy and re-derive the binding energy for this case.In the same manner we derive analytical expressions forthe isotropic and anisotropic binding energies in the op-posite limit where the 2D interaction potential behaveslogarithmically. We also compare our analytical expres-sions with the numerically solved variational problem aswell as with the exact numerical solution. In Sec. IVwe propose an alternative variational approach based ongaussian orbitals. In Sec. V, after computing the 2Dpolarizability with density functional theory (DFT), ouranalytical approach is applied to the case phosphorene.Finally we present our conclusions in Sec. VI. a r X i v : . [ c ond - m a t . m e s - h a ll ] M a r II. BINDING PARTICLE-HOLE POTENTIALSIN 2D
As originally derived by Keldysh , the Coulomb po-tential energy created by a point charge at the origin thatelectrons feel in 2D layers follows the expression: V ( r ) = − e (cid:15) ¯ (cid:15)r (cid:20) H (cid:18) rr (cid:19) − Y (cid:18) rr (cid:19)(cid:21) , (1)where r ≡ d(cid:15)/ ( (cid:15) + (cid:15) ) and ¯ (cid:15) = ( (cid:15) + (cid:15) ) /
2. Here d isthe thickness of the 2D material, (cid:15) is its bulk dielectricconstant, and (cid:15) and (cid:15) are the dielectric constants ofthe surrounding media, typically substrate and vacuum.Here r plays the role of a screening length and sets theboundary between two different behaviours of the poten-tial. For r < r the potential diverges logarithmically,as if created by line charges. In this limit, the potentialtakes the simplified form also given by Keldysh : V ( r (cid:28) r ) ≈ e π(cid:15) ¯ (cid:15) r (cid:20) ln (cid:18) r r (cid:19) + γ (cid:21) , (2)where γ is the Euler constant. For r > r the poten-tial becomes the standard Coulomb potential created bypoint charges which decays as 1 /r : V ( r (cid:29) r ) ≈ − e π(cid:15) ¯ (cid:15) r . (3)A very good approximation to the Keldysh potential,fairly accurate in both limits and simpler to use, wasintroduced by Cudazzo et al. : V C2D ( r ) = e π(cid:15) ¯ (cid:15) r (cid:20) ln (cid:18) rr + r (cid:19) + ( γ − ln(2)) e − r/r (cid:21) . (4)It is interesting to compare these four expressions asa function of the distance r in a range of several ordersof magnitude both above and below r . We present sucha comparison in Fig. 1. There it can seen the range ofvalidity of each approximation, the Cudazzo et al. ex-pression being remarkably accurate for all distances. III. VARIATIONAL WAVEFUNCTIONAPPROACH
For generic 2D crystals with electrons and holes pre-senting anisotropic effective masses, m e ( h ) x (cid:54) = m e ( h ) y , weconsider variational solutions for the exciton wave func-tion of the type : φ ( x, y ) = (cid:18) a x λπ (cid:19) / exp (cid:16) − (cid:112) ( x/a x ) + ( y/λa x ) (cid:17) , (5)where λ is the variational anisotropy scaling factor relat-ing the exciton extension along the x direction, a x (whichis also a variational parameter) and the one along the y FIG. 1: (color online). Keldysh, 3D Coulomb,logarithmic, and Cudazzo et al. potentials in log-logscale in a range of distances spanning several decadesaround r .direction ( a y = λa x ). With this variational wavefunc-tion we can evaluate the expectation value of the kineticenergy: E kin ( a x , λ ) = (cid:126) (cid:90) (cid:90) φ (cid:20) µ x ∂ φ∂x + 1 µ y ∂ φ∂y (cid:21) dxdy = (cid:126) a x (cid:18) µ x + 1 λ µ y (cid:19) 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 E pot ( a x , λ ) = (cid:90) (cid:90) V ( x, y ) φ ( x, y ) dxdy (6)and the variational exciton binding energy is obtainedfrom the addition of these two quantities, E X ( a x , λ ) = E kin + E pot . (7)Upon minimization with respect to a x and λ , one obtainsthe optimal parameters defining the extension and shapeof the exciton and the actual binding energy E X . Resultsfrom three minimization procedures, one analytical andtwo numerical, are presented in next section. A. Analytical Results
The integral for the potential energy in Eq. (6) turnsout to be too difficult for an exact variational analyticalsolution. The main goal of this section is to make useof the asymptotic behaviour of the Keldysh potential toget analytical expressions for E X in the limits r (cid:29) r and r (cid:28) r , namely, valid for large and small excitons,respectively. r (cid:29) r limit We begin by evaluating E X in the isotropic case ( λ = 1, a x = a y = a ), considering only the long-range behaviourof the Keldysh potential (see Eq. 3). The contributionof the potential energy to E X is given in this limit by E pot = − e π(cid:15) ¯ (cid:15) a . (8)Now minimizing E X ( a ) with respect to the variationalexciton radius one obtains E X = − e π(cid:15) ¯ (cid:15) a , (9)where the minimal exciton radius ˜ a is given by˜ a = a ¯ (cid:15)m µ , (10)and m and a = π(cid:15) e (cid:126) m are the free electron mass andthe Bohr radius, respectively.For the anisotropic case ( λ (cid:54) = 1) the exciton extensionalong the x -direction is now given by˜ a x ( λ ) = a ¯ (cid:15)m (cid:18) µ x + 1 λ µ y (cid:19) I ( λ ) . (11)In the previous expression we find a function of λ definedthrough the elliptic integral I ( λ ) ≡ π (cid:90) π dθ (cid:112) λ −
1) cos θ = 1 π (cid:18) K (1 − /λ ) λ + K (1 − λ ) (cid:19) , (12)where the function K is the complete elliptic integral ofthe first kind. Defining now µ xy as µ xy ( λ ) ≡ (cid:18) µ x + 1 λ µ y (cid:19) − I ( λ ) , (13)the exciton extension along the x axis can now be writtenas ˜ a x ( λ ) = a (cid:15)mµ xy ( λ ) . (14)The exciton extension along the y direction is thus˜ a y ( λ ) = a ¯ (cid:15)m (cid:18) µ x + 1 λ µ y (cid:19) λI ( λ ) = a (cid:15)m λµ xy ( λ ) (15)and the λ -dependent binding energy of the exciton nowbecomes E X ( λ ) = − e π(cid:15) ¯ (cid:15) I ( λ )˜ a x ( λ ) . (16) We now define I E ( λ ) ≡ ( λ − dI ( λ ) /dλ + λI ( λ )= 1 π (cid:18) E(1 − /λ ) + E(1 − λ ) λ (cid:19) , (17)where E is the complete elliptic integral of second kind.We can see that the minimal λ , ˜ λ , satisfies in general thefollowing equation: µ x µ y = λ I E ( λ ) − λI ( λ ) I ( λ ) − λI E ( λ ) , (18)which has no analytical solution for λ . However, it canbe shown that for λ (cid:46) λ ≈ (cid:18) µ x µ y (cid:19) / . (19)Finally, notice that the results obtained in this subsectionwill be valid as long as the exciton extension in both x and y directions is much larger than r . The consistencyof this approximation for given experimental parameters( r , µ x , µ y , and ¯ (cid:15) ) has to be checked a posteriori. r (cid:28) r limit As r → a = (cid:114) ¯ (cid:15)mµ a r (20)and the binding energy of the exciton is E X = e π(cid:15) ¯ (cid:15) r (cid:20)
32 + ln (cid:18) ˜ a r (cid:19)(cid:21) . (21)For an anisotropic system the λ -dependent exciton ex-tension along the x direction is given by˜ a x ( λ ) = (cid:115) a r ¯ (cid:15)m (cid:18) µ x + 1 λ µ y (cid:19) . (22)Using now a different definition for µ xy µ xy ( λ ) ≡ (cid:18) µ x + 1 λ µ y (cid:19) − , (23)the exciton x -extension now becomes˜ a x ( λ ) = (cid:115) a r ¯ (cid:15)mµ xy ( λ ) . (24)Again, taking into account that a y = λa x , the λ -dependent minimal exciton extension along the y direc-tion is ˜ a y ( λ ) = (cid:115) a r ¯ (cid:15)mλ µ xy ( λ ) . (25)Note that ˜ a x ( µ x , µ y ) = ˜ a y ( µ y , µ x ).Finally we obtain the exciton energy for this case: E X ( λ ) = e π(cid:15) ¯ (cid:15) r (cid:20)
32 + ln (cid:18) ˜ a x ( λ )4 r λ + 12 (cid:19)(cid:21) , (26)where the minimal λ is˜ λ = (cid:18) µ x µ y (cid:19) / (27)for all µ x and µ y . Again, notice that this result will bevalid as long as the x and y minimal extensions of theexcitonic wave function are small compared to r .Note that Eq. (26) can be written in a more symmet-rical way as a function of both ˜ a x and ˜ a y , E X ( λ ) = e π(cid:15) ¯ (cid:15) r (cid:20)
32 + ln (cid:18) ˜ a x ( λ ) + ˜ a y ( λ )8 r (cid:19)(cid:21) , and that it is also symmetrical under exchange of µ x and µ y , as it should be. However, we find that the bindingenergy is not only a function of µ x /µ y , but depends onboth their values.Finally, for completeness, we present an analytical ex-pression for the exciton binding energy using the Cudazzopotential in the isotropic case: E X = e π(cid:15) ¯ (cid:15) (cid:20) a µ ˜ a + 4( γ − ln(2)) r (˜ a + 2 r ) − r (cid:18) γ + ln (cid:18) r ˜ a (cid:19)(cid:19) + ˜ a − r ˜ ar e r / ˜ a Ei (cid:18) − r ˜ a (cid:19)(cid:21) , (28)where Ei is the exponential integral function. We haveonly been able to obtain a working analytical expressionfor ˜ a (too cumbersome to be shown here) in the limit r (cid:38) a where the above expression is actually useful. B. Numerical Optimization and exact solution
To validate and test the accuracy of the limiting ana-lytical expressions given in the previous section, we nowuse the wavefunction in Eq. (5) to numerically com-pute the potential energy given by Eq. (6) for the exactKeldysh potential. We also solve, numerically as well, the2D Schr¨odinger equation for the same potential, whichwill give us the exact value of E X (down to the requirednumerical precision). Exciton binding energies for theisotropic case ( λ = 1) are presented in Fig. 2 as a func-tion of the screening length r . For comparison’s sake,we take ¯ (cid:15) = 1, i.e., the 2D crystal is suspended in vac-uum, and µ x = µ y = m/
2. Thus, according to Eq. (9), E x = − . r = 0. The numerical variational result com-pares very well with the exact numerical value in the large FIG. 2: (color online). Binding energies (in Ryd) forthe isotropic exciton computed in different ways:Numerical optimization of the variational wavefunctionin Eq. 5 (solid line), numerical solution of theSchr¨odinger equation (circles), analytical variationalexpression valid for large values of r (Eq. 21) (dashedline), and expression in Eq. 28 obtained for theCudazzo et al. potential (dotted-dashed line)range of explored screening lengths. For r (cid:29) a the an-alytical solution in Eq. (21) works fairly well. There,the size of the exciton is smaller than r and the 1 /r contribution to the Keldysh potential is negligible. Asexpected, the analytical solution starts to fail as r → a since there the size of the exciton becomes comparableto r and the long-range 1 /r contribution to the Keldyshpotential becomes dominant. (One should keep in mindthat the limit of validity of the analytical result, as shownin Eq. (20), depends on the values of ¯ (cid:15) and µ .) We alsocompare with the result given by Eq. (28), obtained us-ing the approximate expression to the potential in Eq.(4). This expression, although not as friendly as the pre-vious one, extends the limit of validity of our analyticalresults down to r ≈ a .The results for the anisotropic case are presented inFig. 3 for r = 20 a and r = 400 a as a function of theanisotropy ratio µ y µ x with µ x = m/ µ x µ y with µ y = m/
2. Once again there is closeagreement between the analytical solution [Eq. (26)], thenumerical optimization, and the exact numerical solutionfor large r , while for the smaller value, the analyticalsolution visibly deviates from the other two.An important prediction of our analytical results is therelation between the anisotropy in the exciton extensionand the effective masses: ˜ λ = (cid:16) a y a x (cid:17) ∼ (cid:16) µ x µ y (cid:17) / , whichbecomes exact in the limit of small excitons. To test thisrelation we fitted the optimal value of the variational (a)(b) FIG. 3: (color online). Exciton binding energy as afunction of the anisotropy µ y /µ x (for fixed µ x = m/ r = 20 a (a) and r = 400 a (b).parameter ˜ λ to the law˜ λ = C ( r ) (cid:18) µ x µ y (cid:19) α ( r ) (29)for a large range of r . The results of this fit are presentedin Fig. 4. They confirm our analytical results and recoverthe exact 1/3 exponent in the limit of large r .We finally provide a comparison between the exactand variational wave functions for several values of r in the isotropic limit (see Fig. 5). Note that the distanceis rescaled with the optimal radius and the amplitudeof the wave function with the normalization constant A = (cid:16) a x π (cid:17) / .This representation illustrates to whatextent the exact and variational wave functions satisfy FIG. 4: (Color online). Parameters of the fit inexpression (29) as a function of r : C ( r ) (left axis) and α ( r ) (right axis).FIG. 5: (Color online). Exact (dashed lines) andvariational (solid line) wave functions for r = 0 . , , , variational extension of the exciton and the amplitudewith the normalization constant of the variational wavefunction A = (cid:16) a x π (cid:17) / .similar scaling relations. Note that at r = 0 the exactwave functions do not show the prominent cusp of a 1sSlater-type orbital. This softened behavior at the originsuggests than a combination of gaussian functions maycapture more accurately this feature of the wave function,as shown in next section. IV. GAUSSIAN-BASIS VARIATIONALMETHOD
We have found that in the limit of very small r , thebinding energy is very sensitive to small changes in r .Furthermore, the numerical solution of the Schr¨odingerequation requires a very fine mesh to reproduce thebound state in such a limit. On the other hand, theanalytical result found in the 1 /r limit of the potentialconstitutes an isolated point an thus cannot be easilyextended to small but finite values of r . It is there-fore interesting to find an alternative numerical methodto study anisotropic excitons in the limit of small r .Moreover, as we have presented in Fig. 5, the behaviorof the exciton exact wave functions for different valuesof the screening length r resembles more a 1s Gaussianthan a Slater-type orbital. Gaussian-type orbitals (GTO)are very efficient basis sets used intensively in quantumchemistry and solid state calculations.The gaussian basis functions { χ p } follow the expres-sion: χ p = e − ( α px x + α py y ) . (30)The index p is an integer that here has been chosen torun from 1 to 4 and the exponents α are coefficients thathave to be optimized to minimize the ground state energyobtained by the variational method. Unlike the conven-tional gaussian approach, the anisotropy of the problemintroduces two coefficients α px , α py per GTO. We limit thevariational freedom assuming that the anisotropy is iden-tical for the four basis functions: α py = κα px . (31)In this equation, κ is a constant that does not dependon p so that we reduce the number of exponents that wehave to optimize from eight to four. The variational wavefunction in the GTO basis is given by φ G ( x, y ) = (cid:88) p =1 C p χ p . (32)For fixed values of α p , the energy is computed by gen-eralized diagnonalization of a 4x4 matrix. The matrixelements of the kinetic energy are H kin pq = πM pq (cid:18) µ x α px α qx α px + α qx + 1 µ y α py α qy α py + α qy (cid:19) . (33)The matrix elements of the potential energy are com-puted by numerical integration, H pot pq = (cid:90) (cid:90) V K ( x, y ) χ p χ q dxdy, (34)and the overlap matrix is S pq = πM pq , (35) FIG. 6: (Color online). Ground state energies obtainedusing the gaussian variational method (continuous lines)and the numerical solution of Schr¨odinger equation for r = 10 a (circles) as a function of the effective massratio µ y /µ x for µ x = m .where M pq = (cid:112) ( α px + α qx )( α py + α qy ).The binding energy and the optimal wave function areobtained by minimizing numerically the energy with re-spect to the five variational parameters. Efficient opti-mization of the energy requires a careful choice of the ini-tial guess for the values of the exponents α p . In our case,we choose the optimal values for a 1s orbital of a “2D hy-drogen atom”, taking µ x = µ y = m and r approachingzero. Once these optimal exponents are obtained, r ischanged slightly and the problem is solved again, usingthis time the exponents α obtained in the previous step.The procedure continues until the ground state energyfor the desired r is reached.For example, fixing r = 10 a and µ x = m , the effec-tive mass along the y axis µ y is varied from 1 to 40 m .In this case the initial guess for the exponents is the lastset of coefficients α obtained when changing r . Theresult for the ground state energy of this calculation isshown in Fig. 6 in comparison to the numerical solutionof the Schr¨odinger equation with the Keldysh potential.They match perfectly. Two other values of r /a ob-tained with the same gaussian-basis variational methodare also shown.In Fig. 7, the length scales a px = 1 / √ α px of the wholeset of optimal GTO’s for r = 10 a are plotted vs µ y /µ x for µ x = m . κ is also plotted in Fig. 8 against the samequantity. In Fig. 9 we show a log-log representation of κ versus the asymmetry ratio for different values of r where a linear fitting has been made. V. APPLICATION TO PHOSPHORENE
A paramount example of an anisotropic 2D crystal isphosphorene , where the effective masses along x andFIG. 7: (Color online). Coefficients a px = 1 / √ α px for r = 10 a obtained versus µ y /µ x for µ x = m .FIG. 8: (Color online). Coefficient κ = α py /α px for r = 10 a obtained versus µ y /µ x for µ x = m .FIG. 9: (Color online). Log-log representation of thecoefficient κ = α py /α px obtained for different values of r (circles) and its linear fitting (continuous line) versus µ y /µ x for µ x = m . y directions can differ by even an order of magnitude.From the start we have chosen to express the 2D potentialconstant in terms of the bulk dielectric constant (cid:15) andthe effective thickness of the 2D crystal d (see Eq. 1).Equivalent expressions for the 2D potential, which relyon the evaluation of the actual 2D polarizability of the2D crystal, χ , have recently been proposed . These areprobably more appropriate for actual crystals, althoughit has also been shown that Eq. (1) works well as long as (cid:15) is taken as the in-plane component of the bulk dielectrictensor of the 3D crystal. Here we will compare bothpossibilities.As shown in Ref. 9, the screening length r dependson the polarizability χ as r ≡ πχ . The polarizabilityfor 2D materials can be computed using the expression (cid:15) ( L ) = 1 + 4 πχL , (36)where L is the distance between layers in a 3D lay-ered structure. As can be seen, the dielectric func-tion (cid:15) tends to unity as the inter layer distance L tends to infinity. We have computed the dielectricfunction at different inter layer distances within thedensity functional theory framework using the Perdew-Burke-Ernzerhof (PBE) functional and norm conserv-ing Troullier-Martins (TM) pseudopotentials as availablein the SIESTA package . The atomic and electronicstructures have been duly converged on all parameters.SIESTA calculates the imaginary part of the dielectricfunction from which the real part of it is obtained usingthe Kramers-Kronig relations. In order to account for theunder-estimated band gap, the scissor approximation, asimplemented in SIESTA, has been utilized. The scissorshift of 1.2505 eV was made to match our previously re-ported band gap value of 2.15 eV . While more elegantapproaches to the gap problem of phosphorene have beenreported in the literature , the scissor approximation suf-fices to our purpose here.Using the plane-averaged static dielectric function cal-culated with SIESTA (see Fig. 10), a value in the vicinityof χ of 3.8 ˚A is obtained. This value for χ yields a screen-ing length of r = 23 . a x = 11 . a y = 5 . r , the use of the an-alytical expressions obtained in the logarithmic limit ofthe potential is justified. Also this was expected fromFig. 2 and, in particular, from the comparison shown inFig. 3 for anisotropic cases. There it can be seen that al-ready for r = 20 a the deviation between the analyticalresult and the numerical ones is less than 10% for a ratio µ y /µ x ≈ E X = 0 .
61 eV, while thenumerical variational value is E X = 0 .
78 eV. This resultis remarkably close to a recently reported experimentalvalue of ≈ . x and y components of the dielectric constant (evaluated at zerofrequency) for different values of the interlayer distance L in a 3D black phosphorous structure. (b)Polarizability as obtained from Eq. 36 for differentvalues of L after averaging on the plane.substrate . A more recent experiment, however, reportsa smaller value for E X , which is maybe more expecteddue to the screening of the substrate .Similarly, we can use the value of r obtained from thereal part of the bulk (cid:15) (at zero frequency) and the thick-ness d of the monolayer. Since the thickness is somewhatundetermined, we have computed the binding energy forvarying d , as shown in Fig. 11. It can be observed thatthe binding energy of the monolayer computed using themicroscopically derived χ matches the binding energy ob-tained using r ≡ d(cid:15)/ ( (cid:15) + (cid:15) ) for d ≈ VI. CONCLUSIONS
We have shown that a variational approach to the exci-ton binding energy in anisotropic 2D crystals can give ex-cellent results when compared to numerical approaches.FIG. 11: (Color online). Exciton binding energy (in eV)as a function of the thickness when all the otherparameters (effective mass, bulk dielectric constant, andscreening length) to correspond to phosphorene. Thehorizontal line marks the value obtained using theactual 2D polarizability of phosphorene.Furthermore, we have studied the range of validity of an-alytical solutions to the variational approach and foundthat these can give highly satisfactory results in a rangeof values of screening lengths which is relevant for actual2D crystals such as phosphorene. We have computedthe exciton binding energy in this case and found a verygood agreement with a recently reported experimentalresult . As long as the screening length to exciton sizeis large, our analytical results can be trivially used topredict the exciton binding energy of any 2D crystal. ACKNOWLEDGMENTS
This work was supported by MINECO under GrantsNos. FIS2013-47328 and FIS2012-37549, by CAM un-der Grants Nos. S2013/MIT-3007, P2013/MIT-2850,and by Generalitat Valenciana under Grant PROME-TEO/2012/011. The authors thankfully acknowledge thecomputer resources, technical expertise, and assistanceprovided by the Centro de Computaci´on Cient´ıfica of theUniversidad Aut´onoma de Madrid. E.P. also acknowl-edges the Ram´on y Cajal Program. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science , 666 (2004). J. N. Coleman, M. Lotya, A. O’Neill, S. D. Bergin, P. J.King, U. Khan, K. Young, A. Gaucher, S. De, R. J. Smith,I. V. Shvets, S. K. Arora, G. Stanton, H.-Y. Kim, K. Lee,G. T. Kim, G. S. Duesberg, T. Hallam, J. J. Boland,J. J. Wang, J. F. Donegan, J. C. Grunlan, G. Moriarty,A. Shmeliov, R. J. Nicholls, J. M. Perkins, E. M. Grieveson,K. Theuwissen, D. W. McComb, P. D. Nellist, and V. Ni-colosi, Science, Volume 331, Issue 6017, pp. 568- (2011). , 568 (2011). M. Kawaguchi, S. Kuroda, and Y. Muramatsu, Journal ofPhysics and Chemistry of Solids , 1171 (2008), 14th In-ternational Symposium on Intercalation Compounds ISIC14. C. Sabater, D. Gos´albez-Mart´ınez, J. Fern´andez-Rossier,J. G. Rodrigo, C. Untiedt, and J. J. Palacios, Phys. Rev.Lett. , 176802 (2013). L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H.Chen, and Y. Zhang, Nature Nanotechnology, Volume 9,Issue 5, pp. 372-377 (2014). , 372 (2014). H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tom´anek,and P. D. Ye, ACS Nano , 4033 (2014). V. Tran, R. Soklaski, Y. Liang, and L. Yang, Physical Re-view B, Volume 89, Issue 23, id.235319 , 235319 (2014). A. Castellanos-Gomez, L. Vicarelli, E. Prada, J. O. Island,K. L. Narasimha-Acharya, S. I. Blanter, D. J. Groenendijk,M. Buscema, G. A. Steele, J. V. Alvarez, H. W. Zandber-gen, J. J. Palacios, and H. S. J. v. d. Zant, 2D Materials , 025001 (2014). P. Cudazzo, I. V. Tokatly, and A. Rubio, Phys. Rev. B , 085406 (2011). T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman,Phys. Rev. B , 045318 (2013). A. S. Rodin, A. Carvalho, and A. H. Castro Neto, PhysicalReview B, Volume 90, Issue 7, id.075429 , 075429 (2014). X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia,H. Zhao, H. Wang, L. Yang, X. Xu, and F. Xia, Nat.Nanotechnology , 517 (2015). L. V. Keldysh, JETP Lett. , 658 (1979). A. Schindlmayr, European Journal of Physics, Volume 18,Issue 5, pp. 374-376 (1997). , 374 (1997). Y. Zhang and W. Yang, Physical Review Letters, Volume80, Issue 4, January 26, 1998, p.890 , 890 (1998). J. M. Soler, E. Artacho, J. D. Gale, A. Garc´ıa, J. Junquera,P. Ordej´on, and D. S´anchez-Portal, Journal of Physics:Condensed Matter , 2745 (2002).17