Magnetoelectronic properties of normal and skewed phosphorene nanoribbons

TL;DR

This paper investigates the electronic properties of phosphorene nanoribbons under magnetic and electric fields, providing analytical models and numerical validation for their energy spectra, edge states, and Landau level formation.

Contribution

It introduces a new continuous model for analyzing phosphorene nanoribbons and derives analytical expressions for Landau levels and flat-band segments, enhancing understanding of their magnetoelectronic behavior.

Findings

Analytical expressions for Landau level formation are derived.

Wavefunctions for flatbands are obtained and validated.

Edge states shift with magnetic field as analytically described.

Abstract

The energy spectrum and eigenstates of single-layer black phosphorous nanoribbons in the presence of perpendicular magnetic field and in-plane transverse electric field are investigated by means of a tight-binding method and the effect of different types of edges is analytically examined. A description based on a new continuous model is proposed by expansion of the tight-binding model in the long-wavelength approximation. The wavefunctions corresponding to the flatband part of the spectrum are obtained analytically and are shown to approach agree well with the numerical results from the tight-binding method. Analytical expressions for the critical magnetic field at which Landau levels are formed and the ranges of wavenumbers in the dispersionless flat-band segments in the energy spectra are derived. We examine the evolution of the Landau levels when an in-plane lateral electric field is applied and determine analytically how the edge states shift with magnetic field.

Figures

Tables

Table 1. Table 1: The effective masses along the main axis for the electron and the hole, where m 0 subscript 𝑚 0 m_{0} is the free electron mass.

	$x$-ZZ	$y$-AC
$m_e^{\ast }$ (the electron)	$1.1547m_0$	$0.1951m_0$
	$0.848m_0^{♣}$	$0.167m_0^{♣}$
$m_h^{\ast }$ (the hole)	$3.2279m_0$	$0.1651m_0$
	$1.142m_0^{♣}$	$0.184m_0^{♣}$
♣results from Ref. [zhou15, ]

Equations

$$H=\sum _{i\ne j}{t}_{ij}{e}^{i{\phi }_{ij}}{c}_i^{†}{c}_j,$$

$$H_{k}^{4\times 4}=\left(\begin{array}[]{cccc}U_{A}&H_{AB}&H_{AD}&H_{AC}\\ H^{*}_{AB}&U_{B}&H_{BD}&H_{BC}\\ H^{*}_{AD}&H^{*}_{BD}&U_{D}&H_{DC}\\ H^{*}_{AC}&H_{BC}^{*}&H_{DC}^{*}&U_{C}\\ \end{array}\right),$$

$$U=U_{A-D}=2t_3\cos (a_x{k}_x)+2t_7\cos (a_y{k}_y)$$

$$+4t_{10}\cos (a_x{k}_x)\cos (a_y{k}_y).$$

$$H_{AB}$$

$$+2t_9\cos (a_x{k}_x)e^{-i{k}_y(a_y-h_{2y})},$$

$$H_{AD}$$

$$H_{AC}$$

$$+2\left(t_1\cos \left(\frac{a_x{k}_x}{2}\right)+t_8\cos \left(\frac{3a_x{k}_x}{2}\right)\right)e^{-i{k}_y{h}_{1y}}.$$

$$H_{\mathbf{k}}^{2\times 2}=\left[\begin{array}[]{cc}U+H_{AD}&H_{AB}+H_{AC}\\ H^{*}_{AB}+H^{*}_{AC}&U+H_{AD}\end{array}\right],$$

$$\psi=\left[\begin{array}[]{cc}(\phi_{A}+\phi_{D})/2\\ (\phi_{B}+\phi_{C})/2\end{array}\right]=\left[\begin{array}[]{cc}\phi_{1}\\ \phi_{2}\end{array}\right].$$

$$E=U+H_{AD}\pm |H_{AB}+H_{AC}|,$$

$$H_{{\mathbf{k}}\rightarrow\Gamma}^{2\times 2}=\left[\begin{array}[]{cc}f&g\\ g^{\dagger}&f\end{array}\right],$$

$$E_0$$

$$E_g$$

$${\chi }_x$$

$${\chi }_y$$

$${\alpha }_x$$

$${\alpha }_y$$

$$-\left(t_6/2+t_9\right){\left(a_y-h_{2y}\right)}^2=2.026\mathrm{eV}{\mathring{\mathrm{A}}}^2,$$

$$\gamma$$

$$+(t_6+2t_9)(a_y-h_{2y})=-5.952\mathrm{eV}\mathring{\mathrm{A}},$$

$$E$$

$$U=\frac{1}{\sqrt{2}}\left[\begin{array}[]{cc}1&1\\ 1&-1\end{array}\right].$$

$$H_{{\mathbf{k}}\rightarrow\Gamma}^{U}=\left[\begin{array}[]{cc}E_{c}&i\gamma k_{y}\\ -i\gamma k_{y}&E_{v}\end{array}\right],$$

$$E_c$$

$$E_v$$

$$E_c^{eff}$$

$$E_v^{eff}$$

$$m_{e(h),x}^{\ast }=m_{0e(0h),x}^{\ast }=\frac{{\hslash }^2}{2({\alpha }_x\pm {\chi }_x)},$$

