Improvements in the GW and BSE calculations on phosphorene

# Improvements in the $GW$ and Bethe-Salpeter-equation calculations on phosphorene

## Abstract

Phosphorene is a bidimensional material that has properties useful for semiconductor devices. In this work we studied the electronic and optical properties of this material using the approximation and the Bethe-Salpeter equation (BSE) methods. We stress the importance of a careful convergence study of the most relevant parameters, and we show how they affect the result of the calculations. A comparison with previous results is given. The QP band gap obtained was 2.06 eV and it is in good agreement with experimental results. BSE calculations were performed on top of in order to include excitonic effects. The absorption spectrum was analyzed and an optical gap of 1.22 eV was obtained. The calculated excitonic binding energy is 0.84 eV, also in good agreement with experimental results.

## I Introduction

Since its recent synthesisLiang et al. (2014), phosphorene has been attracting interest due to its peculiar and distinct properties. Unlike graphene, this material is a semiconductor with a band gap of approximately 2.0 eV,Liang et al. (2014) which is much higher than the bulk band gap of 0.3 eV.Keyes (1953) A wide band gap makes it very useful for semiconductor applications. Studies have reported high mobilityQiao et al. (2014) and high current on-off ratio for field effect transistors.Liu et al. (2014) Phosphorene can be stacked in two or more layers due to the weak Van der Walls interlayer interactions.Li et al. (2014) It was shown that the band gap value can be tuned by the number of stacked layers.Tran et al. (2014) This stacking can be used to control the electronic and optical properties by changing the number of the stacked layers. The infrared properties of Phosphorene have been recently measuredZhang et al. (2017) and found to depend on uniaxial strain and useful for infrared photonics. Band alignment on heterostructures involving phosphorene has been theoretically studiedÖzçelik et al. (2016); Ganesan et al. (2016) and a high power conversion efficiency for solar cells was found.Ganesan et al. (2016) Another interesting property is the strong in-plane anisotropy of transport and optical properties. This is due to its anisotropic crystal structure, shown in Figure 1. Phosphorene has a puckered lattice due to sp hybridization, which also allows band gap engineering by applying strain in-plane or out-plane.Çak ır et al. (2014) It was also shown that the deformation of phosphorene in a specific direction can be used to induce a semiconductor-metal transition.Fei and Yang (2014) Phosphorene was shown to have negative Poisson’s ratioJiang and Park (2014) and mechanical flexibilityWei and Peng (2014) that allows its use under extreme mechanical conditions. Also it can be used for gas sensorsKou et al. (2014) and solar cell applications.Deng et al. (2014)

Density Functional Theory (DFT) has been extensively used to calculate electronic and optical properties. Yet, it is well known that DFT, in the most common exchange-correlation approximations used, does not provide a good description of the electronic band structures, namely the band gap is underestimated. Also, because phosphorene has low dimensionality, the excitonic effects can not be neglected and must be taken into account if we want to study its optical properties. DFT does not include excitonic effects.

In this paper we study the electronic properties of phosphorene using the quasiparticle (QP) approximation.Hedin (1965); Hedin and S. (1969) To calculate the spectrum we have to determine the electron self-energy operator , which is equivalent to solve a many-body problem. The self-energy operator is approximated and expressed in terms of the screened Coulomb potential as , with being the single-particle Green’s function. To include the excitonic effects and obtain the optical properties we use a two-particle Green’s function formalism, through the Bethe-Salpeter equation (BSE)Salpeter and Bethe (1951); Albrecht et al. (1998) on top of the approximation.

