# Enhanced thermoelectric performance of phosphorene by strain-induced band convergence

###### Abstract

The newly emerging monolayer phosphorene was recently predicted to be a promising thermoelectric material. In this work, we propose to further enhance the thermoelectric performance of phosphorene by the strain-induced band convergence. The effect of the uniaxial strain on the thermoelectric properties of phosphorene was investigated by using the first-principles calculations combined with the semi-classical Boltzmann theory. When the zigzag-direction strain is applied, the Seebeck coefficient and electrical conductivity in zigzag direction can be greatly enhanced simultaneously at the critical strain of 5% where the band convergence is achieved. The largest value of 1.65 at 300 K is then achieved conservatively estimated by using the bulk lattice thermal conductivity. When the armchair-direction strain of 8% is applied, the room-temperature value can reach 2.12 in the armchair direction of phosphorene. Our results indicate that strain-induced band convergence could be an effective method to enhance the thermoelectric performance of phosphorene.

###### pacs:

73.50.Lw, 73.61.Cw, 73.22.-f, 71.15.Mb## I Introduction

Thermoelectric materials, which can directly convert heat into electricity and vice versa, have attracted much interest from the science community due to the current critical energy and environmental issues. The performance of a thermoelectric material is quantified by the dimensionless figure of merit , where is the Seebeck coefficient, is the electrical conductivity, is the absolute temperature, and are the electronic and lattice thermal conductivity, respectively. A good thermoelectric material should have large Seebeck coefficient and electrical conductivity, and/or low thermal conductivity. It is challenging to achieve a high value since optimizing one transport coefficient often leads to another adversely affected, which hinders the wide applications of thermoelectric materials. It is very important to find methods to solve such problem.

NanostructuringMS1-1993 (); MS2-1993 () and band convergenceO.Rabin-2001 () were suggested to be two promising solutions. The two- or one-dimensional structures could have much larger values than their bulk counterparts, due to the enhanced power factor () caused by a sharper density of states near the Fermi level, or the reduced lattice thermal conductivity caused by the increased phonon scattering. On the other hand, it was found that the conduction or valence band extrema can be modulated and converged by tuning the doping and composition.Y.Pei-Nature (); W.Liu-PRL (); X.J.Tan-PRB () Such so called band convergence can significantly enhance the Seebeck coefficient without detrimental effects on the electrical conductivity. If we can choose an appropriate material and apply the two methods in it, good thermoelectric performance can be expected.

Recently, the single layer of black phosphorus (black-P), the so called phosphorene, was successfully synthesized.Y.B.Zhang-2014 (); H.Liu-ACS nano (); W.Lu-Nano-research () Black-P is an elemental solid which is the most stable form among the phosphorus allotropes under normal condition.T.Nishii-1987 () Similar to graphite, black-P crystallizes in a layered structure, namely, each phosphorus atom is covalently connected to three neighboring phosphorus atoms to form a puckered layer. It is a direct-gap semiconductor with a band gap of about 0.33 eV.R.W.Keyes-1953 () When the black-P is exfoliated into few or even single layer, extraordinary optoelectronic and electronic properties emerge.Y.B.Zhang-2014 (); H.Liu-ACS nano (); F.Xia-arXiv (); J.Qiao-arXiv (); D.F.Shao-arXiv (); V.Tran-arXiv () A much larger direct band gap of about 2 eV was found in monolayer phosphorene.V.Tran-arXiv () In particular, the nanostructuring can largely enhance the thermoelectric performance: the value of monolayer phosphorene is much larger than that of bulk black-P.H.Y.Lv-arXiv (); R.Fei-arXiv () This implies the nanoscale phosphorene-based material can be a good candidate of the applicable thermoelectric material. Very recently, it was reported that the electronic band structure of phosphorene can be tuned by the in-planeR.Fei-Nanolett (); X.Peng-arXiv () or out-of-plane strains.A.S.Rodin-PRL () This inspires us to investigate whether the band convergence can be introduced into the system by the method of strain. If so, the thermoelectric performance should be further optimized.