$$m_{e(h),y}^{\ast }=\frac{{\hslash }^2}{2\left({\alpha }_y\pm {\chi }_y+\frac{{\gamma }^2}{E_g}\right)}.$$

$$\mathbf{d}=c_x{a}_x{\mathbf{e}}_x+c_y{a}_y{\mathbf{e}}_y.$$

$$\frac{1}{m_{e(h),x^{\prime }}^{\ast }}=\frac{1}{m_{e(h),x}^{\ast }}{\cos }^2\phi +\frac{1}{m_{e(h),y}^{\ast }}{\sin }^2\phi ,$$

$$\frac{1}{m_{e(h),y^{\prime }}^{\ast }}=\frac{1}{m_{e(h),x}^{\ast }}{\sin }^2\phi +\frac{1}{m_{e(h),y}^{\ast }}{\cos }^2\phi .$$

$$k_x$$

$$k_y$$

$$\frac{{\hslash }^2}{2m_{e(h),y^{\prime }}^{\ast }}{\left[-i\frac{\partial }{\partial \tilde{y}}-\tilde{y}\left(\frac{m_{e(h),y^{\prime }}^{\ast }}{m_{e(h),x^{\prime }}^{\ast }}-1\right)\frac{\tan 2\phi }{2l_B^2}\right]}^2$$

$$+\frac{1}{2}{m}_{e(h),y^{\prime }}^{\ast }{\omega }_{e(h)}^2{\tilde{y}}^2=\Delta E_{c(v)},$$

$${\varphi }_{c(v)}(\tilde{y})=\exp \left(-\frac{m_{e(h),y^{\prime }}^{\ast }{\omega }_{e(h)}}{\hslash }\frac{{\tilde{y}}^2}{2}\right)$$

Full text

I Introduction

Phosphorene, a single layer of black phosphorus, is a recently isolated material,liu14 ; lu14 providing some extraordinary advantages over other two-dimensional (2D) materials. Most importantly, it has a large band gap of about $1.84$ eV,rudenko15 thus circumventing the main drawback of graphene with its zero band gap. Additionally, because of the puckered structure liu14 ; fei14 ; aierken15 it has highly anisotropic electronic and thermal properties, allowing for the emergence of new functionalities. On top of this, higher carrier mobilities can be achieved as compared to those 2D materials consisting of transition metal dichalcogenides.qiao14

However, phosphorene also has two major downsides: i) the quality of the samples degrades with time when exposed to air,favron14 and ii) unlike few-layer phosphorene the band gap is not electrically tunable. The first issue can be evaded successfully by for instance encapsulating phosphorene with hexagonal boron nitride, which can has the added benefit of increasing the crystal quality and carrier mobility, as in the case of graphene.cao15 ; li15 ; chen15 Besides strain,rodin14 ; cakir14 fashioning phosphorene into nanoribbons can modify the band gap through the quantum confinement effect. Indeed, several theoretical papers have studied various types of phosphorene nanoribbons, and all of them show potentials for certain applications.ezawa14 ; sisakht15 ; grujic16 Note that due to the anisotropy of phosphorene there are two distinct types of zigzag (ZZ) and armchair (AC) edges, normal and skewed, with skewed edges intersecting the ridges of phosphorus atoms at a sharp angle.grujic16 Intriguingly, normal zigzag and skewed armchair (sAC) nanoribbons are metallic, while skewed zigzag (sZZ) and normal armchair nanoribbons are insulating, and all of them have an electrically tunable band structure, featuring metal-insulator transitions for a sufficiently strong electric field.

While it was predicted that ”bulk phosphorene” displays linearly dispersing Landau levels (LLs) for low-energy quasiparticles in low magnetic fields,pereira15 ; ostahie15 ; zhou15 which was confirmed by recent experiments,li15 ; chen15 the results concerning the magnetic response of phosphorene nanoribbons are relatively scarce. In Ref.  [zhou15, ] the formation of LLs in normal nanoribbons was reported, while in Ref.  [ostahie15, ] the impact of magnetic field on the quasi-flat bands (QFBs) of zigzag nanoribbons were investigated. The purpose of this paper is to examine the influence of magnetic field on phosphorene nanoribbons in more details. In particular, we will examine the band structure of phosphorene nanoribbons with various edge types (including skewed). We will discuss the impact of dispersion anisotropy on LL formation from both the qualitative and quantitative point of view. We will also briefly consider the effect of in-plane electric field on the LLs.

II Theoretical model

We model the phosphorene nanoribbons using the tight-binding model which is defined by $10$ hopping parameters rudenko15 which are indicated in Fig. 1. Within the tight-binding approximation the Hamiltonian reads

[TABLE]

where the summation runs over all the lattice sites of phosphorene, $t_{ij}$ are the hopping parameters, $\varphi_{ij}=\frac{e}{\hbar}\int_{\mathbf{r_{j}}}^{\mathbf{r_{i}}}\mathbf{A}\cdot d\mathbf{l}$ denotes the Peierls phase picked up while hopping in the presence of the magnetic field, and $c_{j}^{\dagger}(c_{i})$ is the creation (annihilation) operator of an electron on the site j(i).

The unit cell of the single phosphorene sheet, framed by the solid square line in Fig. 1, contains four atoms labeled by A (blue), B (green), C (yellow) and D (black). Using the tight binding model, the four-band Hamiltonian ezawa14 ; pereira15 for single-layer phosphorene sheet is given by:

[TABLE]

with eigenvectors represented by spinor $[\phi_{A}\phi_{B}\phi_{D}\phi_{C}]^{T}$, and diagonal terms are given by

[TABLE]

The interaction terms between sublattice sites are

[TABLE]

Here ${h}_{1y}=d_{1}\cos(\alpha_{hex}/2)$, ${{h}_{2y}}=d_{2}\cos\beta$ (see Fig. 1). The parameter values are taken from Ref.  [rudenko15, ]. The equality of certain terms of the Hamiltonian (2) indicates the ”equivalence” of certain atomic sites. This is actually due to the $D_{2h}$ point group invariance ezawa14 ; voon15 , which is preserved even when a perpendicular magnetic field is applied. Atoms in the upper and lower lattice, namely $A\equiv D$ and $B\equiv C$, are ”indistinguishable” and the unit cell is reduced to a single dimmer, framed by the dashed square in Fig. 1. Using this symmetry argument, the simplified two-band Hamiltonian reads

[TABLE]

which acts upon the spinors

[TABLE]

The eigenvalue problem for the Hamiltonian (7) can be solved analytically

[TABLE]

where the upper (lower) sign is for the conduction (valence) band.

Recent theoretical analysis ezawa14 ; pereira15 based on the tight binding model with five hopping parameters shows that the Maclaurin series of analytical functions in the vicinity of the $\Gamma$-point up to the quadratic term give satisfactory accuracy for the phosphorene band structure. Using a similar argument we derive the effective two-band Hamiltonian

[TABLE]

where $f=E_{0}+\chi_{x}k_{x}^{2}+\chi_{y}k_{y}^{2}$, and $g=E_{g}/2-i\gamma k_{y}+\alpha_{x}k_{x}^{2}+\alpha_{y}k_{y}^{2}$, with

[TABLE]

and obtain the simplified dispersion for electrons and holes

[TABLE]

From Eq. (18) it is obvious that, phosphorene has a non-elliptical dispersion due to the parameter $\gamma$.

In order to transform the low-energy Hamiltonian given by Eq. (10) into a more convenient form we perform the unitary transformation given by the Hadamard matrix

[TABLE]

The resulting Hamiltonian has the same form as the one obtained in Ref.  [rodin14, ] by fitting to results from ab-initio calculations

[TABLE]

where diagonal terms

[TABLE]

are related to the conduction and valence band dispersions, respectively, while the off-diagonal term introduces interband coupling. Here, $E_{c(v)0}$ are the conduction (valence) band edges in monolayer phophorene, $m^{*}_{0e,x(y)}$ and $m^{*}_{0h,x(y)}$ are the electron and hole effective masses when the coupling is neglected. Furthermore, spinors given by Eq. (8) are transformed to $\psi^{U}=\left[\phi_{c}\phi_{v}\right]^{T}$, where $\phi_{c(v)}=(\phi_{1}\pm\phi_{2})/\sqrt{2}$.

In the long-wavelength limit $E_{c}-E_{v}>E_{g}\gg|\gamma k_{y}|$. Thus, to lowest order the Eq. (18) is expanded around the $\Gamma$-point and the energy dispersions of conduction and valence bands are approximated by $E^{eff}_{c}=E_{c}+\gamma^{2}k_{y}^{2}/E_{g}$ and $E^{eff}_{v}=E_{v}-\gamma^{2}k_{y}^{2}/E_{g}$, respectivelyzhou15 . Therefore, the Hamiltonian given by Eq. (20) is reduced to the diagonal form $H_{k}^{eff,U}\rightarrow{\rm diag}(E^{eff}_{c},E^{eff}_{v})$ where

[TABLE]

The effective masses along the $\Gamma-X$ direction remain unchanged,

[TABLE]

while the effective masses along the $\Gamma-Y$ direction are modified by the term that perturbatively takes into account the interband coupling,

[TABLE]

The upper ”+” (lower ”$-$”) sign is for the electron (hole). The effective masses along the main axes (see Fig. 1) are given in Tab. I. Note that there is a difference between our results and the one obtained from the previously derived continuum model zhou15 based on 5 hopping parametersrudenko14 .

Our goal is to analyze ribbons with arbitrary edges, where the translation vector of a unit cell is

[TABLE]

Here, $c_{x}$ and $c_{y}$ are mutually prime integers, ${\mathbf{e}}_{x(y)}$ are the unit vectors, while $d=|{\mathbf{d}}|$ is the length of the unit cell. Since ribbon edges are not flat, we define the effective ribbon width as the total square of hexagonal plaquettes inside the unit cell divided by the unit cell width $d$.

It is straightforward to show that the effective masses of electrons and holes in the vicinity of the $\Gamma$-point along the $x^{\prime}$ and $y^{\prime}$ axes, which are directed along and perpendicular to the ribbon edge, respectively, can be estimated as:

[TABLE]

Here $\varphi=\arctan(c_{y}a_{y}/c_{x}a_{x})$ is the angle of rotation with respect to the $(x,y)$ coordinate system shown in Fig. 1. In order to write the Hamiltonian in the rotated frame we substitute

[TABLE]

in Eqs. (23) and (24). For infinitely long nanoribbons with arbitrary edges the Hamiltonian should not depend on $x^{\prime}$, therefore $\phi_{c(v)}(x^{\prime},y^{\prime})=\phi_{c(v)}(y^{\prime})\cdot e^{ik_{x^{\prime}}x^{\prime}}$. For convenience we set $y^{\prime}=0$ in the middle of the ribbon.

We include the magnetic field perpendicular to the structure (${\mathbf{B}}=B{\mathbf{e}}_{z}$). One may show that in the continuum model $k_{x^{\prime}}\rightarrow k_{x^{\prime}}-y^{\prime}/l_{B}^{2}$, where we adopt the Landau gauge ${\mathbf{A}}=(-By^{\prime},0,0)$, with $l_{B}=\sqrt{\hbar/eB}$ being the magnetic field length and operators $k_{x^{\prime}}=-i\partial/\partial x^{\prime}$ and $k_{y^{\prime}}=-i\partial/\partial y^{\prime}$. Using Eqs. (28,29) and after some elaborate algebra we obtain a set of decoupled differential equations for the conduction and valence band,

[TABLE]

where the subscripts $c$ and $e$ are used for the conduction band, and $v$ and $h$ denote the valence band. Here, $\tilde{y}=y^{\prime}-k^{\prime}_{x}l^{2}_{B}$ is the shifted $y^{\prime}$ coordinate, $\omega_{e(h)}=eB/\sqrt{m^{*}_{e(h),x}m^{*}_{e(h),y}}$ are the cyclotron frequencies for the electron (hole) in bulk phosphorene, and $\Delta E_{c(v)}=\pm(E_{c(v)}-E_{c(v)0})$ are eigenenergies measured from the bottom (top) of the conduction (valence) band edges. One may notice that the above equation does not include the term related to the band offset at the edges of the ribbon. This is because we want to analyze the conditions under which LLs exist in a ribbon. Namely, we insist that, either the ribbon is sufficiently wide, or the magnetic field is sufficiently strong, so that particle localization is governed by the effective parabolic potential $V_{e(h)}(k_{x^{\prime}},y^{\prime},B)=\frac{1}{2}m_{e(h),y^{\prime}}\omega_{e(h)}^{2}\left(y^{\prime}-l^{2}_{B}k_{x}^{\prime}\right)^{2}$ given in Eq. (30), rather than by the ribbon edges.

We seek the eigenvector components in the form

[TABLE]

which leads to the differential equation

[TABLE]

where $\xi=\xi_{e(h)}=\sqrt{m^{*}_{e(h),y^{\prime}}\cdot\omega_{e(h)}/\hbar}\cdot\tilde{y}$ is the dimensionless coordinate and $\varepsilon_{e(h)}=2\Delta E_{c(v)}/\hbar\omega_{e(h)}$ is the dimensionless energy. Finally, to get eigenvalues denoted by the integer quantum number we impose the condition $\varepsilon_{c(v)}-1=2n_{e(h)}$. Thus, the solutions of Eq. (32) are Hermite polynomials

[TABLE]

where the principle quantum number $n\triangleq n_{e(h)}=\{0,1,2,...\}$ is the LL number, the normalization constant is $C_{n}=\sqrt[4]{m^{*}_{e(h),y^{\prime}}\cdot\omega_{e(h)}}/\sqrt{{n!}2^{n}\pi^{1/2}\hbar^{1/2}}$, and the eigenvalues of the LLs follow the quantization of a one-dimensional quantum harmonic oscillator (QHO)

[TABLE]

It is obvious that the separation between adjacent energy levels ($|\Delta n_{e(h)}|=1$) is approximately $\Delta E_{c(v)}=\hbar\omega_{e(h)}$. This energy difference is the same for ribbons independent of edges and is equal to the one found in single-layer phosphorous sheetzhou15 ; pereira15 .

It should be pointed out that the exact treatment of LLs in bulk phosphorene shows that the off-diagonal terms, that accounts for the conduction and valence band coupling, are similar in form as the Rasshba (Dresselhaus) spin-orbit interaction in conventional semiconductorswinkler03 . However, these terms do not contribute significantly to the spectra of the lowest energy states and the spatial density distribution corresponding to the first few LLs are found to have elliptical shapezhou15 . In a quasi-classical picture we might infer that, when the magnetic field is turned on perpendicular to the bulk phosphorene sheet, electrons and holes undergo elliptical cyclotron orbits.

Let us now consider the impact of a perpendicular magnetic field on the formation of LLs. When the magnetic field is turned on, the particle is essentially confined in the transversal direction by a restricted effective parabolic potential $V_{e(h)}(k_{x^{\prime}},y^{\prime},B)$. One may notice that on increase (decrease) of $k_{x^{\prime}}$ shifts the effective parabolic potential towards the upper (lower) ribbon edge, for both the electron and the hole. Therefore, states occupying positive and negative momenta reside on opposite sides of the ribbon. Also, the magnetic field does not split oppositely charged particles in the transversal direction, and the electron and the hole will have similar localization in space even when the magnetic field is turned on. Furthermore, the effective potential is proportional to the transverse effective mass $m^{*}_{e(h),y^{\prime}}$. Thus, the confinement will be stronger in ribbons with higher effective mass in the transversal direction.

For relatively low values of $B$, the influence of the parabolic potential is small and the wave function is essentially determined by edges of the ribbon, namely, the potential in transversal direction resembles an infinite potential well. However, when eigenvalue energy of the $n_{e(h)}$-th electron (hole) state in bulk phosphorene is smaller (larger) than the effective potential at the ribbon edge closer to the potential extrema, i.e. when

[TABLE]

the confinement along the $y^{\prime}$ direction is dominantly QHO-like. We should note that this criteria is not so rigid. Namely, the effective potential at the closer edge should be sufficiently larger than the QHO eigenenergy so that the corresponding eigenfunction decreases sufficiently before it reaches the ribbon edge. More comprehensive criteria are established in Ref.  [topalovic16, ], but since we have already made a few approximations it would not improve the analytical results significantly.

The smallest value of the magnetic field for LL formation in the $n$-th electron(hole) state is found when equality is replaced by $<$ in Eq. (36) and $k_{x^{\prime}}=0$

[TABLE]

Note that the minimal field $B_{n_{e(h)}}^{\min}$ is proportional to the number of the LL. In the above equation only the transverse effective mass $m^{*}_{e(h),y^{\prime}}$ depends on the edge orientation, and $B_{n_{e(h)}}^{\min}$ monotonically decreases with it. In phosphorene the largest value of the effective mass is in the zigzag direction, while the smallest one is in the armchair direction. Therefore, the magnetic field required for the formation of the LLs is the smallest in the case of AC ribbons and the largest for ZZ ribbons. Similar conclusion could be drawn intuitively from the analysis of the effective potential, based on the argument related to the confinement strength.

Furthermore, when $B>B^{\min}_{n_{e(h)}}$ we deduce from Eq. (36) that in the range $k_{x^{\prime}}\in\left(-k_{n_{e(h)},x^{\prime}}^{FB},k_{n_{e(h)},x^{\prime}}^{FB}\right)$ eigenvalues become independent on $k_{x^{\prime}}$ where

[TABLE]

is the flat-band boundary wavenumber. Namely, the band structure of states which satisfy the condition (36) should appear flat. It is straightforward to show that the flat-band boundary wavenumber increases with transverse mass. Therefore, the flat-band range is the widest for a ribbon with AC edges and the smallest for ZZ ribbon. Furthermore, the width of the flat-band decreases with increasing LL number, which is due to the fact that higher states extend more along the width of the ribbon, so the condition that the wavefunction ”touches” the ribbon edge given by Eq. (36) is satisfied for smaller values of the longitudinal momentum.

Let us also briefly discuss the influence of an in-plane electric field $E_{y^{\prime}}$ on the LLs. Such an external electric field modifies the on-site energy. Therefore, in the continuum model we add the potential $V^{ext}(E_{y}^{\prime},y^{\prime})=eE_{y^{\prime}}y^{\prime}=eE_{y^{\prime}}(\tilde{y}+l^{2}_{B}k_{x^{\prime}})$ to the diagonal terms. It is straightforward to show that the differential equations (30) are modified so that the effective potential becomes $V_{e(h)}(k_{x^{\prime}},y^{\prime},B)=\frac{1}{2}m_{e(h),y^{\prime}}\omega_{e(h)}^{2}\left(y^{\prime}-l^{2}_{B}k_{x}^{\prime}\pm eE_{y^{\prime}}/(m_{e(h),y^{\prime}}\omega^{2}_{e(h)})\right)^{2}$ while two terms are added to the right-side of the equations $\Delta E_{c(v)}\rightarrow\Delta E_{c(v)}+eE_{y^{\prime}}k_{x^{\prime}}l^{2}_{B}\mp e^{2}E^{2}_{y^{\prime}}/(2m_{e(h),y^{\prime}}\omega^{2}_{e(h)})$. Finally, we obtain similar solutions for the eigenfunctions as for the case when only the magnetic field is applied. As expected, the wavefunctions are shifted along the ribbons width in the direction of transverse electric field for holes, and in the opposite direction for electrons. Also, eigenenergies of LLs are modified

[TABLE]

We infer that energy gap, between the states in the conduction and valence bands with the same $k_{x^{\prime}}$, is reduced by $e^{2}E^{2}_{y^{\prime}}/(m_{e(h),y^{\prime}}\omega^{2}_{e(h)})$, while formerly flat bands adopt linear dispersion ($\mathtt{\sim}eE_{y^{\prime}}k_{x^{\prime}}l^{2}_{B}$) turning the band-gap to indirect. This behavior is expected since $B$ leads to a shift of states with positive (negative) momenta to the upper (lower) side of the ribbon, and therefore electrons and holes experience opposite potential shifts. Moreover, the wavenumbers that determine the boundaries of these linear regions can be found from Eq. (36) when $k_{x^{\prime}}$ is substituted by $k_{x^{\prime}}\mp\Delta k(E_{y^{\prime}},B)$, where $\Delta k(E_{y^{\prime}},B)=eE_{y^{\prime}}/(m_{e(h),y^{\prime}}\omega^{2}_{e(h)}l^{2}_{B})=(E_{y^{\prime}}/B\hbar)\cdot m_{e(h),x}m_{e(h),y}/m_{e(h),y^{\prime}}$. One may note that these spectral shifts depend linearly on $E_{y^{\prime}}$. Moreover, based on the argument regarding the value of the effective mass in a certain direction, it is easy to conclude that these shifts are smallest in a ribbon with AC edges and are largest for the one with ZZ edges. Consequently, the band dispersions become tilted when $E_{y^{\prime}}$ is applied.

The presented continuous model does not account for the edge states. The presence and the origin of these states in ribbons with various edges were discussed in detail in Refs. [grujic16, ,ezawa14, ]. Recently proposed continuous model with proper boundary conditions results in adequately modeled edge states for ZZ ribbonssousa16 . In former approach ezawa14 the edge states are treated as quasi-flat bands and are determined as the zero-energy states in the anisotropic honeycomb lattice model.

We note that for chosen directions of the magnetic and electric fields, the wavefunctions are symmetric along $x^{\prime}$ direction of the ribbon’s unit cell. Therefore, the effective translation vector that is twice shorter $({\mathbf{d}}/2)$ can be introduced. The same argument is used for the derivation of the two-band Hamiltonian. As a consequence, the First Brillouin zone (FBZ) becomes twice wider and the energy spectrum unfolds. The calculation times become more that twice shorter.

III Nanoribbons in external perpendicular magnetic and transverse electric field

Our goal is to compare nanoribbons with equal width having different edges. Therefore, we chose the number of dimmers along the ribbon cross-section to be 61, 38, 46 and 74, so we have approximately 10 nm wide nanoribbons with AC, sZZ, ZZ and sAC edges, respectively. In order to establish conditions for LL formation in all nanoribbons, we chose the value of perpendicular applied magnetic field $B=500$ T which is rather large in order to have significant spacing between LL’s which improves its visualisation. Similar results are obtained for smaller $B$ if we increase the width $W$.

In Fig. 2 we show the band structures of AC nanoribbons: (a) in the absence of fields, (b) with applied magnetic field only, and (c) with both $B$ and $E_{y^{\prime}}$ turned on. In the absence of fields these nanoribbons are insulating, grujic16 and host a series of parabolic bands, which are effectively ”sampled” from the band structure of bulk phosphorene, as might be inferred from Fig. 2(a). When the magnetic field is turned on LLs are formed and flat-bands appear, as shown in Fig. 2(b). Quasiclassically, for states away from the edges the magnetic field enforces closed elliptical cyclotron orbits. These states in turn quantize into LLs and form the flat parts of the bands displayed in Fig. 2(b). Note that these segments get narrower at higher energies and lower magnetic fields; the reason being that the cyclotron radius enlarges, and therefore fewer orbits are uninterrupted by the edges.

As we discussed in the theory part, due to the large effective mass in the transversal ZZ direction confinement is strong, and LLs might be found almost as soon as $B$ is turned on. Therefore, a large number of LLs is supported by nanoribbons with AC edges. Dashed magenta and blue curves denote theoretically calculated edges of the flat-bands in the conduction and valence band, respectively. There is good agreement between these edges and our numerical results. One of the reasons for certain small differences is explained in Sec. II, and is related to the flexible condition for these edges, given in Eq. (36). The other reasons are finite width of the ribbon and interband coupling, zhou15 due to which the spacing between LLs decreases as we move from the top (bottom) of the valence (conduction) band. Therefore, $\Delta E^{LL}_{c(v)}$ are overestimated and slopes of the theoretically calculated edges are higher than the actual ones. Furthermore, we must bear in mind that the continuum approximation is valid in a relatively narrow range of momenta around the $\Gamma$-point. Consequently, the discrepancy between the analytical and numerical results is larger for the AC ribbon.

When an electric field $E_{y^{\prime}}$ is applied the flat bands become linearly dependent on $B$ and appear as tilted. As predicted by theory, due to the linear term the bandgap switches to an indirect one, as soon as the in-plane electric field is applied. For the value of the electric field $E_{y^{\prime}}=20{\rm{mV/\AA}}$ the band-gap closes, as might be observed from Fig. 2(c). Furthermore, shifts of these linear segments in momentum space seem to be almost equal for the conduction and valence band. Earlier, we analyzed the spectral shift with respect to the effective masses and found that $\Delta k\sim m_{e(h),x}m_{e(h),y}/m_{e(h),y^{\prime}}$. In the case of AC ribbons $\Delta k^{AC}\sim m_{e(h),y}$, where $m_{e,y}=0.195m_{0}$ and $m_{h,y}=0.165m_{0}$, thus the spectral shifts have similar values for electrons and holes. We note that only for AC ribbons the shift is slightly larger for electrons than for holes. For all the other ribbons discussed in the paper, the mass-dependent part of the shift for the holes is approximately 1.5-3 times larger than for the electrons. Therefore, we expect that for ribbons other than AC, the holes are much more sensitive to the in-plane field than the electrons.

In Fig. 3 we show the real-space probability density of the states at the top (bottom) of the valance (conduction) band for $k_{x^{\prime}}=0$. The states are marked with correspondingly colored circles in Fig. 2. Numerical results are compared to the analytical solutions. In order to compare probability densities found by the tight-binding model to those determined by means of the continuum approximation, we divide the probability for occupying a certain state by the distance between the adjacent atoms in each sublattice that is $\Delta y^{\prime}_{AC}=a_{r}/2$ for AC and $\Delta y^{\prime}=\frac{c_{r}}{2}\cos\varphi=\frac{|a_{t}|}{d}\frac{a_{r}c_{r}}{2}$ for all the other. When magnetic and electric fields are not present, confinement along transversal direction is as in an infinite potential well. Wavefunctions that are obtained numerically (denoted by open circles in Fig. 3) are compared to the theoretical expression

[TABLE]

displayed by solid lines in Figs. 3(a) and (b), for the hole and the electron, respectively. We found that the wavefunction is spread along the ribbons with similar distribution for both, the electron and the hole, as predicted by the theoretical expression. According to expression (31) for the case when both magnetic and electric fields are applied, the ground state wavefunction at $k_{x^{\prime}}=0$ reads

[TABLE]

When only the magnetic field is present we found excellent agreement between our numerical results and the theoretically predicted Gaussian function from our simplified model. Namely, almost perfect bell-shaped functions are formed in the middle of the ribbon, as may be inferred from Figs. 3(c) and (d) for the hole and electron, respectively. Good agreement is also found in the presence of electric field as shown in Figs. 3(d) and (e). The wavefunctions remain Gaussian-like but shift toward upper (lower) edge for the hole (electron).

Next, we display the evolution of the band structure of sZZ nanoribbons in magnetic and electric fields in Fig. 4. In the absence of the fields, similarly for AC, these ribbons are insulating, as inferred from Fig. 4(a). When sufficiently large $B$ is applied, LLs are formed in the middle of the FBZ, i.e. the parabolic bands evolve, and start featuring dispersionless segments as we observed earlier in Fig. 4(b). By inspection of Fig. 4(b) we note that there is a smaller number of flat segments that are somewhat narrower than in the case of the AC ribbon with approximately the same width.

In order to elucidate these effects, in Fig. 5 we plot the real-space probability density of the ground states marked with correspondingly colored circles in Fig. 4, where solid and dashed curves depict $|\Psi|^{2}$ on opposite sublattices. Both sublattices exhibit the same functional variation of $|\Psi|^{2}$, with a slight lateral offset, which is a general property, independent of the applied fields. The left (right) panel column corresponds to the valence (conduction) band, while the first, second, third and fourth row in the panel correspond to $k_{x^{\prime}}=-0.8\pi/d$, $k_{x^{\prime}}=-0.2\pi/d$, $k_{x^{\prime}}=0.2\pi/d$, and $k_{x^{\prime}}=0.8\pi/d$, respectively. In the absence of fields the ground electron and hole states fully extend across the ribbon width, regardless of the $k_{x^{\prime}}$, as might be observed from the light blue curves in Figs. 5(a-h). For nonzero $k_{x^{\prime}}$ energy increases as described by $\hbar^{2}k^{2}_{x^{\prime}}/2m_{e(h),x^{\prime}}$, while confinement along ribbon width remains unchanged. The probability density of the ground states corresponding to the dispersionless segments in the valence (conduction) band are shown by the red curves in Figs. 5(c-f). We infer that these probability densities have almost unchangeable Gaussian shape for any $k_{x^{\prime}}$ in the flat-band range, as can be concluded by comparing the red curves in Figs. 5(c) and (e), for the valence band, as well as from Figs. 5(d) and (f), for the conduction band.

Next, we analyze effects of applied transverse electric field. By comparing the green curves in Figs. 5(b-g) to the red ones we found that the states are almost the same as in the case when only the magnetic field is applied. As expected, calculated wavefunctions in linear segments of the spectrum shown in Fig. 4(c) have the shape that is almost identical to the corresponding eigenfunctions in the flat-bands. However, these states are shifted along the ribbon width due to the electric field. For negative values of the longitudinal momenta where the dispersion is linear (see linear segments limited by dashed blue line in Fig. 4(c)), the hole wavefunction is localized around the center of the ribbon. Note that the valence band ground state at $k_{x^{\prime}}=0.2\pi/d$ is not in the linear region, as shown in Fig. 4(c). In fact, since $\Delta k^{sZZ}_{h}/\Delta k^{sZZ}_{e}=2.352$ linear segments are much more shifted in the valence than in the conduction band, as displayed in Fig. 4(c). Therefore, the probability density displayed in Fig. 5(e) interferes with the upper ribbon edge, and differs from the bell-shaped function shown in Fig. 5(c).

However, the most intriguing behavior occurs near the left zone edge, where crossing is found for the hole ground state. The crossing involves the states localized at the opposite edges of the ribbon (see Fig. 4(c)). Therefore, the localization of the hole ground state can be abruptly changed by the electric field, as shown in Fig. 5(a). The same effect occurs in the conduction band, but near the upper edge of the FBZ (see Fig. 4(c)). Due to the crossing the electron localization switches to the opposite edge, as shown by the green curves in Fig. 5(h). Additionally, the position of the crossings can be tuned by magnetic field.

The band structures of zigzag nanoribbons are shown in Fig. 6. As discussed in the theoretical part, the interband coupling has a smaller influence on the motion along the $\Gamma$-X direction and the dispersion is almost perfectly parabolic in most of the bands (see Fig. 6(a)). Since ZZ ribbons have weakest confinement in transversal direction even when a high $B$ is applied, only a few LLs evolve (see small regions encircled by dashed lines in Fig. 6(b)). When the electric field is applied, bands become immediately tilted, and newly formed linear bands are shifted. By comparing these shifts in the conduction and valence bands in Fig. 6(c), we realize that the hole is much more sensitive to the electric field, which is supported by the fact that the ratio $\Delta k^{ZZ}_{h}/\Delta k^{ZZ}_{e}=2.795$ for ZZ nanoribbons is the highest out of all considered ribbons.

In Fig. 7, we show the band structure of skewed armchair nanoribbons. Note that there are less flat segments in the spectrum of Fig. 7(b) corresponding to LLs forming around the middle of the ribbon than in the case of AC and sZZ ribbons. Also, these segments are much narrower, which can be explained by having in mind that the width of FBZ for sAC is 1.5-3 times narrower than for the other ribbon types. Since $\Delta k^{sAC}_{h}/\Delta k^{sAC}_{e}=2.734$ linear segments are much more shifted in the valence than in the conduction band, and dashed blue line delimiting linear bands starts at the band edge, as shown in Fig. 7(c).

Finally, we investigate the behavior of edge states. These states do not undergo Landau quantization, as can be seen in Fig. 6(b) and Fig. 7(b), which is a consequence of their exponential localization near the edges.grujic16 Instead, each QFB is effectively shifted along the $k_{x^{\prime}}$ axis, but in opposite directions. This behavior was absent in the case of edge states in zigzag graphene ribbons, since the corresponding bands were not fully detached from the bulk bands,ezawa14 as is the case with edge states in phosphorene nanoribbons.

In Fig. 8 we take a closer look at QFBs in (a) $N_{ZZ}=46$ and (b) $N_{sAC}=74$ under a range of magnetic field values. In the absence of magnetic field, the bands are degenerate for both, ZZ and sAC, as depicted by solid blue curves in Fig. 8. It is clear that the QFBs get progressively shifted in opposite directions with increasing magnetic field (see red dashed lines in Fig. 8). An explanation for this behavior is very similar to the one given in Ref. [grujic15, ]; exponentially localized states ”sample” the vector potential only in a small area where it is effectively constant, and therefore the effect of magnetic fields amounts only to a phase shift. Additionally, the two QFB states are localized at opposite edges (at $y^{\prime}=\pm W/2$, for a ribbon of width $W$), where the local vector potential $\mathbf{A}=B\left(-y,0,0\right)$ has different signs, which explains the opposite shifts.

A quantitative description can be obtained by employing the principle of minimal coupling $k_{x^{\prime}}\rightarrow k_{x^{\prime}}+\frac{e}{\hbar}A_{x^{\prime}}$. In other words, the bands get shifted left and right by $\frac{eBW}{2\hbar}$, depending on which edge they are localized at. We have confirmed that this is indeed a very good approximation for smaller fields. This also suggests that when $\pi/d=eB_{cr}W/2\hbar$, where $d$ is the unit cell length, QFBs get completely ”out of phase”, so that maxima and minima interchange places (see green dotted curves in Fig. 8). This critical magnetic field is then given as

[TABLE]

which in fact means that the magnetic flux through the unit cell ($\Phi_{cr}=B_{cr}dW$) is equal to one flux quantum ($\Phi_{0}=h/e$).

The exact expressions for the critical magnetic field in ZZ and sAC ribbons reads

[TABLE]

respectively. Here, $P_{hex}=a_{r}c_{r}/2$ is the area of plaquet projection into the $x,y$ plane. The green dotted curves in Fig. 8 show the band structure at $B_{cr}$. However, these two bands are not degenerate and thus not exactly in opposition with respect to the solid blue curves. There are two reasons for this, both having to do with the fact that edge states have some finite spread towards the ribbon center. On the one hand, this means that the edge states effectively experience somewhat smaller ribbon widths (thus increasing true $B_{cr}$). On the other hand, for larger magnetic fields the vector potential has a stronger spatial variation, so that the simple picture of phase shifts (assuming relatively constant $\mathbf{A}$ in the narrow space occupied by the edge state) gradually losses validity.

IV Summary and Conclusions

In summary, we derived a new Hamiltonian to describe the electronic structure of phosphorene nanoribbons with arbitrary edges in transverse electric and perpendicular magnetic fields. We found that when a magnetic field is turned on the states of positive and negative momenta split to the opposite sides of the nanoribbon. An analytical expression for the minimal magnetic field when the bands become flat is obtained. The boundaries of these flat segments in momentum space are also described by analytical functions. We show that both the minimal field and the extension of the flat bands depend on the type of nanoribbon edge. Furthermore, when magnetic field increases, the Landau levels spectrum emerges first in an AC nanoribbon and last in the ZZ nanoribbon of equal width. Thus, the transversal confinement of electrons is found to be weakest in the ZZ nanoribbons and strongest in the AC nanoribbons. Moreover, an application of in-plane electric field causes the band gap to decrease and turns the previously flat segments into ranges of linear variation of the energy spectra in momentum space. We found that the electric field gives rise to shifts of crossings toward the center of the phosphorene Brillouin zone. For all ribbon types we found good agreement between the numerical results obtained by means of the tight binding and the results derived from an analytical model. Finally, the analytical expression for the critical magnetic when the edge states are localized at opposite edges acquire counter phases is derived for the cases of ZZ and sAC nanoribbons.

Acknowledgements.

This work was supported by the Serbian Ministry of Education, Science and Technological Development, and the Flemish Science Foundation (FWO-Vl).