Phosphorene has been studied before using the approximationLiang et al. (2014); Rudenko et al. (2015); Rudenko and Katsnelson (2014) and also the BSE method,Tran et al. (2014); Çak ır et al. (2014); Tran et al. (2015) but the results of these works do not agree (see tables 1 and 2). One possible reason may be different criteria for convergence. Ref. Çak ır et al., 2014 calculates with 112 and 326 empty bands and a vacuum size of 28.34 bohr, while Ref. Rudenko et al., 2015 and Rudenko and Katsnelson, 2014 used 360 bands (90 bands per atom). Ref. Tran et al., 2014 used 10 times the valence bands to converge the dielectric matrix. These numbers of bands are small when compared to the ones that we found good for convergence in our calculations. Ref. Tran et al., 2015 does not mention the number of bands used for the calculation. The truncation technique, which is needed to avoid interactions with periodic images, is not used in some of these works.Liang et al. (2014); Çak ır et al. (2014); Rudenko et al. (2015); Rudenko and Katsnelson (2014)

Concerning the BSE results, we found that we have to use a grid of -points in order to achieve convergence (see subsection III.3 below), a value larger than some of the works published.Çak ır et al. (2014); Tran et al. (2014, 2015)

In this work we study phosphorene properties with BSE, since to our knowledge very few works have been done. We also expect to contribute to clarify the convergence studies. The BSE approach is much more computationally expensive than DFT, because of the interdependence of the convergence parameters and higher sensitivity to the number of -points and number of bands, as will be explained in Sec. II. A common error may be to not consider that convergence parameters are interdependent and so can result in a false convergence behavior. For instance, Ref. Shih et al., 2010, shows that one of the causes of underestimation of the QP band gap of ZnO was the result of an inadequate treatment of the convergence parameters. The demanding resources that are needed to obtain rigorous results in BSE calculations may enhance this problem.

This paper is organized as follows. In Sec. II, the computational details are described, with an emphasis on the parameters that may affect convergence. Section III includes the results of the calculations done with the approximation and the BSE and the discussion. The conclusions are given in Sec. IV.

## Ii Method

This work was done in three steps. In the first step, DFT calculations were performed using the open source software package Quantum ESPRESSO.Giannozzi et al. (2009) In the second step, single-shot () calculations were done using the DFT results. Lastly, BSE was applied on calculations. Both the and BSE calculations were performed with the software package BerkeleyGW.Deslippe et al. (2012)

### ii.1 DFT calculations

DFT calculations were done using a scalar-relativistic norm-conserving pseudopotential. The exchange-correlation functional used was the generalized gradient approximation of Perdew-Burke-Ernzerhof (GGA-PBE).Perdew et al. (1996) The plane-wave cut-off used was 70 Ry and the integration over the Brillouin-zone was performed using the scheme proposed by Monkhorst-Pack Monkhorst and Pack (1976) with a grid of k-points. Both the plane-wave cut-off and the k-points grid were chosen after a convergence analysis. A vacuum size () between the layer images equal or greater than 20 bohr was enough to avoid interactions between the periodic images.

### ii.2 G0w0 calculations

calculations were performed on top of DFT. We used the Generalized-Plasmon-Pole (GPP) model proposed by Hybertson and Louie,Hybertsen and Louie (1986) where we compute the static dielectric matrix and then extend it to finite frequencies. This is done by considering the static polarizability matrix

 χGG′(q,0)=∑kocc∑nemp∑n′M∗nn′(k,q,G)Mnn′(k,q,G′)×1Enk+q−En′k, (1)

where

 Mnn′(k,q,G)=⟨nk+q|ei(q+G).r|n′k⟩ (2)

are the plane-wave matrix elements, is the momentum transfer, is the reciprocal lattice vector, and are the band numbers of occupied () and empty () bands respectively, and are the DFT eigenvectors and eigenvalues respectively. The plane-wave matrix elements are evaluated up to , where defines the dielectric matrix cut-off . The dielectric matrix is given by

 ϵGG′(q,0)=δGG′−v(q+G)χGG′(q,0), (3)