In this work, we demonstrate that we can direct the conduction band convergence of phosphorene under a simple tensile strain condition, which will greatly enhance the thermoelectric performance of phosphorene. The effect of the uniaxial strain on the thermoelectric properties of phosphorene is investigated using the first-principles calculations combined with the semi-classical Boltzmann transport theory. Our results show that with a moderate tensile strain, the conduction band extrema of phosphorene can be converged, which results in an increase in the Seebeck coefficient. At the same time, the electrical conductivity at a particular direction is dramatically increased and therefore the largely increased power factor is obtained. When the zigzag-direction strain is applied, the largest value of 1.65 at 300 K is obtained in the zigzag direction of phosphorene at the critical strain of 5%, conservatively estimated by using the bulk lattice thermal conductivity. The room temperature value can reach 2.12 in the armchair direction of phosphorene under an 8% armchair-direction strain.

## Ii Computational details

The structural and electronic properties of phosphorene are investigated using the first-principles pseudopotential method as implemented in the ABINIT code.X.Gonze-2009 (); X.Gonze-2002 (); X.Gonze-2005 () The Brillouin zones are sampled with a Monkhorst-Pack -mesh. The cutoff energy for the plane wave expansion is set to be 800 eV. For the structural relaxation, the exchange-correlation energy is in the form of Perdew-Burke-Ernzerhof (PBE)PBE-1996 () with generalized gradient approximation (GGA). Both the geometries and atomic positions are fully relaxed until the force acting on each atom is less than . In the calculations of the electronic structures, the TB-mBJ potentialTB-mBJ () is used, which can reproduce accurate band gaps for many semiconductors. Based on the electronic structure, the electronic transport coefficients are derived by using the semi-classical Boltzmann theory within the relaxation time approximationBoltz () and doping is treated by the rigid band model.Rigid-band () The electronic thermal conductivity is calculated using the Wiedemann-Franz law , where is the Lorenz number. In this work, we use the Lorenz number of .Franz-law ()

## Iii Results and Discussion

Figure 1 shows the structure of phosphorene, with the top and side views illustrated in Figs. 1(a) and (b), respectively. The dashed line denotes the primitive cell used in our calculation, with the corresponding first Brillouin zone inserted in the figure. To confirm the reliability of our method, we first do the calculation for bulk black-P, which has the experimental results to compare with. The van der Waals interactions between the neighboring layers of bulk black-P are treated by the vdW-DFT-D2 functional.vdW-DFT () The calculated lattice parameters are =3.34 Å, =10.51 Å and =4.43 Å, which are very close to the experimental values.L.Cartz-1979 () The bulk black-P is semiconducting with a direct band gap of 0.34 eV at the point, in good agreement with those reported experimentally.R.W.Keyes-1953 (); D.Warschauer-1963 () Our calculation confirms that the TB-mBJ potential can accurately predict the band gap of our investigated system, which is an important factor in determining the electronic transport properties. In the following, we use the same method to deal with the calculations for phosphorene. The lattice parameters of phosphorene are calculated to be =3.32 Å and =4.63 Å. The strain is applied along the zigzag or armchair direction of phosphorene, as indicated in Fig. 1(a).

### iii.1 Strain applied along the zigzag direction

#### iii.1.1 Energy band structure

First, we discuss the case when the tensile strain is applied along the zigzag direction.

Figure 2 shows the evolution of the band structure of phosphorene under the uniaxial strain along the zigzag direction. When no strain is applied (see Fig. 2(a)), the phosphorene is a direct-band-gap semiconductor with a gap of 1.80 eV located at the point. The band gap is much larger than that calculated using DFT-PBE method (0.92 eV)R.Fei-Nanolett () but is very close to the value calculated by the method (2.0 eV).V.Tran-arXiv () The conduction band minimum (CBM) of phosphorene is highlighted in the red color. Note that there exist the other three band extrema located at the point and along the - and - directions, respectively, highlighted in the blue color. The four band extrema are denoted by the symbols “\@slowromancapi@”, “\@slowromancapii@”, “\@slowromancapiii@”, and “\@slowromancapiv@” respectively. The applied zigzag-direction strain mainly affects the conduction bands of phosphorene. As increasing the strain, the band extremum “\@slowromancapiii@” is gradually lowered while the band extremum “\@slowromancapi@” is elevated gradually. As for the band extrema “\@slowromancapii@” and “\@slowromancapiv@”, there exists a critical strain, that is, 5%. When the strain is smaller than 5%, the band extremum “\@slowromancapii@” is elevated while “\@slowromancapiv@” is lowered as increasing the value of strain. The band extrema “\@slowromancapii@” and “\@slowromancapiv@” reach their maximum and minimum respectively at the strain of 5%. When further increasing the strain, the band extremum “\@slowromancapii@” is however lowered and “\@slowromancapiv@” is elevated gradually. The different behavior of the four band extrema leads to the band convergence at the strain of 5%, which is denoted by the dashed green line in Fig. 2(c). The convergence of the band extrema will in turn lead to the increase in the Seebeck coefficient, which will be discussed later. Moreover, when the strain is smaller than 5%, the CBM of phosphorene locates at the point (see Figs. 2(a) and (b)); when the strain reaches 5%, the band extremum “\@slowromancapiii@” becomes energetically lower than band extrema “\@slowromancapi@”, “\@slowromancapii@”, and “\@slowromancapiv@”, and thus the transition of direct-indirect band gap occurs. When the doping level in the system is low, the electrical conductivity is dominated by the CBM, so we can expect that at the critical strain of 5%, the electrical transport property will be changed dramatically. We will come back to this point later.

