Revised Periodic Boundary Conditions (RPBC)

In this study the periodicity has a central role. In general many structures that are studied in solid state physics are periodic. This way computations are much easier, because it is enough to simulate a single unit cell with suitable boundary conditions.

As described in the introduction, the unit cell of the phosphorene contains four atoms and it is periodic in in-plane directions. The method to handle periodicity of the structure was RPBC [15]. This approach is used by Hotbit calculator. Basically the theory is a new form of Bloch’s theorem. It represents symmetries as operators.

Usually Bloch’s theorem is written as [30, pp. 163-164]

ψ_a~k(~r−~T) = e^−i~k·T~ψ_a~k(~r). (54) This means that in an initial unit cell lies a point~r and a translation T~ is added to it. Translation moves the point in the periodic direction. This causes the original wavefunction to be multiplied with additional phase factor. In the end the wavefunction is the same as in position~r but its phase is changing due the periodicity From now on the derivation follows the references [15] and [31]. Next Kit et al. introduce another way to represent this transformation [15]. Translation can be denoted with operator τⁿ~r =~r+T~_n and corresponding inverse transformation τ⁻ⁿ~r = ~r +T~−n = ~r −T~_n. Then they use action ˆD(τⁿ) to describe the Bloch’s theorem as [15]

D(τˆ ⁿ)ψ_a~k(~r) = ψ_a~k(τ⁻ⁿ~r) =e^−i~k·T~nψ_a~k(~r). (55) In order to proceed into actual revised Bloch’s theorem authors point out that the electrons within a potential should be invariant in a general symmetry operation r~⁰ =Sⁿ~r. In other words the potential is symmetric (invariant in symmetry operation) as [15, 31]

D(Sˆ ⁿ)V(~r) =V(S⁻ⁿ~r) =V(~r). (56) The revised Bloch’s theorem works with any symmetries not just with translations.

Next they use a set of commuting symmetry operations S_iⁿⁱ by notation S^~n ≡ S₁ⁿ¹S₂ⁿ². . . . If symmetry operation is used for wavefunction, one gets

D(Sˆ ^~n)ψ_aκ(~r) =ψ_aκ(S^−~n~r) =e^−iκ·~nψ_aκ(~r). (57)

that tell the number of corresponding symmetry operationsS_iⁿⁱ. Symmetry operations D(Sˆ ^~n) commute with each other and with Hamiltonian [15]. Because of the symmetry, if the operation is done several times, at some point it returns to its original state.

This is represented as

D(Sˆ ^M~) = ˆD(S^M~1S^M~2. . .) = 1 (58) and as before with ~n now M~ = (M₁⁰M₂⁰. . .). On the other hand nowM_j⁰ = 0 or Mj. HereM_j is the number, which leads the operations back to the original system.

Following the original derivation next the unitarity of the symmetry operation is needed. In the reference [32, pp. 16-17] Rose points out that unitary operator can be represented uniquely using exponent of the Hermitian operator. This yields

D(Sˆ ^~n) = e^−iˆκ·~n, (59) which is essentially a constant. These leads to consider an eigenvalue problem

D(Sˆ ^~n)|ψi=e^−iˆκ·~n|ψi=C_~n|ψi. (60) Here C~n is a constant. It is convenient to operate equation (60) from the left with hκ|, which leads to κ-representation

hκ|D(Sˆ ^~n)|ψi=hκ|e^−iˆκ·~n|ψi=e^−iκ·~nψ(κ) =C_~nψ(κ). (61) Next the eigenvalue equation (60) is operated from left with h~r|. The aim is to form position representation of ψ same way as its κ representation was formed. In that case

h~r|D(Sˆ ^~n)|ψi=hS^−~n~r|ψi=e^−iκ·~nh~r|ψi (62) and then in position presentation

D(Sˆ ^~n)ψ(~r) = ψ(S^−~n~r) = e^−iκ·~nψ(~r). (63) Next step in the original derivation [15] is to use the cyclic conditions from equation (58). This leads to restrict the values ofκ which is seen as

e^iκjMj = 1 =e^i2πmj (64)

and then

κ= 2π

m₁ M₁, m₂

M₂, . . .

, (65)