where is the bare Coulomb potential. To construct we have to do a summation on the empty bands and on the vectors which are given by the . These parameters have to be consistent with each other: if we choose for instance 100 empty bands for the summation, then the energy of the 100th empty band is the energy. So we have just one convergence parameter for the construction of the : either or the number of empty bands used. Here we will usually use the number of empty bands for this convergence criterion. The GPP model allows us to extend the static dielectric matrix to a finite frequency dependent matrix and then calculate . It is divided in two terms , where is the screened-exchange term and the is the Coulomb-hole term and they have the following expression

 ⟨nk|ΣSX(ω)|n′k⟩=−occ∑n′′∑qGG′M∗n′′n(k,−q,−G)Mn′′n′(k,−q,−G′)×[δGG′+Ω2GG′(q)(1−i tanϕGG′(q))(ω−En′′k−q)2−~ω2GG′(q)]v(q+G′), (4)
 ⟨nk|ΣCH(ω)|n′k⟩=12∑n′′∑qGG′M∗n′′n(k,−q,−G)Mn′′n′(k,−q,−G′)×Ω2GG′(q)(1−i tanϕGG′(q))~ωGG′(q)(ω−En′′k−q−~ωGG′(q))v(q+G′). (5)

The is the effective bare plasma frequency and it is given by

 Ω2GG′(q)=ω2p(q+G)⋅(q+G′)∣∣q+G∣∣2ρ(G−G′)ρ(0), (6)

where is the electron charge density in reciprocal space and is the classical plasma frequency defined by . The is the GPP mode frequency and is related to by the following equations

 ~ω2GG′(q)=|λGG′(q)|cosϕGG′(q), (7)
 |λGG′(q)|eiϕGG′(q)=Ω2GG′(q)δGG′−ϵ−1GG′(q,0), (8)

where and are the amplitude and phase of the renormalized respectively. For computing the self-energy operator matrices it is necessary to construct the plane-wave matrix just like in the case of the dielectric matrix. For the terms and we have to choose the screened-Coulomb cut-off. This cut-off has the same role of the dielectric cut-off, which is to define an energy truncation for the evaluation of the elements of . This energy has to be less or equal to . We chose it to be equal to in all calculations. The convergence of the number of bands included in the summation of eq. 5 has to be studied.

So we have to consider two convergence parameters: (or number of bands ) of and the number of bands included in .

### ii.3 BSE calculations

A detailed explanation of the BSE theory can be found in Ref. Rohlfing and Louie, 2000. For each excitonic state , BSE can be written

 (EQPck−EQPvk)ASvck+∑v′c′k′⟨vck|Keh|v′c′k′⟩=ΩSASvck, (9)

where and are the QP energies for the conduction () and valence bands () respectively, is the electron-hole kernel, and the excitonic wavefunction and energy respectively for an excitonic state . Tamm-Dancoff approximationDancoff (1950) is considered, where only valence to conduction band transitions are included. The first step to solve the BSE is to compute the electron-hole kernel. It can be separated in two terms , where is the screened direct interaction term and is the bare exchange interaction term. They are defined in the following way

 ⟨vck|Kd|v′c′k′⟩=∑GG′M∗c′c(k,q,G)WGG′(q,0)Mv′v(k,q,G′), (10)
 ⟨vck|Kx|v′c′k′⟩=∑G≠0M∗vc(k,q,G)v(q,0)Mv′c′(k,q,G) (11)

where the static screened Coulomb potential has the following expression

 WGG′(q,0)=ϵ−1GG′(q,0)v(q+G′). (12)

We also have to choose a cut-off for the kernel matrix construction because of matrices . The excitonic properties are very sensitive to the grid of -points that is used because the contributions of are very important. We then have to compute the kernel in very fine grids, which is computationally prohibitively expensive. The BerkeleyGW package computes the kernel with a coarse grid of -points and then interpolate it in a very fine grid of -points. To perform this interpolation we have to choose a grid finer than the coarse grid and include a number of valence and conduction bands. After solving the BSE, the imaginary part of the dielectric function, which is proportional to absorption spectrum, can be calculated by using the expression

 ϵBSE2(ω)=16π2e2ω2∑S|e⋅⟨0|v|S⟩|2δ(ω−ΩS). (13)