#### iii.1.2 Electronic transport coefficients

Based on the calculated band structure, the electronic transport coefficients of phosphorene can be evaluated by using the semi-classical Boltzmann theory and rigid band model. To get reliable results, a very dense mesh up to 840 points in the irreducible Brillouin zone (IBZ) was used. Within this method, the Seebeck coefficient can be calculated independent of the relaxation time ; however, the electrical conductivity can only be calculated with inserted as a parameter, that is, what we obtain is /. The relaxation time is determined by the formula , where is the carrier mobility and is the band effective mass. The effective mass tensor is defined as .

The mobility of phosphorene is calculated using the deformation potential (DP) theory on the basis of the effective mass approximation:J.Bardeen-1950 (); P.J.Price-1981 (); J.Xi-2012 ()

(1) |

where is the temperature. is the elastic modulus and for the 2 system, the in-plane value is defined as , where , and are, respectively, the total energy, the applied uniaxial strain and the area of the investigated system. The DP constant along a certain direction is obtained by , where is the energy of the band edges (valence band maximum for the holes and conduction band minimum for the electrons) and is the applied small strain (by a step of 0.5%).

The elastic modulus can be obtained by fitting the curve of energy-strain relationship. The calculated value of in the zigzag direction is 106.18 N/m. The value is smaller than that of (about 120 N/m)R.C.Cooper-2013 () and graphene (about 335.0 N/m).C.Lee-2008 () The DP constant is obtained by the linear fitting of the energy shift of band edges as a function of the strain. The calculated is 3.98 eV for electrons. When the uniaxial zigzag-direction strain is applied, the effective mass as well as the corresponding mobility at 300 K in the zigzag and armchair directions are shown in Figs. 3(a) and (b), respectively. Since the strain mainly modulates the bands above the Fermi level, we only focus on the properties of electrons. We can see that for the zigzag direction (see Fig. 3(a)), there is a sharp drop of the effective mass at the critical strain (5%), which in turn leads to the dramatically increased carrier mobility, with the values about two orders of magnitude larger than those before the critical strain. Note that here the effective mass and mobility are calculated with respect to the CBM (highlighted in red color in Fig. 2). However, at the critical strain of 5%, the energy difference of the four conduction band extrema is very small, thus except for the CBM, the other three band extrema will also contribute to the electrical conductivity. The values of the carrier mobility of band extrema “\@slowromancapi@”, “\@slowromancapii@” and “\@slowromancapiv@” are denoted by the red stars in the Fig. 3. The larger the size of the star, the more it will contribute to the electrical conductivity. After the critical strain, the greatly enhanced mobility and therefore the electrical conductivity will be very beneficial to the thermoelectric application. The trend of the mobility along the armchair direction (see Fig. 3(b)) is just reversed. However, although there is a sharp drop of the mobility at the critical strain, the mobilities after that point are still considerable, larger than 1000 . The different behavior of the effective mass and carrier mobility in the two different directions upon the applied zigzag-direction strain will in turn influence their thermoelectric properties, which will be discussed later.

Strain | 0% | 2% | 4% | 5% | 6% | 7% | 8% | 10% | |
---|---|---|---|---|---|---|---|---|---|

Zigzag | () | 1.246 | 1.274 | 1.302 | 0.148 | 0.150 | 0.151 | 0.160 | 0.170 |