where m_j = 0,1,2, . . . Mj−1. In the end κ is a good quantum number and revised Bloch’s theorem is proved [15]. This representation of Bloch’s theorem makes it possible to handle easily several different symmetries such as planes, tubes and spheres just by using suitable operators. These geometries are implemented into the Hotbit and all of them can be found in references [15] and [19].

3 Parametrization for DFTB and the basic prop-erties of phosphorene

The parametrization process has two main parameter types to determine. The first type is electronic parameters and the second type is repulsive energy parameters [33]. Some of these parameters are acquired directly from the computations. On the other hand repulsive energy parameters are adjusted mostly by hand. Theory behind the parameters is described in the theory part. Here the parameter determination process and actual parameters are presented. The parameters are listed in table 1.

The first electronic parameters are acquired directly in the beginning of the process. They are calculated by Hotbit according to its pre-set information about different atoms. They are the energies of the highest occupied orbitals (3s and 3p), Hubbard U and full width half maximum (FWHM) of the Gaussian profile of the Coulombic interaction term. After these we added the energy of the 3d orbitals which are unoccupied. This value was adjusted several times. Most importantly it affects to the size of the band gap.

After determining the parameters to represent individual atoms, the task was to define the orbital overlap into the Slater-Koster tables. Here the parameter r₀ was set. As explained in the theory part r₀ is used to create spacially restricted pseudoatoms and it is usually kept about twice the size of the covalent bond length (r_cov) [18, 26]. Length of the covalent bond is acquired directly from the calculations

so only the multiplier is needed.

The third step is to find suitable repulsive potential between the atoms. First some simple structures with constant bond length were scaled using full DFT and the repulsive forces were measured. This produced several data points where force was presented as a function of bond length. The structures were dimer, white phosphorus (tetrahedron), red phosphorus (chain of tetrahedrons), and black phosphorus. Next task was a curve fitting to these data points. Giving them suitable weights the position and the form of the curve was adjusted. The smoothness of the curve was controlled with parameter s. The distance where the repulsion became zero was

determined by r_cut parameter. The repulsive force curve was used to integrate the repulsive energy curve. These are presented in the figure 2.

In the figure 2 dimer has two different values to weight. The literature value for the bond length was set to 1.895 Å and the repulsion was calculated by Hotbit. The result was weighted with w_dimer. Another parameter for dimer was wdimer_curve. It weights the repulsion data points formed by stretching the bond. The repulsion in stretching was calculated with DFT.

Finding the parameters was done by repeating the same steps several times. First electronic parameters were set and then repulsive potential was formed. After this the unit cell was optimized and some basic properties were calculated. They were compared to the literature values. The properties were the dimensions of the unit cell, bond lengths, Young’s moduli, shear moduli, Poisson ratios, bending moduli, band gap and the overall band structure. After the comparison parameters were adjusted and the process was repeated.

Table 1. Only_3s, _3p, U and FWHM were acquired directly from the compu-tations. Here s determines the smoothness of the repulsion curve. Weights of different data points are w_i. Bohr radius is r_Bohr

Parameters Values

_3s −0.512163 eV _3p −0.206132 eV

_3d 0.262 eV

U 0.293863 Ha

FWHM 4.52019 r_Bohr r₀ 1.873· r_cov

r_cut 2.58 Å

s 5.9

wdimer 0.6

wdimer_curve 1.6

w_{W P} 0.7

w_RP 0.6

w_BP 1.6

1.9 2.0 2.1 2.2 2.3 2.4 2.5

Figure 2. Right side presents the repulsive forces. The datapoints are acquired from scaling of the corresponding structures with constant bond lengths using DFT. On the left there is a repulsive energy curve which is integrated from the repulsive force curve.

Gaus et al. have also formed a DFTB parametrization for phosphorus [33].

Our electronic parameters are in agreement. It has to be noted that 3d differs significantly, but this is not a concern. Gaus et al. have done parametrization for phosphorus in general case, but ours is designed especially for phosphorene. Also our other parameters have effects that overlap with the ones of _3d. The parametrization is always a compromise between different properties to model.

The parametrization is done so that it could be used to simulate phosphorene as well as possible in both mechanical and electronic cases. On the other hand the broadness causes some errors. The most significant one is the direct band gap. The gap should be essentially direct in Γ point but because of the parametrization the peak in the Γ-Y path is slightly higher than the peak in Γ. In this study we focus on the direct band gap in Γ. The other phenomena are observed qualitatively in

comparison to this. The initial band structure can be seen in the figure 4a).

The properties that were acquired from the final parametrization are listed in the tables 2, 3 and 4 with reference results from other studies. These are the same properties which were described earlier. We studied just a single layer of phosphorene, so the z dimension of the unit cell is not relevant. DFTB cannot be used to determine the z dimension of the unit if there are no interactions in that direction. Furthermore the interaction between two planes are van der Waals interactions, which would need a different parametrization [18].

In table 2 there are some reference values for the physical dimensions of the phosphorene. There are some clear trends in the values. Our parametrization yields results in the correct scale but they don’t fully agree with previous studies.

This happens, because we emphasized the elastic and electronic properties of the phosphorene in the parametrization. Bond angle α is large andβ is small compared to previous studies. Also bond lengths and unit cell dimensions differ as seen in the table 2. If these dimensions would be perfect the problem is that some more important elastic or electronic properties, such as the size of the band gap, would be flawed.

The most suitable reference for the elastic properties is the paper written by Verma et al. [14]. They used objective molecular dynamics (OMD) [16], which is related to revised periodic boundary conditions (RPBC) [15] used in this study.

They have computed the bending moduli, Young moduli, shear moduli and Poisson ratios in every direction in the plain of the phosphorene. Also some other results about elastic properties were used as a reference. These are presented in the table 3.

Results are in a good agreement. There are no such problems as with bond lengths and angles. Interestingly shear moduli of this study differ little bit depending on the direction but values computed by Verma et al. do not differ [14].

We computed that the direct band gap was 1.872 eV. Earlier it was mentioned that the band gap of phosphorene is determined to be up to about 2.0 eV [5] so the result is reasonable. It is worth mentioning that the band gap of single layer phosphorene is always larger than the band gap of the bulk black phosphorus, which is around 0.3 eV [5, 6, 34]. There are also numerous other studies where the band gap has been determined computationally. Some of them are listed in the table 4.

Most of the DFT methods seem to underestimate the size of the band gap. They are nearly always considerably smaller than experimental results. Despite this variation

results in the table 4.

Table 2. Here are listed the bond lengths, bond angles and the unit cell dimensions from the different references and from this study. Comparing these values it can be seen that the parametrization yields reasonable physical measures for the phosphorene.

Method d₁ [Å] d₂ [Å] α[^◦] β [^◦] x[Å] y [Å]

DFTB (this study) 2.220 2.236 98.903 101.309 3.374 4.236 DFT (PBE, PAW) [35] 2.22 2.26 95.9 104.1 3.298 4.627 DFT (GGA, PBE, DFT-D2)[36] 2.22 2.25 96.28 103.75 --

--TB [37] 2.224 2.244 96.34 102.09 3.314 4.376

DFT [38] 2.24 2.28 96.00 103.51 3.32 4.58

DFT (PBE, HSE06)[6] -- -- -- -- 3.35 4.62

Table 3. Elastic properties of phosphorene have been studied in several papers.

Here Y is Young’s modulus, S shear modolus, σ Poisson’s ratio and K is a bending modulus. Armchair direction is marked with ac and zigzag direction with zz.

Method Y[_m^N] S[_m^N] σ K_ac[·10⁻¹⁸J]

DFTB (this study)

ac 32.135 38.258 0.267 0.444

zz 108.131 37.583 0.898 1.662

OMD [14]

ac 33.3 31.2 0.31 0.524

zz 79.1 31.2 0.73 1.311

DFT (PBE, PAW) [35]

ac 48.796 − − −

zz 184.094 − − −

VFFM* [39]

ac 52.2 − − −

zz 85.4 − − −

analytic [40]

ac − − − 0.76949

zz − − − 1.3781

*Valence Force Field Model

Table 4. Here are several values for the band gap of the phosphorene. The band gap vary depending on the used method. Especially many DFT studies have acquired relatively small band gaps compared to the experimental results.

Method Band gap [eV]

DFTB (This study) 1.872

G₀W₀ [5] 2.0

G₁W₁ [5] 1.4

DFT (PBE, HSE06) [6] 1.0

DFT (GGA, PBE, DFT-D2)[36] 0.88

DFT [38] 1.51

