Strain and electric field tuning of 2D hexagonal boron arsenide

Group theory and density functional theory (DFT) methods are combined to obtain compact and accurate k · p Hamiltonians that describe the bandstructures around the K and Γ points for the 2D material hexagonal boron arsenide predicted to be an important low-bandgap material for electric, thermoelectric, and piezoelectric properties that supplements the well-studied 2D material hexagonal boron nitride. Hexagonal boron arsenide is a direct bandgap material with band extrema at the K point. The bandgap becomes indirect with a conduction band minimum at the Γ point subject to a strong electric field or biaxial strain. At even higher electric field strengths (approximately 0.75 V Å−1) or a large strain (14%) 2D hexagonal boron arsenide becomes metallic. Our k · p models include to leading orders the influence of strain, electric, and magnetic fields. Excellent qualitative and quantitative agreement between DFT and k · p predictions are demonstrated for different types of strain and electric fields.


Introduction
Two-dimensional (2D) materials have opened a new field of research initiated by the micromechanical cleavage of graphene and the characterization of the unique properties of 2D materials [1,2]. Since the discovery and isolation of graphene in 2004, a vast number of 2D materials has been synthesized using various fabrication techniques [3][4][5]. With the increasing number of materials and possibilities of functionalization and combining them in heterostructures [6], they constitute a versatile class of materials with a large range of possible applications within nanotechnology. Strain engineering is a well established tool to optimize the electronic properties of semiconductors [7], and is particularly relevant for 2D materials due to their ability to sustain much larger strain magnitudes than bulk crystals [8,9]. Strain magnitudes above 10% can be achieved in graphene without damaging the material [10] and the effect of strain on the bandstructure has been studied theoretically [11,12] and experimentally [13,14]. This gives potential for use in strain-sensing applications [15][16][17].
The k·p method has been used to derive effective models to describe the electronic structure of graphene, silicene, phospherene, MoS 2 etc [18][19][20][21]. Strain is known to strongly affect the electronic properties of 2D transition metal dichalcogenides (TMDC) [22][23][24][25][26], and a recent study shows that for photodetector applications strain can tune the photoresponsivity of MoS 2 by 2-3 orders of magnitude [27]. Another way to engineer the electronic properties of 2D materials is by application of a perpendicular electric field [28,29].
Computational studies using density functional theory (DFT) reveal that boron arsenide (BAs) can form a stable 2D hexagonal structure similar to hexagonal boron nitride showing a semiconducting nature with a bandgap of approximately 1 eV [30,31]. Boron arsenide forms a cubic 3D structure with a remarkably high thermal conductivity [32,33], a property the 2D structure hexagonal boron arsenide (h-BAs) is expected to share [34,35] with application perspectives as coolant in nanodevices. In [30] they explore the effect of biaxial strain on the electronic structure and find a transition to a metallic state at a biaxial strain of 14%. The application of 2D h-BAs for gas sensing has also been examined [36]. Functionalization of h-BAs has been studied theoretically [37,38] and is found promising for applications in e.g. spintronics. Hexagonal BAs is an interesting 2D candidate since it has the same symmetry (point group D 3h ) as the important and well-studied 2D material h-BN. In contrast to h-BN, which is a high-bandgap semiconductor (5.9 eV), the bandgap of h-BAs is much smaller yet widely tunable through application of strain and external fields.
In this work, we develop analytical and computationally fast models for the band structure of 2D h-BAs based on the k·p method and group theory. Standard DFT are used to determine the relevant model parameters for Hamiltonians describing the near-bandgap states at the K and G points. We carry out a detailed investigation of strain extending the study beyond biaxial strain. An electric field perpendicular to the h-BAs plane is found to have a strong influence on the lowest conduction band that leads to a semiconductor to metal transition at high electric fields.

Effective models
In this section we will derive effective models describing the band structure close to the K and G points, including effects of external electric or magnetic fields and strain. The derivation is based on the k·p method and the Method of Invariants [11,18,[39][40][41][42][43][44] and relies on the symmetry of the crystal as well as data from DFT calculations about the atomic orbitals contributing to the relevant bands. In figure 1 we show the band structure from DFT, the crystal structure and first Brillouin zone of 2D hexagonal boron arsenide. The point group is D 3h , the character table is shown in table 1. In table 2 product tables of irreducible representations are given.
2.1. Single-particle Hamiltonian in the · k p formulation The single-particle Hamiltonian is where  · e r, g 0 μ B σ·B, and H strain are Stark, Zeeman, and strain terms, respectively. The strain Hamiltonian is unspecified but we include any term allowed by symmetry. We note in passing that Landau levels (orbital

K point
From DFT calculations we know that the fundamental band gap occurs at the K point, and the highest valence and lowest conduction bands both consist of a single p z orbital doubly degenerate due to spin. The group of the wave vector at K is C 3h [46], which has only one-dimensional representations as shown in table 1. We define q as the wave vector relative to the K point q=k−K, and q ± =q x ±iq y transform according to Γ 2 and Γ 3 respectively. The p z orbitals transform according to the Γ 4 representation and since Γ 4 ⊗Γ 4 =Γ 1 there is no direct k·p coupling between the valence and conduction bands. We consider perturbation theory to fourth order in q. Only invariant terms in q are allowed for both diagonal and off-diagonal terms since valence and conduction bands have identical representations. The allowed terms up to fourth order in q are 4 .
However, the combination of reflection symmetry σ v in the yz plane and time-reversal (TR) symmetry limits the possible third-order diagonal terms. Since we have that s , , x y x y , the dispersion around K must be even in q y . The only allowed Table 1. Character tables. w = p e i 6 . From [45]. Table 2. Product tables of irreducible representations. From [45]. (a) The k·p Hamiltonian in the basis ñ ¢ñ | | p p , z z is given by: Only second order terms are included in the off-diagonal entries, since higher-order terms do not contribute up to fourth order in the dispersion. Diagonalizing this Hamiltonian and expanding to fourth order in q gives the following dispersions for the valence and conduction bands: . Now we consider the effect of spin-orbit coupling. The SOC Hamiltonian is given by: The operators N ± =N x ±iN y transform according to the representations Γ 5 and Γ 6 , respectively, and we have no nonzero matrix elements between the valence and conduction band. N z is an invariant leading to the following expression for the spin-orbit Hamiltonian (in the basis ñ ¢ñ ñ ¢ñ However, from our DFT calculations we find that the effect is negligible. It is straightforward to extend our analysis to include strain and external electric and magnetic fields. To linear order in strain we have Including only linear strain terms we get the Hamiltonian ( ñ ¢ñ strain for q=0 and expanding to lowest order in  || gives the energies at the K point: For the electric field,  , we have the invariants  z 2 and Again we expand the energies to lowest order in the electric field giving: Finally, in the presence of a magnetic field we have both a Zeeman and a Landau contribution We now consider the valence band at the G point. DFT calculations reveal that the valence band is spanned by a p ± =p x ±ip y pair of orbitals. The group of the wave vector at the G point is D 3h , and the p ± orbitals transform according to the Γ 6 representation, and are therefore degenerate at the G point without spin-orbit coupling. First we derive the k·p Hamiltonian in the absence of spin-orbit coupling. There are no direct k·p terms between the two bands since: To get the correct dispersion and effective masses of the valence bands at the G point, second-order Löwdin perturbation theory including interactions with the s bands must be considered. Along the diagonal we get: ): giving the dispersion: where each subband is doubly degenerate due to spin. Now, consider the effect of spin-orbit coupling. The operators N ± =N x ±iN y transform according to the representation Γ 5 , while N z transforms according to Γ 2 . Since Γ 6 ⊗Γ 5 ⊗Γ 6 =Γ 5 ⊕Γ 5 ⊕Γ 3 ⊕Γ 4 ⊕Γ 5 does not contain the invariant representation there are no non-zero matrix elements between the p ± states and N ± . But the diagonal matrix elements are allowed since N z belongs to Γ 2 and Γ 6 ⊗Γ 2 ⊗Γ 6 =Γ 1 ⊕Γ 2 ⊕Γ 6 contains the invariant representation: Including the SOC Hamiltonian using the basis ñ ñ ñ ñ we find that the degeneracy is broken from 4 to 2 at G, but we still have two doubly degenerate bands, with the dispersion: The effect of strain on the valence band is given by the Hamiltonian ( ñ ñ + - Here the diagonal terms are equal because ñ  |p are two basis functions for a 2D representation in contrast to ñ |p z and ¢ñ |p z having identical, but different representations. Diagonalizing the total Hamiltonian H Γ v and H Γ v strain we get the eigenenergies: showing that for biaxial strain   = xx yy the two states remain degenerate but for shear strain or non-biaxial normal strain, the degeneracy is broken (without incorporating SOC).
For the electric field we get the Hamiltonian ( ñ ñ + - Diagonalizing gives the energies at the G point (without SOC): We see that an in-plane electric field can split the ñ  |p bands, whereas a perpendicular electric field shifts the two bands equally. The Zeeman effect is described by the following Hamiltonian (in the basis states ñ ñ ñ ñ Due to σ v reflection symmetry in the yz plane, the Landau effect vanishes between the states ñ ñ + - . DFT calculations reveal that the conduction band at the G point comprises two s bands ( ñ ¢ñ | | s s , ) described by the Hamiltonian = + + where iä{1, 2} denotes the two subbands. Since N ± transform according to Γ 5 and N z according to Γ 2 there are no spin-orbit interactions between the two conduction band states. The effect of strain on the conduction band is given by the Hamiltonian ( ñ ¢ñ | | s s , ) Here we have included a second-order term proportional to    xx yy xy 2 2 , since there is no firstorder term and from DFT we see an effect of this within reasonable strain magnitudes. The lowest-order contributions to the energies at the G point are then given by: the valence band a p z orbital at the K point. At the G point, the lowest conduction bands are p z and s states while the valence band is a pair of p ± orbitals. Hence, the present Hamiltonian models for the valence bands at the G and K points can be directly applied to hexagonal boron nitride as well.

Computational method and extraction of parameters
We determine parameters of the models derived in the previous section by fitting to ab initio calculations. We have performed standard DFT calculations using the projector augmented wave method as implemented in the GPAW package [47,48] with the Perdew-Burke-Ernzerhof approximation of the exchange-correlation energy [49]. The selfconsistent calculation used a 16x16 Monkhorst-Pack k point grid and a 1200 eV cut-off energy for the plane-wave representation of the wave functions. The structure was relaxed using the Broyden-Fletcher-Goldfarb-Shanno algorithm as implemented in ASE [50], resulting in a lattice constant of 3.39 Å. We use a vacuum layer of 20 Å between adjacent layers to avoid interaction between supercells. The SOC Hamiltonian was diagonalized in a basis of wave functions from scalar-relativistic calculations as documented in [51]. All fitted parameters are listed in table 3.
To extract the k·p parameters we perform bandstructure calculations on a dense k point grid in the the vicinity of either K or G. In figure 2, a contour plot of the valence and conduction bands at the K point is shown using both DFT and k·p results from equation (4b) with parameters fitted to the DFT energies. The model captures the trigonal warping of the bands described by the third-order terms in q. We clearly see that the bands are symmetric under q q y y due to the combination of TR and the σ v reflection as explained in the previous section. The k·p model gives very good agreement with DFT energies up to a distance of 0.1 Å −1 from the K point, due to the inclusion of fourth-order terms in q.
At the G point, k·p parameters are obtained by fitting equation (17) to the DFT result without SOC, and the SOC parameter is simply the energy splitting at the G point. The dispersions given by equation (20) are in good agreement with the DFT results within 0.1 Å −1 of G as seen in figure 3(b). For the conduction band only the diagonal terms in equation (28) contribute to second order in the dispersion, and we find that keeping only those terms gives a relatively good description of the bands as shown in figure 3(a).

Electric field
Using the external potential module implemented in GPAW we perform calculations of the band structure with a constant electric field in the z direction. In accordance with the symmetry analysis in the previous section we see in figure 4(a) a parabolic dependence of the eigenvalues on the strength of the electric field. Interestingly, most bands are only slightly changed, but the lowest conduction band at the G point shows high sensitivity to electric fields. At large electric fields the conduction band minimum changes from the K point to the G point at which point BAs becomes an indirect bandgap material. At electric fields larger than 0.78 V Å −1 , the conduction band decreases below the valence band maximum at K and we get a transition to a metallic state. In figure 4(b) the bandstructure of h-BAs in the metallic state is shown. It should be pointed out that an electric field strength of 0.78 V Å −1 , required for h-BAs to become metallic, may induce instability of the structure. Yet, we point to that fields of the same order of magnitude have been demonstrated experimentally for bilayer graphene [52]. We fit equations (31a)-(31b) to the DFT results to obtain the relevant parameters.The second conduction band at G is not fitted since higher-lying bands are shifted down and mixed with this band.

Strain
To investigate the effect of strain in the linear approximation we modify the unit cell after the relaxation,  ¢ = + a a a n n n where ò is the strain tensor, a n are the equilibrium lattice vectors, and ¢ a n are the strained lattice vectors. After straining the unit cell we perform a relaxation of the atomic positions while fixing the unit cell fixed to the strained value. There are three independent terms in the in-plane strain, biaxial strain with ò xx =ò yy , ò xx =−ò yy , and shear strain ò xy . Only in the biaxial case is the symmetry of the lattice preserved, and we note that for the band structure calculation the Brillouin zone and the high symmetry points are changed. Hence, we consider the point in k space that reduces to K when strain goes to zero and we simply denote it K. We limit our study to the cases where only one of the strain terms is nonzero. At the K point we find that the only linear strain terms that occur in the Hamiltonian, equation (7), are proportional to  || . By fitting equation (8b) to the DFT results we determine the parameters D 1 and D 3 and obtain excellent agreement within a range of  =  || 10% as seen in figure 5(a). In figure 5(b) we observe that the eigenvalues at the K point remain constant with both Figure 3. Dispersion around the G point from DFT (blue solid) and k·p (black dashed) in the k x direction including spin-orbit coupling. The coefficients of the k·p terms are obtained from a fit to the DFT result without spin-orbit coupling and the spin-orbit constant is simply obtained from the eigenvalues at G. ò xx −ò yy and ò xy , showing that the effects of strain is well described by the linear strain terms For the valence band at the G point we fit equation (22) to the DFT results. The resuls are plotted in figure 5, showing good agreement between model and DFT calculations. Close inspection shows that higher-order terms contribute giving a slight curvature of the bands, but the energy shifts are small compared to the contribution from the linear terms The higher-order terms also give a difference between ò xx −ò yy and shear strain ò xy even though they belong to the same representation in the symmetry analysis. This comes from the fact that the shear term is multiplied by the imaginary unit i. The conduction bands are well described by the linear term in  || as well, however, for ò xx −ò yy and ò xy there are no linear terms, but we find a significant parabolic dependence in figure 5(b). Again, higher-order terms show little effect for large strains yet for most cases the second-order term gives a good approximation.

Conclusion
Compact and accurate k·p Hamiltonian models for the 2D material hexagonal boron arsenide (h-BAs) are extracted based on group symmetry discussions. Model coefficients are found by comparison with detailed DFT calculations. 2D hexagonal boron arsenide is expected to be an important material as it is predicted to have an ultrahigh thermal conductivity and provides a low-bandgap candidate to the well-studied high-bandgap 2D semiconductor material h-BN. We derive k · p Hamiltonians and demonstrate excellent agreement with DFT results in the presence of strain and electric fields at two important points in the Brillouin zone, K and G. DFT calculations reveal that the bandgap of 2D boron arsenide, located at the K point, becomes indirect at sufficiently large electric fields with the conduction band minimum located at the G point. At even larger electric fields or a large strain the material becomes metallic. The influence of an external magnetic field is also discussed. The k·p models derived in this work are useful for practical and efficient device simulations in addition to heterostructures, complicated geometries, or accounting for spatially varying external fields where the periodicity is broken and DFT becomes computationally too demanding.