The term is called the velocity matrix element in the direction of the polarization of the light . We can also calculate the absorption spectrum without excitonic interactions with the expression

 ϵRPA2(ω)=16π2e2ω2∑vck|e⋅⟨vk|v|ck⟩|2δ(ω−EQPck+EQPvk), (14)

which is a random phase approximation (RPA) calculation on top of the results.

### ii.4 Coulomb truncation

To eliminate the artificial interactions between periodic images we can use the Coulomb truncation technique. The Wigner-Seitz slab truncation implemented in BerkeleyGW is used. The Coulomb interaction is truncated at the edges of the unit cell in the direction perpendicular to the slab plane. The bare Coulomb potential is then corrected to

 v(q)=4πq2(1−e−qxyzccos(qzzc)), (15)

where is the truncation distance in the direction perpendicular to the phosphorene plane and and are the parallel and perpendicular components of the vector respectively.

## Iii Results and discussion

### iii.1 DFT results

The DFT optimized lattice parameters for the phosphorene structure are bohr and bohr (see Fig. 1). The electronic band structure is shown in Fig. 2, and the band gap is direct at the point with an energy of 0.87 eV. This value is in good agreement with several works that use the GGA as the exchange-correlation functional.Wang et al. (2015b); Qiao et al. (2014); Tran et al. (2014) When HSEHeyd et al. (2003) functionals are used, the band gap increases to about 1.5 eV.Qiao et al. (2014); Wang et al. (2015b) To our knowledge there are two experimental works where the band gap of phosphorene was measured, with values of 2.05 eV (Ref. Liang et al., 2014) and eV (Ref. Wang et al., 2015a). As expected, these values are much larger than the DFT results, and confirm the need of using a QP theory to have better prediction of the electronic properties of phosphorene.

### iii.2 G0w0 results

We calculated the QP band structure of phosphorene using the approach. Fig. 2 shows the band structure of phosphorene with and without the QP corrections. The QP band structure is essentially a shift of the DFT conduction bands and the shape of the band structure changes little with the QP corrections.

First we studied the convergence of the QP band gap with the number of bands and , using a grid of -points and a vacuum size of 40 bohr. Convergence was achieved for bands for (which corresponds to an of 9.73 Ry) and bands for . The band gap calculated is 2.18 eV, much larger than the DFT values and closer to the experimental.

Figure 3 shows that the two convergence parameters are mutually dependent. If we choose to construct with 300 bands, a convergence of 0.01 eV is achieved by including only about 400 bands in the summation, and we would have a 1.95 eV band gap. But if we choose bands to construct the , we will need bands to achieve the same convergence criterion, and would have a band gap of 2.18 eV. These convergence parameters are interdependent because few bands in the calculation means including few vectors, and so we are preventing the contribution from high energy empty bands to the calculation of .Shih et al. (2010)

calculations also depend on the number of or -points. Despite convergence at DFT level being achieved for a small -points sampling, for it usually can only be achieved for a much denser grid. This is due to the large variation of the dielectric function near the point for bidimensional materials.Rasmussen et al. (2016) In order to capture this variation, a denser grid of -points is needed. We repeated the previous calculations with a denser grid of -points. Figure 4 shows the convergence of the QP band gap with that grid.

Convergence was achieved for a QP band gap of 2.11 eV, which is only 0.07 eV different from the calculations with the smaller grid of -points. This shows that for phosphorene calculations are not so dependent on the density of -points. To achieve a convergence of 0.01 eV we used even denser -point grids. The results are shown in Table 3 alongside with results for different vacuum sizes.

