Effective Hamiltonians for phosphorene and 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

$$\begin{eqnarray}\begin{array}{rcl} {{k}_{+}} & = & {{k}_{x}}+i{{k}_{y}}, \\ {{k}_{-}} & = & {{k}_{x}}-i{{k}_{y}}, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

and similarly for other physical quantities. It can be seen that all the irreducible representations are either one-dimensional 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.

Table 1.  Group of wave vector D₃ at the K point (silicene). k_i are components of the wave vector, ${{\epsilon }_{ij}}$ of the strain tensor, B_i of an external magnetic field, E_i of an external electric field, and ${{{\bf J}}_{i}}$ of an angular momentum operator.

D3	E	2C3	$3C_{2}^{\prime }$	Odd	Even	Matrix
${{\Gamma }_{1}}$	1	1	1		$k_{x}^{2}+k_{y}^{2},k_{z}^{2},{{\epsilon }_{xx}}+{{\epsilon }_{yy}},{{\epsilon }_{zz}}$	${\bf 1}$
${{\Gamma }_{2}}$	1	1	−1	${{k}_{z}},{{B}_{z}}$	Ez	${{{\bf J}}_{z}}$
${{\Gamma }_{3}}$	2	−1	0	$({{k}_{-}},{{k}_{+}}),({{B}_{-}},-{{B}_{+}})$	$({{E}_{-}},{{E}_{+}})$	$({{{\bf J}}_{-}},-{{{\bf J}}_{+}})$

Table 2.  D₃ multiplication table of irreducible representations.

	${{\Gamma }_{1}}$	${{\Gamma }_{2}}$	${{\Gamma }_{3}}$
${{\Gamma }_{1}}$	${{\Gamma }_{1}}$	${{\Gamma }_{2}}$	${{\Gamma }_{3}}$
${{\Gamma }_{2}}$	${{\Gamma }_{2}}$	${{\Gamma }_{1}}$	${{\Gamma }_{3}}$
${{\Gamma }_{3}}$	${{\Gamma }_{3}}$	${{\Gamma }_{3}}$	${{\Gamma }_{1}}\oplus {{\Gamma }_{2}}\oplus {{\Gamma }_{3}}$

Table 3.  D₃ coupling constant table (${{\Gamma }_{3}}\;\otimes \;{{\Gamma }_{3}}={{\Gamma }_{1}}\oplus {{\Gamma }_{2}}\oplus {{\Gamma }_{3}}$). The notations $\Gamma _{3}^{1}$ and $\Gamma _{3}^{{{2}^{\prime }}}$ denote the first component of the first ${{\Gamma }_{3}}$ function and the second component of the second ${{\Gamma }_{3}}$ function, respectively, etc.

${{\Gamma }_{3}}\;\otimes \;{{\Gamma }_{3}}$
${{\Gamma }_{1}}$:	$\frac{1}{\sqrt{2}}\left( \Gamma _{3}^{1}\Gamma _{3}^{{{2}^{\prime }}}+\Gamma _{3}^{2}\Gamma _{3}^{{{1}^{\prime }}} \right)$
${{\Gamma }_{2}}$:	$\frac{i}{\sqrt{2}}\left( \Gamma _{3}^{1}\Gamma _{3}^{{{2}^{\prime }}}-\Gamma _{3}^{2}\Gamma _{3}^{{{1}^{\prime }}} \right)$
$\Gamma _{3}^{1}$:	$\Gamma _{3}^{2}\Gamma _{3}^{{{2}^{\prime }}}$
$\Gamma _{3}^{2}$:	$\Gamma _{3}^{1}\Gamma _{3}^{{{1}^{\prime }}}$

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 D₃ that

$$\begin{eqnarray}&&\mathop{\sum }\limits_{g}{{\chi }_{i}}\left( {{g}^{2}} \right)=n,\end{eqnarray} \tag{ 2 }$$

for all representations ${{\chi }_{i}},\;i=1,2,3$ where n = 6 is the order of the group D₃, g is a representation of a group element, and χ is the character of that group element. Further, since a spatial symmetry $\mathcal{R}$ exists such that $-{\bf k}=\mathcal{R}{\bf k}$ where ${\bf k}$ is the K point (and ${\bf k}\ne -{\bf k}$), all irreducible representations belong to the case a2 [12]. In this case, time-reversal symmetry additionally requires all Hamiltonian terms to satisfy

$$\begin{eqnarray}&&{{\mathcal{T}}^{-1}}\mathcal{H}({{\mathcal{R}}^{-1}}\mathcal{K})\mathcal{T}=\mathcal{H}{{(\zeta K)}^{*}},\end{eqnarray} \tag{ 3 }$$

where $\mathcal{T}$, $\mathcal{K}$, and $\mathcal{R}$ are the time-reversal operator, a ${{\Gamma }_{1}}$ invariant term, and a spatial symmetry operator establishing a link between basis functions at ${\bf k}$ and $-{\bf k}$, respectively. In the case of silicene, $\mathcal{T}={\bf 1}$ and $\mathcal{R}={{\mathcal{R}}_{x}}$ can be chosen where ${{\mathcal{R}}_{x}}$ denotes reflection in the x axis using a coordinate system where the x axis passes through 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 time-reversal 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,

$$\begin{eqnarray}&&\mathcal{H}={{\mathcal{H}}_{i}}+{{\mathcal{H}}_{e}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}_{e}}={{\mathcal{H}}^{\epsilon }}+{{\mathcal{H}}^{E}}+{{\mathcal{H}}^{B}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}_{i}} & = & {{a}_{1}}\left( {{k}_{y}}{{{\bf J}}_{x}}-{{k}_{x}}{{{\bf J}}_{y}} \right)+{{a}_{2}}\left( k_{x}^{2}+k_{y}^{2} \right) \\ {} & {} & +{{a}_{4}}{{k}_{y}}\left( 3k_{x}^{2}-k_{y}^{2} \right){{{\bf J}}_{z}}+{{a}_{5}}\left( k_{x}^{2}+k_{y}^{2} \right)\left( {{k}_{y}}{{{\bf J}}_{x}}-{{k}_{x}}{{{\bf J}}_{y}} \right)+\ldots , \\ \end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{\epsilon }} & = & {{e}_{1}}\left( {{\epsilon }_{xx}}+{{\epsilon }_{yy}} \right)+{{e}_{2}}{{\epsilon }_{zz}}+{{e}_{3}}\left( {{\epsilon }_{xx}}+{{\epsilon }_{yy}} \right)\left( {{k}_{y}}{{{\bf J}}_{x}}-{{k}_{x}}{{{\bf J}}_{y}} \right) \\ {} & {} & +{{e}_{4}}\left[ \left( {{\epsilon }_{xx}}-{{\epsilon }_{yy}} \right){{k}_{y}}+2{{\epsilon }_{xy}}{{k}_{x}} \right]{{{\bf J}}_{z}}+{{e}_{5}}{{\epsilon }_{zz}}\left( {{k}_{y}}{{{\bf J}}_{x}}-{{k}_{x}}{{{\bf J}}_{y}} \right)+\ldots , \\ \end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{E}} & = & {{c}_{1}}{{E}_{z}}{{{\bf J}}_{z}}+{{c}_{2}}\left( E_{x}^{2}+E_{y}^{2} \right)+{{c}_{3}}E_{z}^{2}+{{c}_{4}}\left( {{E}_{x}}{{k}_{y}}-{{E}_{y}}{{k}_{x}} \right) \\ {} & {} & +{{c}_{5}}\left[ \left( k_{x}^{2}-k_{y}^{2} \right){{E}_{x}}-2{{k}_{x}}{{k}_{y}}{{E}_{y}} \right]+{{c}_{6}}\left( E_{x}^{2}+E_{y}^{2} \right)\left( {{k}_{y}}{{{\bf J}}_{x}}-{{k}_{x}}{{{\bf J}}_{y}} \right) \\ {} & {} & +{{c}_{8}}\left[ \left( E_{x}^{2}-E_{y}^{2} \right){{k}_{y}}+2{{E}_{x}}{{E}_{y}}{{k}_{x}} \right]{{{\bf J}}_{z}}+\ldots , \\ \end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}^{B}}={{b}_{1}}\left( {{B}_{x}}{{{\bf J}}_{x}}+{{B}_{y}}{{{\bf J}}_{y}} \right)+{{b}_{2}}{{B}_{z}}{{{\bf J}}_{z}}+{{b}_{3}}\left( {{k}_{y}}{{B}_{x}}-{{k}_{x}}{{B}_{y}} \right)+\ldots ,\end{eqnarray} \tag{ 9 }$$

where the dots refer to higher-order terms. a_i, b_i, c_i and e_i are $k\cdot 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 ${{{\bf J}}_{i}}$ matrices represent the pseudospin degree of freedom and are the $(2\times 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₂ and a₃ terms provide quadratic in k contributions while the a₄ and a₅ are of cubic order but only the a₄ 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 ${{\mathcal{H}}_{i}}$ for graphene [8]:

$$\begin{eqnarray}&&{{\mathcal{H}}_{i}}={{a}_{61}}\left( {{k}_{y}}{{{\bf J}}_{x}}+{{k}_{x}}{{{\bf J}}_{y}} \right)+{{a}_{11}}\left( k_{x}^{2}+k_{y}^{2} \right)+{{a}_{62}}\left[ (k_{y}^{2}-k_{x}^{2}){{{\bf J}}_{x}}+2{{k}_{x}}{{k}_{y}}{{{\bf J}}_{y}} \right].\end{eqnarray} \tag{ 10 }$$

The sign difference in the a₆₁ term is due to a different choice of phase while the a₆₂ is the corresponding anisotropy term for graphene. It is seen here to be quadratic in k though it is cubic in k for silicene.

Zoom In Zoom Out Reset image size

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). Its factor group is symmorphic with ${{D}_{2h}}$ whose character table and associated functions are shown in table 4.

Zoom In Zoom Out Reset image size

Table 4.  Group of wave vector ${{D}_{2h}}$ at the Γ point (phosphorene).

${{D}_{2h}}$	E	${{C}_{2}}(z)$	${{C}_{2}}(y)$	${{C}_{2}}(x)$	i	${{\sigma }_{xy}}$	${{\sigma }_{xx}}$	${{\sigma }_{yz}}$	Odd	Even	Matrix
$\Gamma _{1}^{+}$	1	1	1	1	1	1	1	1		$k_{x}^{2},k_{y}^{2},k_{z}^{2},{{\epsilon }_{xx}},{{\epsilon }_{yy}},{{\epsilon }_{zz}}$	${\bf 1}$
$\Gamma _{2}^{+}$	1	−1	1	−1	1	−1	1	−1	By	${{k}_{x}}{{k}_{z}},{{\epsilon }_{xz}}$	${{{\bf J}}_{y}}$
$\Gamma _{3}^{+}$	1	1	−1	−1	1	1	−1	−1	Bz	${{k}_{x}}{{k}_{y}},{{\epsilon }_{xy}}$	${{{\bf J}}_{z}}$
$\Gamma _{4}^{+}$	1	−1	−1	1	1	−1	−1	1	Bx	${{k}_{y}}{{k}_{z}},{{\epsilon }_{yz}}$	${{{\bf J}}_{x}}$
$\Gamma _{1}^{-}$	1	1	1	1	−1	−1	−1	−1
$\Gamma _{2}^{-}$	1	−1	1	−1	−1	+1	−1	1	ky	Ey
$\Gamma _{3}^{-}$	1	1	−1	−1	−1	−1	1	1	kz	Ez
$\Gamma _{4}^{-}$	1	−1	−1	1	−1	1	1	−1	kx	Ex

The multiplication table of irreducible representations for ${{D}_{2h}}$ is given in table 5. Note that the parity of the product representation follows the simple rule: $\Gamma _{1}^{+}\;\otimes \;\Gamma _{2}^{+}=\Gamma _{2}^{+}$, $\Gamma _{1}^{+}\;\otimes \;\Gamma _{2}^{-}=\Gamma _{2}^{-}$, $\Gamma _{1}^{-}\;\otimes \;\Gamma _{2}^{+}=\Gamma _{2}^{-}$, $\Gamma _{1}^{-}\;\otimes \;\Gamma _{2}^{-}=\Gamma _{2}^{+}$ etc. Since the irreducible representations are all one-dimensional, the coupling constants of product representations are trivial. Further, since ${\bf k}=-{\bf 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].

Table 5.  ${{D}_{2h}}$ multiplication table of irreducible representations.

	${{\Gamma }_{1}}$	${{\Gamma }_{2}}$	${{\Gamma }_{3}}$	${{\Gamma }_{4}}$
${{\Gamma }_{1}}$	${{\Gamma }_{1}}$	${{\Gamma }_{2}}$	${{\Gamma }_{3}}$	${{\Gamma }_{4}}$
${{\Gamma }_{2}}$	${{\Gamma }_{2}}$	${{\Gamma }_{1}}$	${{\Gamma }_{4}}$	${{\Gamma }_{3}}$
${{\Gamma }_{3}}$	${{\Gamma }_{3}}$	${{\Gamma }_{4}}$	${{\Gamma }_{1}}$	${{\Gamma }_{2}}$
${{\Gamma }_{4}}$	${{\Gamma }_{4}}$	${{\Gamma }_{3}}$	${{\Gamma }_{2}}$	${{\Gamma }_{1}}$

We are now in a position to write down the most general Hamiltonian for phospherene allowed by symmetry. The result is

$$\begin{eqnarray}&&\mathcal{H}={{\mathcal{H}}_{i}}+{{\mathcal{H}}_{e}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}_{e}}={{\mathcal{H}}^{\epsilon }}+{{\mathcal{H}}^{E}}+{{\mathcal{H}}^{B}}+{{\mathcal{H}}^{mix}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}_{i}}={{a}_{1}}k_{x}^{2}+{{a}_{2}}k_{y}^{2}+\mathop{\sum }\limits_{i\leqslant j}{{a}_{ij}}k_{i}^{2}k_{j}^{2}+\ldots ,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{\epsilon }} & = & \left( {{e}_{1}}+{{e}_{2}}k_{x}^{2}+{{e}_{3}}k_{y}^{2} \right){{\epsilon }_{xx}}+\left( {{e}_{5}}+{{e}_{6}}k_{x}^{2}+{{e}_{7}}k_{y}^{2} \right){{\epsilon }_{yy}} \\ {} & {} & +\left( {{e}_{9}}+{{e}_{10}}k_{x}^{2}+{{e}_{11}}k_{y}^{2} \right){{\epsilon }_{zz}}+\ldots , \\ \end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{E}} & = & \left( {{c}_{4}}+{{c}_{5}}k_{x}^{2}+{{c}_{6}}k_{y}^{2} \right)E_{x}^{2}+\left( {{c}_{8}}+{{c}_{9}}k_{x}^{2}+{{c}_{10}}k_{y}^{2} \right)E_{y}^{2} \\ {} & {} & +\left( {{c}_{12}}+{{c}_{13}}k_{x}^{2}+{{c}_{14}}k_{y}^{2} \right)E_{z}^{2}+\ldots , \\ \end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{B}} & = & {{b}_{4}}{{k}_{x}}{{B}_{x}}+{{b}_{5}}{{k}_{y}}{{B}_{y}} \\ {} & {} & +\left( {{b}_{7}}+{{b}_{8}}k_{x}^{2}+{{b}_{9}}k_{y}^{2} \right)B_{x}^{2}+\left( {{b}_{11}}+{{b}_{12}}k_{x}^{2}+{{b}_{13}}k_{y}^{2} \right)B_{y}^{2} \\ {} & {} & +\left( {{b}_{15}}+{{b}_{16}}k_{x}^{2}+{{b}_{17}}k_{y}^{2} \right)B_{z}^{2}+\ldots , \\ \end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}^{mix}}={{m}_{2}}{{B}_{y}}{{E}_{z}}{{k}_{x}}+{{m}_{3}}{{B}_{z}}{{E}_{x}}{{k}_{y}}+\ldots \end{eqnarray} \tag{ 17 }$$

Our intrinsic Hamiltonian agrees with the result of Li and Appelbaum [11].

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].

3.1.1. Intrinsic

For an arbitrary direction and keeping terms up to cubic in k, the band structure is given by diagonalizing the H_i Hamiltonian. Writing ${{k}^{2}}=k_{x}^{2}+k_{y}^{2}$, one gets

$$\begin{eqnarray}&&E={{a}_{2}}{{k}^{2}}\pm \sqrt{{{\left( {{a}_{1}}+{{a}_{5}}{{k}^{2}} \right)}^{2}}{{k}^{2}}+a_{4}^{2}k_{y}^{2}{{\left( 3k_{x}^{2}-k_{y}^{2} \right)}^{2}}}.\end{eqnarray} \tag{ 18 }$$

Thus, our equation above shows that the intrinsic band dispersion can have quadratic, cubic, ... terms in the wave vector.

The a₁ coefficient is related to the Fermi velocity, a₂ is related to an effective mass, a₄ gives rise to anisotropy, and a₅ leads to cubic terms. Along the $K\to \Gamma$ direction (k_y = 0), the band structure simplifies to

$$\begin{eqnarray}&&E(k)=\pm {{a}_{1}}k+{{a}_{2}}{{k}^{2}}\pm {{a}_{5}}{{k}^{3}}.\end{eqnarray} \tag{ 19 }$$

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₄ is absent along this direction.

Zoom In Zoom Out Reset image size

We have performed a cubic fit of the bands obtained using VASP up to about half-way from K. Rewriting equation (19) as

$$\begin{eqnarray}&&E(k)=\pm \hbar {{v}_{F}}k+\frac{{{\hbar }^{2}}{{k}^{2}}}{2m*}\pm \gamma {{k}^{3}},\end{eqnarray} \tag{ 20 }$$

and performing a fit to

$$\begin{eqnarray}&&y(n)=C+(B1)n+(B2){{n}^{2}}+(B3){{n}^{3}},\end{eqnarray} \tag{ 21 }$$

where n is chosen to be 100 points between K and G, we obtain (with $\Delta k=(2\pi )/(100\sqrt{3}{{a}_{0}})$)

$$\begin{eqnarray}&&{{v}_{F}}=\frac{B1}{\hbar \Delta k}=4.188\times {{10}^{6}}{{a}_{0}}B1\;{\rm m}\;{{{\rm s}}^{-1}},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&m*=\frac{{{\hbar }^{2}}\Delta {{k}^{2}}}{2B2}=\frac{2.29\times {{10}^{-4}}}{B2}{{m}_{0}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\gamma =\frac{B3}{\Delta {{k}^{3}}}={{10}^{6}}\;{{\left( \frac{\sqrt{3}{{a}_{0}}}{2\pi } \right)}^{3}}B3=1.186\times {{10}^{6}}B3\;{\rm eV}{{\overset{\circ}{A} }^{-3}},\end{eqnarray} \tag{ 24 }$$

and where we have used ${{a}_{0}}=3.84\;$Å.

The results of the curve fitting are given in table 6.

Table 6.  Band fitting for silicene along K to Γ.

	B1	vF (${{10}^{5}}\;{\rm m}\;{{{\rm s}}^{-1}}$)	B2	${{m}^{*}}/{{m}_{0}}$	B3	$\gamma ({\rm eV}{{\overset{\circ}{A} }^{-3}})$
VB	−0.040	6.43	$-3.461\times {{10}^{-4}}$	−0.66	$5.608\times {{10}^{-6}}$	6.65
CB	0.035	5.63	$6.315\times {{10}^{-4}}$	0.36	$-7.173\times {{10}^{-6}}$	−8.51

Figure 3 (right) shows the dispersion along $K\to 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.

3.1.2. Strain

In the presence of strain, the Hamiltonian at the K point is given by

$$\begin{eqnarray}&&{{H}_{K}}={{e}_{1}}\left( {{\epsilon }_{xx}}+{{\epsilon }_{yy}} \right)+{{e}_{2}}{{\epsilon }_{zz}}\equiv \Delta (\epsilon ).\end{eqnarray} \tag{ 25 }$$

First, we note that applying strain perpendicular to the plane (${{\epsilon }_{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., ${{\epsilon }_{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₁; 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–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.

For finite k and an arbitrary strain, one can also diagonalize the Hamiltonian to give

$$\begin{eqnarray}&&E(k,\epsilon )=\Delta (\epsilon )\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray*}&&\pm {{\left\{ e_{4}^{2}{{\left[ \left( {{\epsilon }_{xx}}-{{\epsilon }_{yy}} \right){{k}_{y}}+2{{\epsilon }_{xy}}{{k}_{x}} \right]}^{2}}+{{\left[ {{e}_{3}}\left( {{\epsilon }_{xx}}+{{\epsilon }_{yy}} \right)+{{e}_{5}}{{\epsilon }_{zz}} \right]}^{2}}\left( k_{x}^{2}-k_{y}^{2} \right) \right\}}^{1/2}}.\end{eqnarray*}$$

3.1.3. Electric field

In the presence of an electric field, the Hamiltonian at the K point is given by

$$\begin{eqnarray}&&{{H}_{K}}={{c}_{1}}{{E}_{z}}{{{\bf J}}_{z}}+{{c}_{2}}\left( E_{x}^{2}+E_{y}^{2} \right)+{{c}_{3}}E_{z}^{2}.\end{eqnarray} \tag{ 27 }$$

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}_{1}}\ne 0$) 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}_{1}}=0.037\;$ eÅ. No existing calculations allow us to determine c₂ and c₃.

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 ${{E}_{y}}={{E}_{z}}=1$ (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.

Zoom In Zoom Out Reset image size

3.1.4. 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

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{B}} & = & {{a}_{1}}\left( {{k}_{y}}{{{\bf J}}_{x}}-{{k}_{x}}{{{\bf J}}_{y}} \right)+{{b}_{1}}\left( {{k}_{y}}{{B}_{x}}-{{k}_{x}}{{B}_{y}} \right)+{{b}_{2}}{{B}_{z}}{{{\bf J}}_{z}} \\ {} & {} & +{{b}_{3}}\left( {{B}_{x}}{{{\bf J}}_{x}}+{{B}_{y}}{{{\bf J}}_{y}} \right). \\ \end{array}\end{eqnarray} \tag{ 28 }$$

Thus, at the K point, we have

$$\begin{eqnarray}&&{{H}_{K}}={{b}_{2}}{{B}_{z}}{{{\bf J}}_{z}}+{{b}_{3}}\left( {{B}_{x}}{{{\bf J}}_{x}}+{{B}_{y}}{{{\bf J}}_{y}} \right).\end{eqnarray} \tag{ 29 }$$

Similarly to the electric field case, the magnetic-field problem is separable and only the soluble part is plotted in figure 4.

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).

Zoom In Zoom Out Reset image size

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].

3.2.1. 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 $\sim 20$%. The leading terms in the strain Hamiltonian are

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}^{\epsilon }} & = & \left( {{e}_{1}}+{{e}_{2}}k_{x}^{2}+{{e}_{3}}k_{y}^{2} \right){{\epsilon }_{xx}}+\left( {{e}_{5}}+{{e}_{6}}k_{x}^{2}+{{e}_{7}}k_{y}^{2} \right){{\epsilon }_{yy}} \\ {} & {} & +\left( {{e}_{9}}+{{e}_{10}}k_{x}^{2}+{{e}_{11}}k_{y}^{2} \right){{\epsilon }_{zz}}+{{e}_{13}}{{k}_{x}}{{k}_{y}}{{\epsilon }_{xy}}. \\ \end{array}\end{eqnarray} \tag{ 30 }$$

At the Γ point, the energies are

$$\begin{eqnarray}&&{{E}_{i}}(\epsilon )={{E}_{i}}+{{e}_{1}}{{\epsilon }_{xx}}+{{e}_{5}}{{\epsilon }_{yy}}+{{e}_{9}}{{\epsilon }_{zz}},\end{eqnarray} \tag{ 31 }$$

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 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

$$\begin{eqnarray*}&&\mathcal{H}={{E}_{i}}(\epsilon )+({{a}_{1}}+{{e}_{2}}{{\epsilon }_{xx}}+{{e}_{6}}{{\epsilon }_{yy}}+{{e}_{10}}{{\epsilon }_{zz}})k_{x}^{2}\end{eqnarray*}$$

$$\begin{eqnarray}&&({{a}_{2}}+{{e}_{3}}{{\epsilon }_{xx}}+{{e}_{7}}{{\epsilon }_{yy}}+{{e}_{11}}{{\epsilon }_{zz}})k_{y}^{2}.\end{eqnarray} \tag{ 32 }$$

3.2.2. 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

$$\begin{eqnarray}&&{{E}_{\Gamma }}={{E}_{i}}+{{c}_{12}}E_{z}^{2}.\end{eqnarray} \tag{ 33 }$$

DFT results just published appear to confirm our prediction [28].

3.2.3. 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 ${\bf B}=(0,0,{{B}_{z}})$ with a vector potential ${\bf A}=(0,{{B}_{z}}x,0)$, the Hamiltonian is given by [7]

$$\begin{eqnarray}&&H=\frac{p_{x}^{2}}{2m_{x}^{*}}+\frac{1}{2m_{y}^{*}}{{\left( {{p}_{y}}+e{{B}_{z}}x \right)}^{2}}.\end{eqnarray} \tag{ 34 }$$

The eigenfunctions can be separated,

$$\begin{eqnarray}&&f({\bf r})=\frac{1}{\sqrt{{{L}_{y}}}}{{e}^{i{{k}_{y}}y}}\phi (x),\end{eqnarray} \tag{ 35 }$$

and the eigenvalues are equidistant Landau levels

$$\begin{eqnarray}&&{{E}_{n}}=\left( n+\frac{1}{2} \right)\hbar {{\omega }_{c}},\end{eqnarray} \tag{ 36 }$$

where ${{\omega }_{c}}=e{{B}_{z}}/\sqrt{m_{x}^{*}m_{y}^{*}}$. Thus, the prediction is that, to the lowest order, a perpendicular magnetic field will change the energy levels linearly.