() | 61.46 | 58.79 | 56.24 | 4356.09 | 4240.65 | 4184.72 | 3727.14 | 3301.55 | |

(fs) | 43.53 | 42.57 | 41.62 | 366.45 | 361.56 | 359.17 | 338.96 | 319.02 | |

Armchair | () | 0.227 | 0.251 | 0.282 | 0.473 | 0.471 | 0.445 | 0.460 | 0.440 |

() | 5149.24 | 4211.61 | 3336.55 | 1185.97 | 1196.06 | 1339.91 | 1253.95 | 1370.53 | |

(fs) | 664.39 | 600.86 | 534.81 | 318.85 | 320.20 | 338.91 | 327.86 | 342.76 |

Based on the calculated effective mass and carrier mobility , we are now able to obtain the relaxation time . In Table I, we summarize the relaxation time at 300 K along the two different (zigzag and armchair) directions of phosphorene when the uniaxial zigzag-direction strain is applied, with the corresponding effective mass and carrier mobility included. We can see that before the critical strain of 5%, the relaxation times along the armchair direction are much larger than those along the zigzag direction. However, they become comparable when the strain reaches 5%.

Inserting the calculated relaxation time , we plot in Figs. 4(a) and (b) the room-temperature electrical conductivity as a function of the carrier concentration along the zigzag and armchair directions, respectively, when the different zigzag-direction strains are applied. When no strain is applied, the electrical conductivity exhibits obvious anisotropic property, with the value along the armchair direction much larger than that along the zigzag direction. This anisotropic property is due to the different dispersions of the CBM along the - (zigzag direction in the real space) and - (armchair direction in the real space) directions. The band along the - direction is much flatter than that along the - direction, which results in the much larger band effective mass and therefore much smaller carrier mobility and electrical conductivity in the zigzag direction. For the transport along the zigzag direction, after the critical strain (5%), the electrical conductivity is obviously enhanced. If we notice the band structure of phosphorene under the critical strain (Fig. 2(c)), we can see that the conduction band minimum (the red line) which dominates the electrical conductivity moves from the point to the one along the - direction. The effective mass of that band in the zigzag direction is dramatically decreased (see Fig. 3(a) and Table I), which results in the largely increased electrical conductivity. However, in the strain range of , the electrical conductivity under the strain of 5% is the smallest (see Fig. 4(a)), which is consistent with our explanation above that although the mobility of the CBM is the largest at the 5% strain, the other three band extrema “\@slowromancapi@”, “\@slowromancapii@” and “\@slowromancapiv@” will also contribute to the electrical conductivity, which will more or less decrease at that strain. The electrical conductivity along the armchair direction is however decreased by the strain due to the increased effective mass and therefore the decreased carrier mobility in this direction, as shown in Fig. 3(b).

Figures 4(c) and (d) show the Seebeck coefficient at 300 K as a function of the carrier concentration along the zigzag and armchair directions, respectively, under the zigzag-direction strain condition. For both directions, the absolute values of the Seebeck coefficients reach the maxima at the critical strain of 5%, where the band convergence is achieved as indicated in Fig. 2(c). A high Seebeck coefficient is often caused by a high overall density-of-states effective mass . The is related to the band effective mass (i.e., the effective mass of a single carrier pocket) by , where is the number of the degenerated band valleys. A large value of will lead to a large , and therefore a large absolute value of the Seebeck coefficient. When the band extrema of the multiple bands have no or little difference in energy (i.e., the bands are converged), they are considered to be effectively degenerated and is the increased. Here, we obtain the band convergence by applying a uniaxial strain and a significantly increased Seebeck coefficient is achieved at the critical strain. It has been reported that the convergence of many valleys in bulk materials could be achieved by tuning the doping and composition.Y.Pei-Nature (); W.Liu-PRL (); X.J.Tan-PRB () Note here that for the zigzag direction and at the critical strain, not only the Seebeck coefficient reach the maximum value, but also the electrical conductivity is greatly increased, which will benefit the power factor of phosphorene.