It is well known that the QP band gap has a strong dependence on the vacuum size at the level because of the non-local screening effects of the approximation which makes the gap converge as .Inkson (1971) A possible explanation for the increase of band gap with is that the screening becomes weaker due to the electron-electron correlation being stronger. The weaker the screening the larger the band gap. In order to reduce this dependence, we use the truncation method mentioned in Subsection II.4. Table 3 shows that the QP band gap indeed increases with vacuum size, even with truncation. We also made these calculations without truncation and we achieved the much lower value for the band gap of 1.81 eV and a stronger vacuum dependence. This shows that using truncation method is indeed necessary to obtain the right results. We also notice that when we increase the vacuum size, more -points are needed to achieve the convergence criterion of 0.01 eV. The larger the vacuum size the larger the variation of the dielectric function near the point, and so the denser the point sampling needs to be.Hüser et al. (2013)

In Table 1 we summarize the results from other works. The QP band gap we obtained is in excellent agreement with the experimental ones.

### iii.3 BSE results

The kernel matrix of the BSE calculations used a coarse grid of -points. The starting point for our BSE calculations was chosen to be the results with a vacuum size of 40 bohr. For this calculation, we verified the convergence of the bands above and below the band gap, for the bands that were used for the absorption spectrum calculation. We found that most of the energies were converged to better than 0.01 eV, with some cases with higher values, up to 0.014 eV, for the bands farther away from the gap. We deemed this good enough, specially since the most relevant bands were well converged. The first calculation was done with a fine grid for interpolation of -points. We tested for 4 and 5 valence and conduction bands in the kernel matrix and for interpolation, and we found no difference for energies below 4 eV, as can be seen in Fig. 5. So we only included 4 valence bands and 4 conduction bands for the construction of the kernel matrix and for interpolation.

We also tested the convergence with respect to the number of -points used for the interpolation. Fig. 6 shows that both the BSE and the RPA absorption spectrum are converged for a grid of (8100) -points.

In Fig. 7 a comparison between the BSE and the RPA spectra is shown. When excitonic effects are included, the spectrum red-shifts and discrete peaks of the excitonic states appear, as expected. The first excitonic peak is at an energy of 1.22 eV. This is the optical gap of phosphorene. Knowing that the optical gap is 1.22 eV and the QP band gap is 2.06 eV, the excitonic binding energy is 0.84 eV. This value is in agreement with the experimental one of eV (Ref. Wang et al., 2015a). A high excitonic binding energy is expected in this material because of its low dimensionality and low screening. The low dimensionality increases the confinement between the electron and hole which enhances its Coulomb interaction.

Table 2 summarizes the results of our work and three other works that we found that use BSE approach. It also shows the only experimental value we could find. The optical gap we obtained is a bit smaller than the experimental value. This may be due to substrate effects in the experiment; our calculation assumes a suspended monolayer while the experimental measurements were done on a substrate. Our results also differ a bit from other theoretical works. It is not clear to us if calculations were fully converged in other works as we discussed above. And our results showed that a grid of -points is needed in order to achieve convergence, much larger than the ones used by others.

## Iv Conclusion

approximation and the BSE methods were used to study the electronic and optical properties of phosphorene in a careful and systematic numerical convergence study. We found that the QP band gap is 2.06 eV, closer to the experimental values than the results from other works. We also did BSE calculations in order to include the excitonic effects which are important in materials with low dimensionality. An optical gap of 1.22 eV was obtained and an excitonic binding energy of 0.84 eV. These values are consistent with a recent experimental work and results from other theoretical works. We also saw that in order to achieve convergence using those techniques we have to consider the interdependence of the parameters that have to be converged.

## Acknowledgments

R.M. Ribeiro acknowledge support from the European Commission through the project “Graphene-Driven Revolutions in ICT and Beyond” (Ref. No. 696656), COMPETE2020, PORTUGAL2020, FEDER and the Portuguese Foundation for Science and Technology (FCT) through project PTDC/FIS-NAN/3668/2014 and in the framework of the Strategic Financing UID/FIS/04650/2013.

