Electronic structure and tunability of 2D hexagonal boron arsenide

TL;DR

This paper combines group theory and density functional theory to develop accurate models of the electronic bandstructure of 2D hexagonal boron arsenide, revealing its tunable electronic properties under strain and electric fields.

Contribution

It introduces a new $k extit{ extperiodcentered}p$ Hamiltonian model for h-BAs, capturing effects of strain, electric, and magnetic fields with high accuracy.

Findings

h-BAs is a direct bandgap material at the K point

Electric field or biaxial strain can make the bandgap indirect

High electric field or strain induces metallic behavior in h-BAs

Abstract

Group theory and density functional theory methods are combined to obtain compact and accurate $k\cdot p$ Hamiltonians that describe the bandstructures around the $K$ and $\Gamma$ points for the 2D material hexagonal boron arsenide (h-BAs) 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 $\Gamma$ point subject to a strong electric field or biaxial strain. At even higher electric field strengths (approximately 0.75 V/$\r{A}$) or a large strain ($14$~\%) 2D hexagonal boron arsenide becomes metallic. Our $k\cdot p$ models include to leading orders the influence of strain, electric, and magnetic fields. Excellent qualitative and quantitative agreement between density functional theory and $k\cdot p$ predictions are demonstrated for different types of strain and electric fields.

Tables

Table 1. (a)

$D_{3h}$	$E$	${\sigma }_h$	$2C_3$	$2S_3$	$3C_2$	$3{\sigma }_v$
${\Gamma }_1$	1	1	1	1	1	1	$q_{\left|\right|}^2,\mathrm{Re}\left[q_{+}^3\right]$
${\Gamma }_2$	1	1	1	1	$-1$	$-1$	$\mathrm{Im}\left[q_{+}^3\right]$
${\Gamma }_3$	1	$-1$	1	$-1$	1	$-1$
${\Gamma }_4$	1	$-1$	1	$-1$	$-1$	1
${\Gamma }_5$	2	$-2$	$-1$	1	0	0
${\Gamma }_6$	2	2	$-1$	$-1$	0	0	$\{q_{+},q_{-}\},\{q_{-}^2,q_{+}^2\}$

Table 2. (a)

$D_{3h}$	$E$	${\sigma }_h$	$2C_3$	$2S_3$	$3C_2$	$3{\sigma }_v$
${\Gamma }_1$	1	1	1	1	1	1	$q_{\left|\right|}^2,\mathrm{Re}\left[q_{+}^3\right]$
${\Gamma }_2$	1	1	1	1	$-1$	$-1$	$\mathrm{Im}\left[q_{+}^3\right]$
${\Gamma }_3$	1	$-1$	1	$-1$	1	$-1$
${\Gamma }_4$	1	$-1$	1	$-1$	$-1$	1
${\Gamma }_5$	2	$-2$	$-1$	1	0	0
${\Gamma }_6$	2	2	$-1$	$-1$	0	0	$\{q_{+},q_{-}\},\{q_{-}^2,q_{+}^2\}$

Table 3. (b)

$C_{3h}$	$E$	$S_3^{-1}$	$C_3$	${\sigma }_h$	$C_3^{-1}$	$S_3$
${\Gamma }_1$	1	1	1	1	1	1	$q_{\left|\right|}^2,q_{+}^3,q_{-}^3$
${\Gamma }_2$	1	$-{\omega }^2$	${\omega }^4$	1	$-{\omega }^2$	${\omega }^4$	$q_{+},q_{-}^2$
${\Gamma }_3$	1	${\omega }^4$	$-{\omega }^2$	1	${\omega }^4$	$-{\omega }^2$	$q_{-},q_{+}^2$
${\Gamma }_4$	1	$-1$	1	$-1$	1	$-1$
${\Gamma }_5$	1	${\omega }^2$	${\omega }^4$	$-1$	$-{\omega }^2$	$-{\omega }^4$
${\Gamma }_6$	1	$-{\omega }^4$	$-{\omega }^2$	$-1$	${\omega }^4$	${\omega }^2$

Table 4. (a)