The as a function of the carrier concentration in the zigzag and armchair directions of phosphorene at 300 K are plotted in Figs. 4(e) and (f), respectively, when the zigzag-direction strain is applied. We can see that for the zigzag direction (Fig. 4(e)), the of phosphorene is greatly enhanced when the strain reaches 5%. At relatively higher carrier concentration, the phosphorene under 5% zigzag-direction strain has the largest , due to the largest Seebeck coefficient and greatly enhanced electrical conductivity . However, for the armchair direction, although the phosphorene under 5% strain has the largest Seebeck coefficient, the electrical conductivity is greatly decreased, which leads to the decrease of the . The electrical conductivity exhibits the obvious anisotropic property, which results in the different behavior of the between the zigzag and armchair directions upon the applied zigzag-direction strain.

#### iii.1.3 Dimensionless figure of merit

The electronic thermal conductivity of phosphorene is calculated based on the Wiedemann-Franz law as mentioned above. We then estimate the values of phosphorene at 300 K under the zigzag-direction strain by using the experimental lattice thermal conductivity of bulk black-P (12.1 ). The values at 300 K along the zigzag and armchair directions are plotted as a function of the carrier concentration in Figs. 5(a) and (b), respectively. We can see that the value of phosphorene without strain exhibits the strong anisotropic property, with the value along the armchair direction much larger than that along the zigzag direction. The applied zigzag-direction strain has different effect on the values in the two particular directions. For the armchair direction (see Fig. 5(b)), the value is decreased by the strain. The largest value of 1.44 can be obtained in this direction when no strain is applied and the corresponding carrier concentration is . For the transport along the zigzag direction, however, the value is greatly enhanced by the strain. The largest value of 1.65 is achieved under the zigzag-direction strain of 5%, with the corresponding carrier concentration of . The maximum of the value is 50 times larger than that without strain. Applying strain is shown to be an effective way to tune the band structure of phosphorene and the thermoelectric performance in particular direction can be largely optimized. Note here that we use the lattice thermal conductivity of bulk black-P to estimate the values of phosphorene. If the lattice thermal conductivity of the phosphorene can be reduced compared with the bulk value, which can be realized in many low-dimensional structures, the value of the phosphorene can be further enhanced.

### iii.2 Strain applied along the armchair direction

We also apply the strain along the armchair direction of phosphorene. The same as the case when the zigzag-direction strain is applied, the band structure and the thermoelectric properties of phosphorene can also be tuned by the strain. In this case, when the armchair-direction strain reaches 8%, the band extrema convergence can be achieved, as indicated by the dashed green line in the inset in Fig. 6(a). However, the number of the converged bands is less than that when the zigzag-direction strain is applied. As a result, the enhancement of the value is not as large as the case when the stain is applied along the zigzag direction. The largest values at 300 K in the zigzag and armchair directions are respectively 7.4 and 0.5 times larger than the corresponding values when no strain is applied. At the critical strain of 8%, we can obtain the largest room-temperature values of 0.27 and 2.12 in the zigzag and armchair directions of phosphorene, respectively.

## Iv Conclusion

In summary, we have investigated the strain effect on the electronic and thermoelectric properties of phosphorene based on the the first-principles calculations combined with the semi-classical Boltzmann theory. The band structure of phosphorene can be modulated by the uniaxial strain and conduction band extrema are effectively converged at the critical strain. The Seebeck coefficient is largely enhanced due to the band convergence. The electrical conductivity exhibits anisotropic property and behaves differently in the zigzag and armchair directions of phosphorene upon the applied strain. When the zigzag-direction strain is applied, the Seebeck coefficient and electrical conductivity in zigzag direction can be greatly enhanced simultaneously at the critical strain of 5%. The largest value of 1.65 is then achieved, which is 50 times larger than that without strain. When the armchair-direction strain of 8% is applied, we can obtain the largest value of 2.12 in the armchair direction of phosphorene. Our results indicate that strain-induced band convergence could be an effective method to enhance the thermoelectric performance of phosphorene.

## V Acknowledgement

This work was supported by the National Key Basic Research under Contract No.2011CBA00111, the National Nature Science Foundation of China under Contract No.11274311, the Joint Funds of the National Natural Science Foundation of China and the Chinese Academy of Sciences’ Large-scale Scientific Facility (Grand No.U1232139), and Anhui Provincial Natural Science Foundation under Contract No.1408085MA11. The calculation was partially performed at the Center for Computational Science, CASHIPS.

## References