### References

1. L. Liang, J. Wang, W. Lin, B. G. Sumpter, V. Meunier,  and M. Pan, Nano Letters 14, 6400 (2014), pMID: 25343376.
2. R. W. Keyes, Phys. Rev. 92, 580 (1953).
3. J. Qiao, X. Kong, Z.-X. Hu, F. Yang,  and W. Ji, Nat Commun 5 (2014), 10.1038/ncomms5475.
4. H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek,  and P. D. Ye, ACS Nano 8, 4033 (2014), pMID: 24655084.
5. L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen,  and Y. Zhang, Nat Nano 9, 372 (2014), article.
6. V. Tran, R. Soklaski, Y. Liang,  and L. Yang, Phys. Rev. B 89, 235319 (2014).
7. G. Zhang, S. Huang, C. Chaves, A.and Song, V. Özçelik, T. Low,  and H. Yan, Nature Communications 8, 14071 (2017).
8. V. O. Özçelik, J. G. Azadani, C. Yang, S. J. Koester,  and T. Low, Phys. Rev. B 94, 035125 (2016).
9. V. D. S. Ganesan, J. Linghu, C. Zhang,  and Y. P. Feng, Applied Physics Letters 108, 122105 (2016).
10. D. Çak ır, H. Sahin,  and F. M. Peeters, Phys. Rev. B 90, 205421 (2014).
11. R. Fei and L. Yang, Nano Letters 14, 2884 (2014), pMID: 24779386.
12. J.-W. Jiang and H. S. Park, Nature Communications 5, 4727 EP (2014), article.
13. Q. Wei and X. Peng, Applied Physics Letters 104, 251915 (2014).
14. L. Kou, T. Frauenheim,  and C. Chen, The Journal of Physical Chemistry Letters 5, 2675 (2014), pMID: 26277962.
15. Y. Deng, Z. Luo, N. J. Conrad, H. Liu, Y. Gong, S. Najmaei, P. M. Ajayan, J. Lou, X. Xu,  and P. D. Ye, ACS Nano 8, 8292 (2014), pMID: 25019534.
16. L. Hedin, Phys. Rev. 139, A796 (1965).
17. L. Hedin and L. S., Solid State Physics 23 (1969).
18. E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 (1951).
19. S. Albrecht, L. Reining, R. Del Sole,  and G. Onida, Phys. Rev. Lett. 80, 4510 (1998).
20. A. N. Rudenko, S. Yuan,  and M. I. Katsnelson, Phys. Rev. B 92, 085419 (2015).
21. A. N. Rudenko and M. I. Katsnelson, Phys. Rev. B 89, 201408 (2014).
22. V. Tran, R. Fei,  and L. Yang, 2D Materials 2, 044014 (2015).
23. X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia, H. Zhao, H. Wang, L. Yang, X. Xu,  and F. Xia, Nat Nano 10, 517 (2015a).
24. B.-C. Shih, Y. Xue, P. Zhang, M. L. Cohen,  and S. G. Louie, Phys. Rev. Lett. 105, 146401 (2010).
25. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari,  and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009).
26. J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen,  and S. G. Louie, Computer Physics Communications 183, 1269 (2012).
27. J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
28. H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
29. M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
30. M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000).
31. S. M. Dancoff, Phys. Rev. 78, 382 (1950).
32. V. Wang, Y. Kawazoe,  and W. T. Geng, Phys. Rev. B 91, 045433 (2015b).
33. J. Heyd, G. E. Scuseria,  and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003).
34. F. A. Rasmussen, P. S. Schmidt, K. T. Winther,  and K. S. Thygesen, Phys. Rev. B 94, 155406 (2016).
35. J. C. Inkson, Surface Science 28, 69 (1971).
36. F. Hüser, T. Olsen,  and K. S. Thygesen, Phys. Rev. B 88, 245309 (2013).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters