Effective Hamiltonians for phosphorene and silicene

We derived the effective Hamiltonians for silicene and phosphorene with strain, electric ﬁ eld and magnetic ﬁ eld using the method of invariants. Our paper extends the work of Geissler et al 2013 ( New J. Phys. 15 085030) on silicene, and Li and Appelbaum 2014 ( Phys. Rev. B 90 , 115439) on phosphor-ene. Our Hamiltonians are compared to an equivalent one for graphene. For silicene, the expression for band warping is obtained analytically and found to be of different order than for graphene. We prove that a uniaxial strain does not open a gap, resolving contradictory numerical results in the literature. For phosphorene, it is shown that the bands near the Brillouin zone center only have terms in even powers of the wave vector. We predict that the energies change quadratically in the presence of a perpendicular external electric ﬁ eld but linearly in a perpendicular magnetic ﬁ eld, as opposed to those for silicene which vary linearly in both cases. Preliminary ab initio calculations for the intrinsic band structures have been carried out in order to evaluate some of the k p · parameters.


Introduction
Graphene is an interesting electronic material because of its two-dimensionality, zero-gap, and linear dispersion near the Fermi energy [1].These properties differ from standard three-dimensional semiconductors.The study of other two-dimensional materials beyond graphene came about quickly, with BN and MoS 2 being two early choices.More recent candidates include silicene [2] and phosphorene [3], though earlier studies of these two materials exist [4,5].
Nowadays, much of the study of electronic properties relies on ab initio calculations, in particular densityfunctional theory (DFT).DFT is very powerful in its ability to predict ground-state properties such as structures fairly well.It is less reliable and much more computationally intensive to calculate excited-state properties such as optical and transport properties; the inclusion of an external magnetic field is still a difficult problem and not implemented in standard DFT codes.Empirical tight-binding models have been early alternatives to DFT; indeed, the linear dispersions for both graphene [6] and silicene [2] were first identified using tightbinding models, even though the linear dispersion for silicene was probably first obtained (but not identified) in 1994 using ab initio calculations [4].
However, the most versatile and physically-transparent band-structure model is the k p • model [7], often giving energy bands analytically in the vicinity of extrema in terms of meaningful parameters such as effective masses and optical matrix elements.For graphene, this is the Dirac Hamiltonian linear in the wave vector, though nonlinear contributions have also been worked out [8].For silicene, only the Dirac Hamiltonian, with linear electric field terms, has been written down [9].For phosphorene, a couple of k p • models have recently appeared [10,11].Li and Appelbaum [11], in particular, did a very careful study of the band structure using perturbation theory near the Γ point, identifying the contributions to the effective masses; however, they again did not include external fields.
In the current work, we are therefore concerned with the most general k p • Hamiltonian for silicene and phosphorene, particularly in the presence of external fields since this problem has not been considered fully.We employ the method of invariants to derive such Hamiltonians.The primary goal is to identify possible terms not considered in previous work.A secondary goal is to thus obtain a better understanding of the band structures of silicene and phosphorene.

Method of invariants
Here we apply the method of invariants [12] to the materials of interest to us.Since most of the DFT calculations to date have not included relativistic effects and we are considering low atomic number elements, we have left out spin effects in this work; the main consequence in the band structure is the neglect of the small spin-orbit splitting of bands (e.g., for silicene, this is a 1.55 meV splitting at the K point [13]).The symmetries of interest are the space-group symmetries and time-reversal symmetry.Space-group symmetry only requires one to study the character tables of the point-group symmetry, whereas the consideration of time-reversal symmetry requires an application of the Herring test [12].Details of the theory are not repeated here as we have applied the formalism exactly as presented in [12].

Silicene
The point of interest for silicene is the K point in the first Brillouin zone (figure 1) since that is where the valence and conduction bands touch with a linear dispersion.In table 1, the group of the wave vector at point K for silicene is shown.We use the expressions x y x y and similarly for other physical quantities.It can be seen that all the irreducible representations are either onedimensional or two-dimensional, implying that the energy bands at K are all either non-degenerate or doubly degenerate.The Fermi level is at a point of double degeneracy.Tables 1-3 allow us to determine the Hamiltonian from combinations of single group functions based on symmetry principles formed.The next step is to apply the Herring test [12].It follows from the character table of where n = 6 is the order of the group D 3 , g is a representation of a group element, and χ is the character of that group element.Further, since a spatial symmetry  exists such that  − = k kwhere k is the K point (and ≠ − k k), all irreducible representations belong to the case a2 [12].In this case, time-reversal symmetry additionally requires all Hamiltonian terms to satisfy where  , , and  are the time-reversal operator, a Γ 1 invariant term, and a spatial symmetry operator establishing a link between basis functions at k and −k, respectively.In the case of silicene,  = 1 and   = x can be chosen where  x denotes reflection in the x axis using a coordinate system where the x axis passes through Table 1.Group of wave vector D 3 at the K point (silicene).k i are components of the wave vector, ϵ ij of the strain tensor, B i of an external magnetic field, E i of an external electric field, and J i of an angular momentum operator.
the center of an edge in the honeycomb lattice formed by the silicon atoms (referring to the same choice of coordinate system as in [14]).The symbol ζ takes the value +1 and −1 for even and odd functions under timereversal symmetry, respectively.
With the preceding analysis, we can write down the most general Hamiltonian for silicene allowed by symmetry.The result is, to leading orders, ) ) where the dots refer to higher-order terms.a i , b i , c i and e i are k p • parameters which are undetermined in the method of invariants.The H i terms are intrinsic band-structure Hamiltonians while the other ones exist in the presence of external fields.In the above Hamiltonian, the J i matrices represent the pseudospin degree of freedom and are the × (2 2) Pauli spin matrices.The difference between the work of Geissler et al [9] and the current one is that they were focussed on the interplay between spin-orbit coupling and an external electric field (and, hence, in the topological properties) at the linear level (in k and E z ), whereas we are more interested in understanding the basic band structure of the material beyond the linear terms and in the presence of other external fields.Thus, the a 2 and a 3 terms provide quadratic in k contributions while the a 4 and a 5 are of cubic order but only the a 4 term gives rise to an anisotropic term.Hence, one can readily say that the band structure of silicene to linear and quadratic orders is isotropic and anisotropic effects only manifest themselves if cubic terms become important.A discussion of the band dispersion is provided in the next section.
For comparison, we reproduce the corresponding  i for graphene [8]: The sign difference in the a 61 term is due to a different choice of phase while the a 62 is the corresponding anisotropy term for graphene.It is seen here to be quadratic in k though it is cubic in k for silicene.

Phosphorene
Phospherene and its bulk counterpart black phosphorus have the same in-plane translational symmetry and the nonsymmorphic space group is base-centered orthorhombic [11] (figure 2).
3 ).The notations Γ 3 1 and Γ ′ 3 2 denote the first component of the first Γ 3 func- tion and the second component of the second Γ 3 function, respec- tively, etc.
2 etc.Since the irreducible representations are all one-dimensional, the coupling constants of product representations are trivial.Further, since = − = k k 0 at the Γ point and the condition in equation ( 2) is satisfied by the irreducible representations, all irreducible representations belong to the case a1 [12].
We are now in a position to write down the most general Hamiltonian for phospherene allowed by symmetry.The result is ) Our intrinsic Hamiltonian agrees with the result of Li and Appelbaum [11].

Band structures
We now provide some discussion of the band structures that can be inferred from the Hamiltonians.We are also more interested in revealing the properties analytically than purely numerically using DFT calculations.However, preliminary calculations of the intrinsic band structures have been performed.The DFT-based calculations were conducted using the projector augmented-wave method [15] and the PBE-GGA exchangecorrelation functional [16] as implemented in the VASP code [17][18][19].The electronic wavefunctions were described using a plane-wave basis set with an energy cutoff of 500 eV for silicene and 1200 eV for phosphorene.Atomic positions were fully relaxed in Γ-centered 25 × 25 × 1 for silicene (9 × 9 × 1 for phosphorene) supercells for k-point sampling until residual forces were lower than 5 meV Å −1 .For silicene, the lattice constant obtained was 3.867 Å.For phosphorene, the optimized lattice constants were 3.300 and 4.624 Å.In both cases, a large interlayer separation in the z direction was chosen (∼30 Å) to minimize interactions.

Silicene
We are mainly concerned with the doubly-degenerate Dirac band since the Fermi level crosses it; it is made of p z and s orbitals from both atoms in the unit cell of silicene [2].

Intrinsic
For an arbitrary direction and keeping terms up to cubic in k, the band structure is given by diagonalizing the 2 , one gets Thus, our equation above shows that the intrinsic band dispersion can have quadratic, cubic, … terms in the wave vector.
The a 1 coefficient is related to the Fermi velocity, a 2 is related to an effective mass, a 4 gives rise to anisotropy, and a 5 leads to cubic terms.Along the Γ → K direction (k y = 0), the band structure simplifies to This reproduces the well-known result for graphene that the Fermi speed is the same for the two bands (figure 3, left).Note that the anisotropy parameter a 4 is absent along this direction.
We have performed a cubic fit of the bands obtained using VASP up to about half-way from K. Rewriting equation (19) as and performing a fit to where n is chosen to be 100 points between K and G, we obtain (with Å.The results of the curve fitting are given in table 6. Figure 3 (right) shows the dispersion along → K M 2 (symbols).It can be seen that the linear dispersion (line) is a good approximation only up to about a third of the way.

Strain
In the presence of strain, the Hamiltonian at the K point is given by First, we note that applying strain perpendicular to the plane (ϵ zz ) does not change the symmetry at the K point and a zero gap is preserved but there is nevertheless an expected shift in the energy levels.Second, an in-plane uniaxial strain (e.g., ϵ xx ), while changing the symmetry, leads to the same energy change for both bands making up the Dirac point since they both have the same deformation potential e 1 ; thus, the only effect is a shift in energies.Our analysis provides a formal resolution to the controversy from DFT calculations about whether a gap is opened by a uniaxial strain or not [20][21][22][23].Zhao and Mohan et al both predicted a gap opening [20,21].However, Qin et al [24] and Yang et al [23] did not obtain a gap and only obtained a shift of the Dirac point.The latter interpreted the disagreement of Zhao to their using insufficient k points near the crossing.However, the influence of strain does significantly change the dispersion at finite k values, i.e., in the Fermi velocity of those bands.Table 6.Band fitting for silicene along K to Γ.For finite k and an arbitrary strain, one can also diagonalize the Hamiltonian to give

Electric field
In the presence of an electric field, the Hamiltonian at the K point is given by An electric field in the z-direction will change the symmetry for silicene (as opposed to graphene which has all the atoms in the plane) and, therefore, a gap will open.Earlier DFT calculations have shown that the gap opening is initially linear in the field (hence ≠ c 0 1 ) and later becomes nonlinear.This behaviour is reproduced by our symmetry analysis.Indeed, by comparing to the DFT calculations of Drummond et al [25], one obtains = c 0.037 (bottom, in arbitrary units).Note that the plots are relative energies as they do not include the solution of the equations in the y direction which requires numerically solving the Airy equation.

Magnetic field
Finally, we present the behaviour of the energies in the presence of a magnetic field.First, the Hamiltonian linear in the magnetic field is Thus, at the K point, we have Figure 4. Left: two plots of dispersion relations near the K point for an applied electric field using the expression for  E in equation (8)  and linear in k.Parameter values are in arbitrary units: (top ; (bottom) . Right: the same with magnetic field instead: .
Similarly to the electric field case, the magnetic-field problem is separable and only the soluble part is plotted in figure 4.

Phosphorene
For phospherene, since the irreducible representations are all one dimensional, the Hamiltonians are the same as the dispersion relations.Thus, the band structure is simpler to interpret.From equation (13), it can be seen that (in the absence of spin-orbit interaction), only terms in even powers of k i are allowed.There is also an obvious anisotropy in the effective mass due to the orthorhombic symmetry.
Our DFT calculations (figure 5) are consistent with others already published (e.g., [3,10]).The gap is almost direct at the Γ point and the GGA value is about 0.98 eV; Tran et al [26] have shown that a GW correction increases the gap to 2.0 eV.We performed a fit to our DFT calculations (table 7).
While early calculations reported phosphorene to be a direct-gap semiconductor [3,5], subsequent calculations have revealed it to be slightly indirect [10,11,29] with an almost flat dispersion in the ΓY direction.An interesting observation is the fact that some bands merge along the surface of the Brillouin zone, leading to double degeneracy (even though the character tables of the group of wave vector only reveal one-dimensional irreducible representations).The origin of the band sticking is the nonsymmorphic nature of the space group [30].

Strain
We note that a number of DFT calculations of strained phosphorene have recently been reported [27-29, 31, 32], though all for large strains up to ∼20%.At the Γ point, the energies are where E i are the band edges in the absence of strain.Contrary to the case for silicene, now a uniaxial strain in any of the directions or a biaxial strain will change the size of the band gap if the deformation potentials for different bands are different.The particular case of a perpendicular strain has recently been modeled by compressing the  −0.19 [28] atoms within the unit cell in a DFT calculation and shown to lead to a possible semiconductor-metal transition [10].Fei et al [27] found that a biaxial strain (about 4%) or a zigzag uniaxial strain (between 5 and 6%) can change the order of the lowest two conduction bands.For a finite wave vector, the total one-dimensional Hamiltonian is

Electric field
We note that, in contrast to both graphene and silicene, equation (15) shows that the band energies change quadratically with an externally-applied electric field.At k = 0, the band energies for a perpendicular field are DFT results just published appear to confirm our prediction [28].

Magnetic field
In an external magnetic field, the one-band problem is exactly solvable in terms of the standard Landau level problem.Thus, with a static magnetic field = B B (0, 0, ) z with a vector potential = B x A (0, , 0) The eigenfunctions can be separated, where ω = eB m m c z x y * * .Thus, the prediction is that, to the lowest order, a perpendicular magnetic field will change the energy levels linearly.

Summary
The method of invariants was used to obtain analytical forms for the Hamiltonians and band structures of silicene and phosphorene beyond the lowest order terms and in the presence of external fields.Differences in the band structures of graphene, silicene and phosphorene are pointed out.Our method provides a formal analysis that is not possible from ab initio calculations.Specifically, some of the main conclusions derived from the theory include • anisotropy in the band structure of silicene shows up at the cubic order whereas they show up at the quadratic level for graphene.
• an in-plane uniaxial strain does not open a gap for silicene.
• a perpendicular strain does not open a gap for silicene but would change the gap for phosphorene.
• to the lowest-order the energies levels of silicene change linearly with an external magnetic field and a perpendicular electric field, whereas, for phosphorene, the energy levels change quadratically with a perpendicular electric field and linearly with a perpendicular magnetic field.
Ab initio calculations for the intrinsic band structures are included and more extensive calculations with external perturbations will be reported in a future publication.

Figure 1 .
Figure 1.Real (left) and reciprocal (right) space pictures of the structure of silicene.The two atoms (solid and open circles) are at different vertical heights.

Figure 2 .
Figure 2. Real (left) and reciprocal (right) space pictures of the structure of phosphorene.The two atoms (solid and open circles) are at different vertical heights.

Figure 3 .
Figure 3. Left: band dispersion for silicene along the K to Γ direction.Right: band dispersion for silicene along the K to M/2 direction.The black line is keeping only the linear term.

1 eÅ.
No existing calculations allow us to determine c 2 and c 3 .In figure 4 (left), dispersion plots are shown near the K point for an applied electric field E y = 1 (top, in arbitrary units) and for = =

Figure 5 .
Figure 5. Band structure of phosphorene.On the right is a close up of the highest valence band near the Γ point.

Table 4 .
Group of wave vector D h 2 at the Γ point (phosphorene).
The leading terms in the strain Hamiltonian are

Table 7 .
Effective masses (in units of m 0 ) of the highest valence and lowest conduction bands of phosphorene at the Γ point.The effective mass of the VB along ΓY is very sensitive to the computation.