$D_{3h}$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_1$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_2$	${\Gamma }_2$	${\Gamma }_1$	${\Gamma }_4$	${\Gamma }_3$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_3$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_6$	${\Gamma }_5$
${\Gamma }_4$	${\Gamma }_4$	${\Gamma }_3$	${\Gamma }_2$	${\Gamma }_1$	${\Gamma }_6$	${\Gamma }_5$
${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_1\oplus {\Gamma }_2\oplus {\Gamma }_6$	${\Gamma }_3\oplus {\Gamma }_4\oplus {\Gamma }_5$
${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_3\oplus {\Gamma }_4\oplus {\Gamma }_5$	${\Gamma }_1\oplus {\Gamma }_2\oplus {\Gamma }_6$

Table 5. (a)

$D_{3h}$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_1$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_2$	${\Gamma }_2$	${\Gamma }_1$	${\Gamma }_4$	${\Gamma }_3$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_3$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_6$	${\Gamma }_5$
${\Gamma }_4$	${\Gamma }_4$	${\Gamma }_3$	${\Gamma }_2$	${\Gamma }_1$	${\Gamma }_6$	${\Gamma }_5$
${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_1\oplus {\Gamma }_2\oplus {\Gamma }_6$	${\Gamma }_3\oplus {\Gamma }_4\oplus {\Gamma }_5$
${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_3\oplus {\Gamma }_4\oplus {\Gamma }_5$	${\Gamma }_1\oplus {\Gamma }_2\oplus {\Gamma }_6$

Table 6. (b)

$C_{3h}$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_1$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$
${\Gamma }_2$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_1$	${\Gamma }_5$	${\Gamma }_6$	${\Gamma }_4$
${\Gamma }_3$	${\Gamma }_3$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_6$	${\Gamma }_4$	${\Gamma }_5$
${\Gamma }_4$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_6$	${\Gamma }_1$	${\Gamma }_2$	${\Gamma }_3$
${\Gamma }_5$	${\Gamma }_5$	${\Gamma }_6$	${\Gamma }_4$	${\Gamma }_2$	${\Gamma }_3$	${\Gamma }_1$
${\Gamma }_6$	${\Gamma }_6$	${\Gamma }_4$	${\Gamma }_5$	${\Gamma }_3$	${\Gamma }_1$	${\Gamma }_2$

Table 7. (a)

$E^{𝐊c}$ $\left[\mathrm{eV}\right]$	$-2.44$
$a$ $\left[\mathrm{eV}{\AA }^2\right]$	$21.4$
$b$ $\left[\mathrm{eV}{\AA }^3\right]$	$-25.0$
$\overline{c}$ $\left[\mathrm{eV}{\AA }^4\right]$	$-245$
$E^{𝐊v}$ $\left[\mathrm{eV}\right]$	$-3.20$
$a^{\prime }$ $\left[\mathrm{eV}{\AA }^2\right]$	$-22.2$
$b^{\prime }$ $\left[\mathrm{eV}{\AA }^3\right]$	$20.0$
${\overline{c}}^{\prime }$ $\left[\mathrm{eV}{\AA }^4\right]$	$250$
$E^{𝚪v}$ $\left[\mathrm{eV}\right]$	$-4.75$
$a^{𝚪v}$ $\left[\mathrm{eV}{\AA }^2\right]$	$-13.9$
$c^{𝚪v}$ $\left[\mathrm{eV}{\AA }^2\right]$	$7.02$
$\Delta$ $\left[\mathrm{eV}\right]$	$0.074$
$E_1^{𝚪c}$ $\left[\mathrm{eV}\right]$	$0.57$
$E_2^{𝚪c}$ $\left[\mathrm{eV}\right]$	$0.40$
$a_1^{𝚪c}$ $\left[\mathrm{eV}{\AA }^2\right]$	$8.25$
$a_2^{𝚪c}$ $\left[\mathrm{eV}{\AA }^2\right]$	$6.39$

Table 8. (a)

$E^{𝐊c}$ $\left[\mathrm{eV}\right]$	$-2.44$
$a$ $\left[\mathrm{eV}{\AA }^2\right]$	$21.4$
$b$ $\left[\mathrm{eV}{\AA }^3\right]$	$-25.0$
$\overline{c}$ $\left[\mathrm{eV}{\AA }^4\right]$	$-245$
$E^{𝐊v}$ $\left[\mathrm{eV}\right]$	$-3.20$
$a^{\prime }$ $\left[\mathrm{eV}{\AA }^2\right]$	$-22.2$
$b^{\prime }$ $\left[\mathrm{eV}{\AA }^3\right]$	$20.0$
${\overline{c}}^{\prime }$ $\left[\mathrm{eV}{\AA }^4\right]$	$250$
$E^{𝚪v}$ $\left[\mathrm{eV}\right]$	$-4.75$
$a^{𝚪v}$ $\left[\mathrm{eV}{\AA }^2\right]$	$-13.9$
$c^{𝚪v}$ $\left[\mathrm{eV}{\AA }^2\right]$	$7.02$
$\Delta$ $\left[\mathrm{eV}\right]$	$0.074$
$E_1^{𝚪c}$ $\left[\mathrm{eV}\right]$	$0.57$
$E_2^{𝚪c}$ $\left[\mathrm{eV}\right]$	$0.40$
$a_1^{𝚪c}$ $\left[\mathrm{eV}{\AA }^2\right]$	$8.25$
$a_2^{𝚪c}$ $\left[\mathrm{eV}{\AA }^2\right]$	$6.39$

Table 9. (b)

$D_1$ $\left[\mathrm{eV}\right]$	$-2.48$
$D_3$ $\left[\mathrm{eV}\right]$	$-3.50$
$D_1^{𝚪v}$ $\left[\mathrm{eV}\right]$	$-4.61$
$D_2^{𝚪v}$ $\left[\mathrm{eV}\right]$	$-5.06$
$D_1^{𝚪c}$ $\left[\mathrm{eV}\right]$	$-0.517$
$D_2^{𝚪c}$ $\left[\mathrm{eV}\right]$	$-7.47$
$D_5^{𝚪c}$ $\left[\mathrm{eV}\right]$	$-19.3$
$D_6^{𝚪c}$ $\left[\mathrm{eV}\right]$	$-64.5$

Table 10. (c)

$A_1\left[\mathrm{eV}{\AA }^2/V^2\right]$	$-0.259$
$A_5\left[\mathrm{eV}{\AA }^2/V^2\right]$	$-0.271$
$A_1^{𝚪v}\left[\mathrm{eV}{\AA }^2/V^2\right]$	$-0.285$
$A_5^{𝚪c}\left[\mathrm{eV}{\AA }^2/V^2\right]$	$-6.91$

Equations

$$H$$

$$+g_0{\mu }_B\sigma \cdot \mathbf{B}+H_{\text{strain}},$$

$$\mathbf{p}\to \mathbf{p}+e\mathbf{A}.$$

$$H^{\mathbf{K}}=$$

$$\left(\begin{array}{cc}{E}^{\mathbf{K}c}+a{q}_{||}^2+b\mathrm{Re}[q_{+}^3]+c{q}_{||}^4& d{q}_{||}^2\\ d^{\ast }{q}_{||}^2& E^{\mathbf{K}v}+a^{\prime }{q}_{||}^2+b^{\prime }\mathrm{Re}[q_{+}^3]+c^{\prime }{q}_{||}^4\end{array}\right).$$

$$E_c(\mathbf{q})$$

$$E_v(\mathbf{q})$$

$$H_{\text{soc}}=\frac{\hslash }{4m_0^2{c}^2}(\nabla V\times \mathbf{p})\cdot \sigma =\frac{\hslash }{4m_0^2{c}^2}\mathbf{N}\cdot \sigma .$$

$$H_{\text{SOC}}^{\mathbf{K}}=\left(\begin{array}{cccc}{\Delta }_1& {\Delta }_3& 0& 0\\ {\Delta }_3^{\ast }& {\Delta }_2& 0& 0\\ 0& 0& -{\Delta }_1& -{\Delta }_3\\ 0& 0& -{\Delta }_3^{\ast }& -{\Delta }_2\end{array}\right).$$

$$H_{\text{strain}}^{\mathbf{K}}=\left(\begin{array}{cc}{D}_1{ϵ}_{||}& D_2{ϵ}_{||}\\ D_2^{\ast }{ϵ}_{||}& D_3{ϵ}_{||}.\end{array}\right).$$

$$E^{\mathbf{K}c}(\mathbf{ϵ})=E^{\mathbf{K}c}+D_1{ϵ}_{||}$$

$$E^{\mathbf{K}v}(\mathbf{ϵ})=E^{\mathbf{K}v}+D_3{ϵ}_{||}.$$

$$H_{\text{E-field}}^{\mathbf{K}}=\left(\begin{array}{cc}{A}_1{\mathcal{E}}_z^2+A_2{\mathcal{E}}_{||}^2& A_3{\mathcal{E}}_z^2+A_4{\mathcal{E}}_{||}^2\\ A_3^{\ast }{\mathcal{E}}_z^2+A_4^{\ast }{\mathcal{E}}_{||}^2& A_5{\mathcal{E}}_z^2+A_6{\mathcal{E}}_{||}^2\end{array}\right).$$

$$E^{\mathbf{K}c}(\mathcal{E})=E^{\mathbf{K}c}+A_1{\mathcal{E}}_z^2+A_2{\mathcal{E}}_{||}^2$$

$$E^{\mathbf{K}v}(\mathcal{E})=E^{\mathbf{K}v}+A_5{\mathcal{E}}_z^2+A_6{\mathcal{E}}_{||}^2.$$

$$H_{\text{Zeeman}}^{\mathbf{K}}=\left(\begin{array}{cccc}{g}_1{\mu }_B{B}_z& 0& g_1{\mu }_B{B}_{-}& 0\\ 0& g_2{\mu }_B{B}_z& 0& g_2{\mu }_B{B}_{-}\\ g_1{\mu }_B{B}_{+}& 0& -g_1{\mu }_B{B}_z& 0\\ 0& g_2{\mu }_B{B}_{+}& 0& -g_2{\mu }_B{B}_z\end{array}\right),$$

$$H_{\text{Landau}}^{\mathbf{K}}=\left(\begin{array}{cccc}{\alpha }_1{B}_z& {\alpha }_2{B}_z& 0& 0\\ {\alpha }_2^{\ast }{B}_z& {\alpha }_3{B}_z& 0& 0\\ 0& 0& {\alpha }_1{B}_z& {\alpha }_2{B}_z\\ 0& 0& {\alpha }_2^{\ast }{B}_z& {\alpha }_3{B}_z\end{array}\right).$$

$$⟨p_{+}|\mathbf{p}|p_{-}⟩=⟨p_{+}|\frac{m_0}{i\hslash }[H_0,\mathbf{r}]|p_{-}⟩=0,$$

$$\frac{{\hslash }^2{k}_{||}^2}{2m_0}+\frac{{\hslash }^2}{m_0^2}\sum _s\frac{⟨p_{\pm }|\mathbf{k}\cdot \mathbf{p}|s⟩⟨s|\mathbf{k}\cdot \mathbf{p}|p_{\pm }⟩}{E_s-E_p}=a^{\mathbf{\Gamma }v}{k}_{||}^2,$$

$$\frac{{\hslash }^2}{m_0^2}\sum _s\frac{⟨p_{+}|\mathbf{k}\cdot \mathbf{p}|s⟩⟨s|\mathbf{k}\cdot \mathbf{p}|p_{-}⟩}{E_s-E_p}=c^{\mathbf{\Gamma }v}{k}_{-}^2.$$

$$H^{\mathbf{\Gamma }v}=E^{\mathbf{\Gamma }v}+\left(\begin{array}{cc}{a}^{\mathbf{\Gamma }v}{k}_{||}^2& c^{\mathbf{\Gamma }v}{k}_{-}^2\\ c^{\mathbf{\Gamma }v\ast }{k}_{+}^2& a^{\mathbf{\Gamma }v}{k}_{||}^2\end{array}\right),$$

$$E_{\pm }^{\mathbf{\Gamma }v}(\mathbf{k})=E^{\mathbf{\Gamma }v}+(a^{\mathbf{\Gamma }v}\pm |c^{\mathbf{\Gamma }v}|)k_{||}^2,$$

$$⟨p_{+}|N_z|p_{+}⟩=-⟨p_{-}|N_z|p_{-}⟩=\Delta .$$

$$H_{SOC}^{\mathbf{\Gamma }v}$$

$$E_{\pm }^{\mathbf{\Gamma }v}(\mathbf{k})=E^{\mathbf{\Gamma }v}+a^{\mathbf{\Gamma }v}{k}_{||}^2\pm \sqrt{{\Delta }^2+|c^{\mathbf{\Gamma }v}{|}^2{k}_{||}^4}.$$

$$H_{\text{strain}}^{\mathbf{\Gamma }v}=\left(\begin{array}{cc}{D}_1^{\mathbf{\Gamma }v}{ϵ}_{||}& D_2^{\mathbf{\Gamma }v}{ϵ}_{+}\\ {D_2^{\mathbf{\Gamma }v}}^{\ast }{ϵ}_{-}& D_1^{\mathbf{\Gamma }v}{ϵ}_{||}\end{array}\right).$$

$$E_{\pm }^{\mathbf{\Gamma }v}(\mathbf{ϵ})=E^{\mathbf{\Gamma }v}+D_1^{\mathbf{\Gamma }v}{ϵ}_{||}\pm |D_2^{\mathbf{\Gamma }v}|\sqrt{({ϵ}_{xx}-{ϵ}_{yy}{)}^2+4{ϵ}_{xy}^2},$$

$$H_{\text{E-field}}^{\mathbf{\Gamma }v}=\left(\begin{array}{cc}{A}_1^{\mathbf{\Gamma }v}{\mathcal{E}}_z^2+A_2^{\mathbf{\Gamma }v}{\mathcal{E}}_{||}^2& A_3^{\mathbf{\Gamma }v}{\mathcal{E}}_{+}^2+A_4^{\mathbf{\Gamma }v}{\mathcal{E}}_{-}{\mathcal{E}}_z\\ A_3^{\mathbf{\Gamma }v\ast }{\mathcal{E}}_{-}^2+A_4^{\mathbf{\Gamma }v\ast }{\mathcal{E}}_{+}{\mathcal{E}}_z& A_1^{\mathbf{\Gamma }v}{\mathcal{E}}_z^2+A_2^{\mathbf{\Gamma }v}{\mathcal{E}}_{||}^2\end{array}\right).$$

$$E_{\pm }^{\mathbf{\Gamma }v}(\mathcal{E})=E^{\mathbf{\Gamma }v}+A_1^{\mathbf{\Gamma }v}{\mathcal{E}}_z^2+A_2^{\mathbf{\Gamma }v}{\mathcal{E}}_{||}^2\pm \left|A_3^{\mathbf{\Gamma }v}{\mathcal{E}}_{+}^2+A_4^{\mathbf{\Gamma }v}{\mathcal{E}}_{-}{\mathcal{E}}_z\right|.$$

$$H_{\text{Zeeman}}^{\mathbf{\Gamma }v}=\left(\begin{array}{cccc}{g}_1^{\mathbf{\Gamma }v}{\mu }_B{B}_z& 0& g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_{-}& 0\\ 0& g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_z& 0& g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_{-}\\ g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_{+}& 0& -g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_z& 0\\ 0& g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_{+}& 0& -g_1^{\mathbf{\Gamma }v}{\mu }_B{B}_z\end{array}\right).$$

$$H^{\mathbf{\Gamma }c}=\left(\begin{array}{cc}{E}_1^{\mathbf{\Gamma }c}+a_1^{\mathbf{\Gamma }c}{k}_{||}^2& c^{\mathbf{\Gamma }c}{k}_{||}^2\\ c^{\mathbf{\Gamma }c\ast }{k}_{||}^2& E_2^{\mathbf{\Gamma }c}+a_2^{\mathbf{\Gamma }c}{k}_{||}^2.\end{array}\right)$$

$$E_i^{\mathbf{\Gamma }c}=E_i^{\mathbf{\Gamma }c}+a_i^{\mathbf{\Gamma }c}{k}_{||}^2,$$

$$H_{\text{strain}}^{\mathbf{\Gamma }c}=\left(\begin{array}{cc}{D}_1^{\mathbf{\Gamma }c}{ϵ}_{||}+D_2^{\mathbf{\Gamma }c}{ϵ}_{+}{ϵ}_{-}& D_3^{\mathbf{\Gamma }c}{ϵ}_{||}+D_4^{\mathbf{\Gamma }c}{ϵ}_{+}{ϵ}_{-}\\ D_3^{\mathbf{\Gamma }c\ast }{ϵ}_{||}+D_4^{\mathbf{\Gamma }c\ast }{ϵ}_{+}{ϵ}_{-}& D_5^{\mathbf{\Gamma }c}{ϵ}_{||}+D_6^{\mathbf{\Gamma }c}{ϵ}_{+}{ϵ}_{-}\end{array}\right).$$

Full text

I 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 materialsNovoselov et al. (2004, 2005). Since the discovery and isolation of graphene in 2004, a vast number of 2D materials has been synthesized using various fabrication techniques Hao et al. (2018); Li et al. (2018); Shen et al. (2018). With the increasing number of materials and possibilities of functionalization and combining them in heterostructuresLiu et al. (2019), 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 semiconductorsBedell et al. (2014), and is particularly relevant for 2D materials due to their ability to sustain much larger strain magnitudes than bulk crystalsDai et al. ; Akinwande et al. (2017). Strain magnitudes above $10\%$ can be achieved in graphene without damaging the materialPérez Garza et al. (2014) and the effect of strain on the bandstructure has been studied theoreticallyWinkler and Zülicke (2010); Cocco et al. (2010) and experimentallyYan et al. (2012); Li et al. (2015). This gives potential for use in strain-sensing applicationsLee et al. (2010); Wang et al. (2015); Yang et al. (2018).

The $k\cdot p$ method has been used to derive effective models to describe the electronic structure of graphene, silicene, phospherene, MoS2 etc. Voon et al. (2015); Xiao et al. (2012); Beiranvand et al. (2018); Kormányos et al. (2013). Strain is known to strongly affect the electronic properties of two-dimensional transition metal dichalcogenides (TMDC) He et al. (2016); Schmidt et al. (2016); Conley et al. (2013); Wu et al. (2018); Castellanos-Gomez et al. (2013), and a recent study shows that for photodetector applications strain can tune the photoresponsivity of MoS2 by 2-3 orders of magnitudeGant et al. (2019). Another way to engineer the electronic properties of 2D materials is by application of a perpendicular electric fieldZubair et al. (2017); Qi et al. (2013).

Computational studies revealed that boron arsenide (BAs) can form a stable 2D hexagonal structure similar to hexagonal boron nitride and explored the material properties using density functional theory, showing a semiconducting nature with a bandgap of approximately $1\text{\,}\mathrm{eV}$Manoharan and Subramanian (2018); Şahin et al. (2009). Boron arsenide forms a cubic 3D structure with a remarkably high thermal conductivityTian et al. (2018); Kang et al. (2018), a property the 2D structure h-BAs is expected to shareShi and Luo (2018) with application perspectives as coolant in nanodevices. In Ref. Manoharan and Subramanian, 2018 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 examinedRen et al. (2019). Functionalization of h-BAs has been studied theoreticallywu Zhang et al. (2015); Ullah et al. (2018) 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 bandstructure of 2D h-BAs based on the $k\cdot p$ method and group theory. We aim to gain a physical understanding of this material as well as a simple accurate model describing the band structure close to the high symmetry points both qualitatively and quantitatively including the spin-orbit coupling (SOC). Standard density functional theory calculations are used to determine the relevant model parameters. We carry out a detailed investigation of strain effects and symmetry-allowed terms extending the study beyond biaxial strain. A perpendicular electric field is found to have a strong influence on the lowest conduction band, giving rise to a semiconductor to metal transition for strong electric fields.

II Effective models

In this section we will derive effective models describing the band structure close to the $\mathbf{K}$ and $\mathbf{\Gamma}$ point, including effects of external electric or magnetic fields and strain. The derivation is based on the $k\cdot p$ methodLew Yan Voon and Willatzen (2009); Bir and Pikus (1974) and relies on the symmetry of the crystal as well as information from DFT calculations about the atomic orbitals contributing to the relevant bands. In Fig. 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.

II.1 Single-particle Hamiltonian in the $k\cdot p$ formulation

The single-particle Hamiltonian is

[TABLE]

where $e{\bm{\mathcal{E}}}\cdot{\bf r}$, $g_{0}\mu_{B}{\bf\sigma}\cdot{\bf B}$, and $H_{\text{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 magnetic field effects) are invoked by minimal substitution, i.e.,

[TABLE]

II.2 $\mathbf{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 $\mathbf{K}$ is $C_{3h}$ Dresselhaus et al. (2008), which has only one-dimensional representations as shown in Table 1. We define $\mathbf{q}$ as the wave vector relative to the $\mathbf{K}$ point $\mathbf{q}=\mathbf{k}-\mathbf{K}$, and $q_{\pm}=q_{x}\pm iq_{y}$ transform according to $\Gamma_{2}$ and $\Gamma_{3}$ respectively. The $p_{z}$ orbitals transform according to the $\Gamma_{4}$ representation and since $\Gamma_{4}\otimes\Gamma_{4}=\Gamma_{1}$ there is no direct $k\cdot p$ coupling between the valence and conduction bands. We consider perturbation theory to fourth order in $\mathbf{q}$. Only invariant terms in $\mathbf{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 $\mathbf{q}$ are $q_{||}^{2},q_{+}^{3},q_{-}^{3},q_{||}^{4}$. However, the combination of reflection symmetry $\sigma_{v}$ in the $yz$ plane and time-reversal (TR) symmetry limits the possible third-order diagonal terms. Since we have that $\sigma_{v}\mathrel{\mathop{\mathchar 58\relax}}\mathbf{K}+(q_{x},q_{y})\rightarrow-\mathbf{K}-(q_{x},-q_{y})$ and $TR\mathrel{\mathop{\mathchar 58\relax}}-\mathbf{K}-(q_{x},-q_{y})\rightarrow\mathbf{K}+(q_{x},-q_{y})$, the dispersion around $\mathbf{K}$ must be even in $q_{y}$. The only allowed third-order diagonal term is therefore $\operatorname{Re}[q_{+}^{3}]=q_{x}^{3}-3q_{x}q_{y}^{2}$. The $k\cdot p$ Hamiltonian in the basis $|p_{z}\rangle,|p_{z}^{\prime}\rangle$ is given by:

[TABLE]

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 $\mathbf{q}$ gives the following dispersions for the valence and conduction bands:

[TABLE]

where $\bar{c}=\left(c+\frac{|d|^{2}}{E^{\mathbf{K}c}-E^{\mathbf{K}v}}\right)$ and $\bar{c}^{\prime}=\left(c^{\prime}+\frac{|d|^{2}}{E^{\mathbf{K}v}-E^{\mathbf{K}_{c}}}\right)$. Now we consider the effect of spin-orbit coupling. The SOC Hamiltonian is given by:

[TABLE]

The operators $N_{\pm}=N_{x}\pm iN_{y}$ transform according to the representations $\Gamma_{5}$ and $\Gamma_{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 $|p_{z}\uparrow\rangle,|p_{z}^{\prime}\uparrow\rangle,|p_{z}\downarrow\rangle,|p_{z}^{\prime}\downarrow\rangle$):

[TABLE]

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 $\epsilon_{||}=\epsilon_{xx}+\epsilon_{yy}$ which is invariant, and $\epsilon_{+}=\epsilon_{xx}-\epsilon_{yy}+2i\epsilon_{xy}$ ($\epsilon_{-}=\epsilon_{xx}-\epsilon_{yy}-2i\epsilon_{xy}$) transforming according to $\Gamma_{2}$ ($\Gamma_{3}$). Including only linear strain terms we get the Hamiltonian ($|p_{z}\rangle,|p_{z}^{\prime}\rangle$):

[TABLE]

Diagonalizing $H^{\mathbf{K}}+H^{\mathbf{K}}_{\text{strain}}$ for $\mathbf{q}=0$ and expanding to lowest order in $\epsilon_{||}$ gives the energies at the $\mathbf{K}$ point:

[TABLE]

For the electric field, $\bm{\mathcal{E}}$, we have the invariants $\mathcal{E}_{z}^{2}$ and $\mathcal{E}_{||}^{2}=\mathcal{E}_{x}^{2}+\mathcal{E}_{y}^{2}$, thus in the basis ($|p_{z}\rangle,|p_{z}^{\prime}\rangle$):

[TABLE]

Again we expand the energies to lowest order in the electric field giving:

[TABLE]

Finally, in the presence of a magnetic field we have both a Zeeman and a Landau contribution ($|p_{z}\uparrow\rangle,|p_{z}^{\prime}\uparrow\rangle,|p_{z}\downarrow\rangle,|p_{z}^{\prime}\downarrow\rangle$):

[TABLE]

[TABLE]

II.3 $\mathbf{\Gamma}$ Point

We now consider the valence band at the $\mathbf{\Gamma}$ point. DFT calculations reveal that the valence band is spanned by a $p_{\pm}=p_{x}\pm ip_{y}$ pair of orbitals. The group of the wave vector at the $\Gamma$ point is $D_{3h}$, and the $p_{\pm}$ orbitals transform according to the $\Gamma_{6}$ representation, and are therefore degenerate at the $\mathbf{\Gamma}$ point without spin-orbit coupling. First we derive the $k\cdot p$ Hamiltonian in the absence of spin-orbit coupling. There are no direct $k\cdot p$ terms between the two bands since:

[TABLE]

where we used $E_{0}(p_{+})=E_{0}(p_{-})$. Next, we consider second-order Löwdin perturbation theory including interaction terms with $s$ bands. Along the diagonal we get:

[TABLE]

defining the constant $a^{\mathbf{\Gamma}v}$. Similarly we have the off diagonal terms:

[TABLE]

Thus without spin-orbit coupling we have to second order in $\mathbf{k}$ the Hamiltonian ($|p_{+}\rangle,|p_{-}\rangle$):

[TABLE]

giving the dispersion:

[TABLE]

where each subband is doubly degenerate due to spin. Now, consider the effect of spin-orbit coupling. The operators $N_{\pm}=N_{x}\pm iN_{y}$ transform according to the representation $\Gamma_{5}$, while $N_{z}$ transforms according to $\Gamma_{2}$. Since $\Gamma_{6}\otimes\Gamma_{5}\otimes\Gamma_{6}=\Gamma_{5}\oplus\Gamma_{5}\oplus\Gamma_{3}\oplus\Gamma_{4}\oplus\Gamma_{5}$ does not contain the invariant representation there are no non-zero matrix elements between the $p_{\pm}$ states and $N_{\pm}$.. But the diagonal matrix elements are allowed since $N_{z}$ belongs to $\Gamma_{2}$ and $\Gamma_{6}\otimes\Gamma_{2}\otimes\Gamma_{6}=\Gamma_{1}\oplus\Gamma_{2}\oplus\Gamma_{6}$ contains the invariant representation:

[TABLE]

Including the SOC Hamiltonian using the basis $|p_{+}\uparrow\rangle,|p_{-}\uparrow\rangle,|p_{+}\downarrow\rangle,|p_{-}\downarrow\rangle$:

[TABLE]

we find that the degeneracy is broken from 4 to 2 at $\mathbf{\Gamma}$, but we still have two doubly degenerate bands, with the dispersion:

[TABLE]

The effect of strain on the valence band is given by the Hamiltonian ($|p_{+}\rangle,|p_{-}\rangle$)

[TABLE]

Here the diagonal terms are equal because $|p_{\pm}\rangle$ are two basis functions for a two-dimensional representation in contrast to $|p_{z}\rangle$ and $|p_{z}^{\prime}\rangle$ having identical, but different representations. Diagonalizing the total Hamiltonian $H^{\mathbf{\Gamma}v}$ and $H^{\mathbf{\Gamma}v}_{\text{strain}}$ we get the eigenenergies:

[TABLE]

showing that for biaxial strain $\epsilon_{xx}=\epsilon_{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 ($|p_{+}\rangle,|p_{-}\rangle$)

[TABLE]

Diagonalizing gives the energies at the $\mathbf{\Gamma}$ point (without SOC):

[TABLE]

We see that an in-plane electric field can split the $|p_{\pm}\rangle$ bands, whereas a perpendicular electric field shifts the two bands equally.

The Zeeman effect is described by the following Hamiltonian (in the basis states $|p_{+}\rangle\uparrow,|p_{-}\rangle\uparrow,|p_{+}\rangle\downarrow,|p_{-}\rangle\downarrow$):

[TABLE]

Due to $\sigma_{v}$ reflection symmetry in the $yz$ plane, the Landau effect vanishes between the states $|p_{+}\rangle,|p_{-}\rangle$.

DFT calculations reveal that the conduction band at the $\mathbf{\Gamma}$ point comprises two $s$ bands ($|s\rangle,|s^{\prime}\rangle$) described by the Hamiltonian

[TABLE]

Diagonalizing the Hamiltonian and expanding to lowest order in $k$ leads to the dispersion:

[TABLE]

where $i\in\{1,2\}$ denotes the two subbands. Since $N_{\pm}$ transform according to $\Gamma_{5}$ and $N_{z}$ according to $\Gamma_{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\rangle,|s^{\prime}\rangle$)

[TABLE]

Here we have included a second-order term proportional to $\epsilon_{+}\epsilon_{-}=(\epsilon_{xx}-\epsilon_{yy})^{2}+4\epsilon_{xy}^{2}$, since there is no first-order term and from DFT we see an effect of this within reasonable strain magnitudes. The lowest-order contributions to the energies at the $\mathbf{\Gamma}$ point are then given by:

[TABLE]

The effect of an electric field is described by the Hamiltonian

[TABLE]

Again, we diagonalize and expand to lowest order in the electric field, giving:

[TABLE]

Finally for a magnetic field we have both a Zeeman and a Landau contribution ($|s\uparrow\rangle,|s^{\prime}\uparrow\rangle,|s\downarrow\rangle,|s^{\prime}\downarrow\rangle$):

[TABLE]

[TABLE]

The symmetry analysis carried out here relies only on the symmetry of the crystal and the orbital characters of the relevant bands and therefore can be generalized to other materials sharing the same point group $D_{3h}$. In particular, we obtain using DFT that for hexagonal boron nitride (h-BN) the conduction band is a $s$ orbital and the valence band a $p_{z}$ orbital at the $\mathbf{K}$ point. At the $\mathbf{\Gamma}$ point, the lowest conduction bands are $p_{z}$ and $s$ states while the valence band is a pair of $p_{\pm}$ orbitals. Hence, the present Hamiltonian models for the valence bands at the $\mathbf{\Gamma}$ and $\mathbf{K}$ points can be directly applied to hexagonal boron nitride as well.

III 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 density functional theory (DFT) calculations using the projector augmented wave method as implemented in the GPAW packageMortensen et al. (2005); Enkovaara et al. (2010) with the Perdew-Burke-Ernzerhof approximation of the exchange-correlation energy Perdew et al. (1996). The selfconsistent calculation used a 16x16 Monkhorst-Pack $k$ point grid and a $1200\text{\,}\mathrm{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 Larsen et al. (2017), resulting in a lattice constant of $3.39\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. We use a vacuum layer of $20\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ 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 Ref. Olsen, 2016.

To extract the $k\cdot p$ parameters we perform bandstructure calculations on a dense $k$ point grid in the the vicinity of either $\mathbf{K}$ or $\mathbf{\Gamma}$. In Fig. 2, a contour plot of the valence and conduction bands at the $\mathbf{K}$ point is shown using both DFT and $k\cdot p$ results from Eq. (4b) with parameters fitted to the DFT energies. The model captures the trigonal warping of the bands described by the third-order terms in $\mathbf{q}$. We clearly see that the bands are symmetric under $q_{y}\rightarrow-q_{y}$ due to the combination of TR and the $\sigma_{v}$ reflection as explained in the previous section. The $k\cdot p$ model gives very good agreement with DFT energies up to a distance of $0.1\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$ from the $\mathbf{K}$ point, due to the inclusion of fourth-order terms in $\mathbf{q}$.

At the $\mathbf{\Gamma}$ point, $k\cdot p$ parameters are obtained by fitting Eq. (17) to the DFT result without SOC, and the SOC parameter is simply the energy splitting at the $\mathbf{\Gamma}$ point. The dispersions given by Eq. (20) are in good agreement with the DFT results within $0.1\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$ of $\mathbf{\Gamma}$ as seen in Fig. 3b. For the conduction band only the diagonal terms in Eq. (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 Fig. 3a.

III.1 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 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 $\mathbf{\Gamma}$ point shows high sensitivity to electric fields. At large electric fields the conduction band minimum changes from the $\mathbf{K}$ point to the $\mathbf{\Gamma}$ point at which point BAs becomes an indirect bandgap material. At electric fields larger than $0.75\text{\,}\mathrm{V}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$, the conduction band decreases below the valence band maximum at $\mathbf{K}$ and we get a transition to a metallic state. In comparison, an electric field induced bandgap of $0.25\text{\,}\mathrm{eV}$ has been demonstrated experimentallyZhang et al. (2009) in bilayer graphene subject to an electric displacement field of $0.3\text{\,}\mathrm{V}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$. We fit Eqs. (31a)-(31b) to the DFT results to obtain the relevant parameters.The second conduction band at $\mathbf{\Gamma}$ is not fitted since higher-lying bands are shifted down and mixed with this band.

III.2 Strain

To investigate the effect of strain in the linear approximation we modify the unit cell after the relaxation, $\mathbf{a}_{n}^{\prime}=\mathbf{a}_{n}+\bm{\epsilon}\mathbf{a}_{n}$ where $\mathbf{\epsilon}$ is the strain tensor, $\mathbf{a}_{n}$ are the equilibrium lattice vectors, and $\mathbf{a}_{n}^{\prime}$ 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 $\epsilon_{xx}=\epsilon_{yy}$, $\epsilon_{xx}=-\epsilon_{yy}$, and shear strain $\epsilon_{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 $\mathbf{k}$ space that reduces to $\mathbf{K}$ when strain goes to zero and we simply denote it $\mathbf{K}$. We limit our study to the cases where only one of the strain terms is nonzero. At the $\mathbf{K}$ point we find that the only linear strain terms that occur in the Hamiltonian, Eq. (7), are proportional to $\epsilon_{||}$. By fitting Eq. (8) to the DFT results we determine the parameters $D_{1}$ and $D_{3}$ and obtain excellent agreement within a range of $\epsilon_{||}=\pm 10\%$ as seen in Fig. 5a. In Fig. 5b we observe that the eigenvalues at the $\mathbf{K}$ point remain constant with both $\epsilon_{xx}-\epsilon_{yy}$ and $\epsilon_{xy}$, showing that the effects of strain is well described by the linear strain terms. For the valence band at the $\Gamma$ point we fit Eq. (22) to the DFT results. The resuls are plotted in Fig. 5, showing good agreement betwen 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 $\epsilon_{xx}-\epsilon_{yy}$ and shear strain $\epsilon_{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 $\epsilon_{||}$ as well, however, for $\epsilon_{xx}-\epsilon_{yy}$ and $\epsilon_{xy}$ there are no linear terms, but we find a significant parabolic dependence in Fig. 5b. Again, higher-order terms show little effect for large strains yet for most cases the second-order term gives a good approximation.

IV Conclusion

Compact and accurate $k\cdot 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 density functional theory 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\cdot p$ Hamlitonians and demonstrate excellent agreement with DFT results in the presence of strain and electric fields at two important points in the Brillouin zone, $\mathbf{K}$ and $\mathbf{\Gamma}$. DFT calculations reveal that the bandgap of 2D boron arsenide, located at the $\mathbf{K}$ point, becomes indirect at sufficiently large electric fields with the conduction band minimum located at the $\mathbf{\Gamma}$ 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\cdot 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.

Acknowledgments

MRB and MW gratefully acknowledge financial support from the Danish Council of Independent Research (Natural Sciences) grant no.: DFF-4181-00182. MRB is grateful for financial support from Otto Mønsted Foundation and Augustinus Foundation, as well as Beijing Institute of Nanoenergy and Nanosystems, Chinese Academy of sciences.