DFT (PBE, GGA) [41] 0.91

DFT [42] 1.0

experimental [6] 1.45

experimental [8] 1.24

experimental [43] 2.05

4 Effects of elastic deformations to band struc-ture

The main focus of the study was the electronic band structure of the phosphorene and its behavior during stretching, shearing and bending. All of these were done in both armchair and zigzag direction but bending was also done in other chiral directions. The curvature radii in these computations are large. The bending is much more modest than in nanotubes.

Electronic structure of phosphorene is sensitive to physical distortions which is seen as changes in the band structure. Although one has to note the weaknesses of the used DFTB method. As mentioned earlier, due to the parametrization the initial band gap is not direct as it should be. This is not a big concern because we measured only the direct band gap along the Γ point. With this kind of an approximate method the relative changes in the band structure are more important than the absolute ones.

Studying the band structure the orthorhombic lattice was used. The observed path in the Brillouin zone was S-X-Γ-Y-S. It is visualized in the figure 3. The notations follow the paper of the Setyawan and Curtarolo [44]. In the figures, where band structures are imaged, G corresponds to Γ. The zero level of the band structure plots is set around the Fermi level. In some cases computations did not yield exactly the same Fermi level but values varied. In order to make plots easier to analyze the zero energy level was set to be about 4.0 eV. This value was kept constant.

4.1 Stretching scales direct band gap

The band structures, that were computed during the stretching, are presented in the figure 4. When phosphorene was stretched or compressed its band gap started to change from direct to indirect and the size of the direct band gap changed. When compression in armchair direction increased, the band gap moved away from direct Γ. HOMO state was on the peak close to Y point in Γ-Y path and LUMO was still on Γ point. The change in HOMO state is so large that it can be verified even with

X

S Γ

Y

k_x k_y

Figure 3. The path for observing the band structure within the orthorhombic Brillouin zone is similar to the one used in reference [44].

this parametrization. At the same time the actual band gap shrunk significantly as did the measured direct band gap. The phenomena are visible with just −6.87%

strain. The same type of a phenomena was also observed by Peng, Wei and Copple [41] but they report that strain, which causes these changes, is between −12% and

−8%. On the other hand, the actual stretching in ac-direction caused band gap to grow nearly as greatly as it shrunk in compression. Stretching kept band gap direct on the Γ point.

In zz-direction stretching caused band gap to change into an indirect one. In this case stretching in zz-direction made the lowest energy peak in the unoccupied states to move from Γ point towards the Y point. This happened when about 7% stretch was applied. This is the same behavior as reported by Peng, Wei and Copple [41].

Compression in zz-direction causes LUMO state to get closer to the HOMO. The band gap became smaller, but wether the transition from direct to indirect band gap happens is not clear with our parametrization.

The behavior of the direct band gap at Γ point of interest is visualized as a function of strain in figure 5. While studying the direct gap the maximum strain was ±11%. Compression made the gap to shrink and stretching made it to grow in both directions. That’s the only similarity. In armchair direction the gap behaves so that data points form a convex curve but in zigzag case they form concave curve as seen in figure 5. After about 6% stretching in zz-direction the gap began to shrink instead of growing. This behavior has been observed also by Peng, Wei and Copple for the real band gap [41] and by Qin et al. for the bulk phosphorus [34]. Later when the mulliken populations are studied one realizes that this is caused by the

S

Figure 4. Stretching and compression cause significant changes in the band structure. In a) there is the initial band structure of the optimized structure, in b) unit has been compressed 6.87 % in ac-direction, in c) it has been stretched 6.88 % in ac-direction. The rest are for the zz-direction: d) sheet is compressed 6.87 %, e) sheet is stretched 5.50 % and f) stretching is 6.88 %. In zz-direction it is seen that between 5.50% and 6.88% stretching the band gap transforms into indirect one. Density of states (DOS) is imaged for angular quantum numbers l= 0, 1 and 2. Blue line is a total DOS.

10 5 0 5 10 strain (%)

1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6

direct band gap (eV)

armchair zigzag

Figure 5. Stretching in ac- and zz-directions cause very different changes in the direct band gap. When strain is close to zero, both has overall the same trend about increase and decrease of the gap but the curving is not similar.

In document Density-functional tight-binding modeling of electromechanics of phosphorene (sivua 28-42)