- (1) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 12727 (1993).
- (2) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 16631 (1993).
- (3) O. Rabin, Y.-M. Lin, and M. S. Dresselhaus, Appl. Phys. Lett. 79, 81 (2001).
- (4) Y. Pei, X. Shi, A. LaLonde, H. Wang, L. Chen, and G. J. Snyder, Nature 473, 66 (2011).
- (5) W. Liu, X. J. Tan, K. Yin, H. J. Liu, X. F. Tang, J. Shi, Q. J. Zhang, and C. Uher, Phys. Rev. Lett. 108, 166601 (2012).
- (6) X. J. Tan, W. Liu, H. J. Liu, J. Shi, X. F. Tang, and C. Uher, Phys. Rev. B 85, 205212 (2012).
- (7) L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nat. Nanotechnol. 9, 372 (2014).
- (8) H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, ACS Nano 8, 4033 (2014).
- (9) W. Lu, H. Nan, J. Hong, Y. Chen, C. Zhu, Z. Liang, X. Ma, Z. Ni, C. Jin, and Z. Zhang, Nano Research (in press, 2014).
- (10) T. Nishii, Y. Maruyama, T. Inabe, and I. Shirotani, Synth. Met. 18, 559 (1987).
- (11) R. W. Keyes, Phys. Rev. 92, 580 (1953).
- (12) F. Xia, H. Wang, and Y. Jia, arXiv: 1402.0270v2.
- (13) J. Qiao, X. Kong, Z.-X. Hu, F. Yang, and W. Ji, arXiv:1401.5045v3.
- (14) D. F. Shao, W. J. Lu, H. Y. Lv, and Y. P. Sun, arXiv: 1405.0092.
- (15) V. Tran, R. Soklaski, Y. Liang, and L. Yang, arXiv:1402.4192v3.
- (16) H. Y. Lv, W. J. Lu, D. F. Shao, and Y. P. Sun, arXiv:1404.5171.
- (17) R. Fei, A. Faghaninia, R. Soklaski, J.-A. Yan, C. Lo, and L. Yang, arXiv:1405.2836.
- (18) R. Fei and L. Yang Nano Lett. 14, 2884 (2014).
- (19) X. Peng, A. Copple, and Q. Wei, arXiv:1403.3771v2.
- (20) A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Phys. Rev. Lett. 112, 176801 (2014).
- (21) X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, T. Deutsch, L. Genovese, Ph. Ghosez, M. Giantomassi, S. Goedecker, D. R. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, M. J. T. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, D. Sangalli, R. Shaltaf, M. Torrent, M. J. Verstraete, G. Zerah, and J. W. Zwanziger, Comput. Phys. Commun. 180, 2582 (2009).
- (22) X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty, and D. C. Allan, Comput. Mater. Sci. 25, 478 (2002).
- (23) X. Gonze, G.-M. Rignanese, M. Verstraete, J.-M. Beuken, Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zerah, M. Mikami, Ph. Ghosez, M. Veithen, J.-Y. Raty, V. Olevano, F. Bruneval, L. Reining, R. Godby, G. Onida, D. R. Hamann, and D. C. Allan, Z. Kristallogr. 220, 558 (2005).
- (24) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (25) F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
- (26) G. K. H. Madsen and D. J. Singh, Comput. Phys. Commun. 175, 67 (2006).
- (27) E. A. Stern, Phys. Rev. 157, 544 (1967).
- (28) R. Venkatasubramanian, E. Siivola, T. Colpitts, and B. OQuinn, Nature 413, 597 (2001).
- (29) S. Grimme, J. Comput. Chem. 27, 1787 (2006).
- (30) L. Cartz, S. R. Srinivasa, R. J. Riedner, J. D. Jorgensen, and T. G. Worlton, J. Chem. Phys. 71, 1718 (1979).
- (31) D. Warschauer, J. Appl. Phys. 34, 1853 (1963).
- (32) J. Bardeen and W. Shockley, Phys. Rev. 80, 72 (1950).
- (33) P. J. Price, Ann. Phys. 133, 217 (1981).
- (34) J. Xi, M. Long, L. Tang, D. Wang, and Z. Shuai, Nanoscale 4, 4348 (2012).
- (35) R. C. Cooper, C. Lee, C. A. Marianetti, X. Wei, J. Hone, and J. W. Kysar, Phys. Rev. B 87, 035423 (2013).
- (36) C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science 321, 385 (2008).