Effect of long-range structural corrugations on magnetotransport properties of phosphorene in tilted magnetic field

TL;DR

This study investigates how long-range structural corrugations influence the electronic and magnetotransport properties of phosphorene under tilted magnetic fields, revealing significant deviations in Landau levels and resistivities, especially at large tilt angles.

Contribution

It provides analytical insights into the effects of surface corrugations on phosphorene's Landau levels and magnetotransport in tilted magnetic fields, a novel aspect in 2D material research.

Findings

Landau levels deviate from linear field dependence due to corrugations

Surface corrugations significantly affect electronic spectrum at large tilt angles

Nonlinear magnetotransport effects are observable at moderate magnetic fields

Abstract

Rippling is an inherent quality of two-dimensional materials playing an important role in determining their properties. Here, we study the effect of structural corrugations on the electronic and transport properties of monolayer black phosphorus (phosphorene) in the presence of tilted magnetic field. We follow a perturbative approach to obtain analytical corrections to the spectrum of Landau levels induced by a long-wavelength corrugation potential. We show that surface corrugations have a non-negligible effect on the electronic spectrum of phosphorene in tilted magnetic field. Particularly, the Landau levels are shown to exhibit deviations from the linear field dependence. The observed effect become especially pronounced at large tilt angles and corrugation amplitudes. Magnetotransport properties are further examined in the low temperature regime taking into account impurity scattering. We calculate magnetic field dependence of the longitudinal and Hall resistivities and find that the nonlinear effects reflecting the corrugation might be observed even in moderate fields (\mbox{$B<10$ T}).

Equations

$$\mathcal{H}_{0}=\left[\begin{array}[]{cc}u_{0}+\bar{\eta}_{x}\pi_{x}^{2}+\bar{\eta}_{y}\pi_{y}^{2}&\delta+\bar{\gamma}_{x}\pi_{x}^{2}+\bar{\gamma}_{y}\pi_{y}^{2}+i\bar{\chi}\pi_{y}\\ \delta+\bar{\gamma}_{x}\pi_{x}^{2}+\bar{\gamma}_{y}\pi_{y}^{2}-i\bar{\chi}\pi_{y}&u_{0}+\bar{\eta}_{x}\pi_{x}^{2}+\bar{\eta}_{y}\pi_{y}^{2}\end{array}\right].$$

$$\mathcal{H}_{0}^{\mathrm{quad}}=\left[\begin{array}[]{cc}u_{0}+\bar{\eta}_{x}\pi_{x}^{2}+\bar{\eta}_{y}\pi_{y}^{2}&\bar{\gamma}_{x}\pi_{x}^{2}+\bar{\gamma}_{y}\pi_{y}^{2}\\ \bar{\gamma}_{x}\pi_{x}^{2}+\bar{\gamma}_{y}\pi_{y}^{2}&u_{0}+\bar{\eta}_{x}\pi_{x}^{2}+\bar{\eta}_{y}\pi_{y}^{2}\end{array}\right],$$

$$\mathcal{H}_{0}^{\mathrm{lin}}=\left[\begin{array}[]{cc}0&\delta+i\bar{\chi}\pi_{y}\\ \delta-i\bar{\chi}\pi_{y}&0\end{array}\right].$$

$$p_i$$

$$r_i$$

$$E_n^{\mathrm{quad}}$$

$$⟨n_x{n}_y|{\pi }_x^2|n_x{n}_y⟩$$

$$⟨n_x{n}_y|{\pi }_y^2|n_x{n}_y⟩$$

$$(E_n^{\mathrm{lin}}{)}^2$$

$$E_{n\lambda }^0$$

$$\mathbf{B}=\left[\frac{B_{\parallel }+B_{\perp }\sin \theta }{\sqrt{2}},\frac{B_{\parallel }+B_{\perp }\sin \theta }{\sqrt{2}},B_{\perp }\cos \theta \right],$$

$$\mathbf{A}$$

$${\pi }_x^2$$

$${\pi }_y^2$$

$$E_{n\lambda }$$

$$\Delta E_n$$

$$\Delta E_n^i$$

$${\mathbb{G}}_n$$

$$⟨\mathbf{r}|n_x{n}_y⟩$$

$${\mathbb{G}}_n$$

$$D(E)$$

$${\sigma }_{xy}$$

$${\sigma }_{xx}$$

$$E_{n\lambda }^0$$

$${\hat{n}}_x|n_x{n}_y⟩$$

$${\hat{n}}_y|n_x{n}_y⟩$$

$$E_{n\lambda }^0$$

$${\mathbb{G}}_n$$

$${\int }_0^{\infty }d\alpha \exp \left(-{\alpha }^2\right){\left[H_n\left(\alpha \right)\right]}^2\cos \left(\sqrt{2}\beta \alpha \right)$$

$$t_{{\mathbf{r}}_1,{\mathbf{r}}_2}\to t_{{\mathbf{r}}_1,{\mathbf{r}}_2}\exp \left(ie/\hslash {\int }_{{\mathbf{r}}_1}^{{\mathbf{r}}_2}\mathbf{A}\cdot d\mathbf{l}\right),$$

Full text

I Introduction

After the first synthesis of graphene Novoselov et al. (2004), the interest in two-dimensional (2D) materials has grown considerably over the past decade. The gapless energy spectrum of graphene has stimulated the search for 2D semiconductors, more suitable for traditional electronic and optoelectronic applications. Besides the other group IV materials Vogt et al. (2012); Dávila et al. (2014); Bampoulis et al. (2014); Zhang et al. (2016); Zhu et al. (2015) and a variety of transition-metal dichalcogenidesWang et al. (2012); Jariwala et al. (2014) fabricated in recent years, new elemental materials appear in the focus of attention. In this context, few-layer black phosphorus is one of the most promising 2D materials potentially interesting for practical applications Ling et al. (2015); Castellanos-Gomez (2015); Carvalho et al. (2016) because of its relatively high carrier mobilityKoenig et al. (2014); Liu et al. (2014); Li et al. (2014), tunable energy gap of 0.3–2.0 eV Tran et al. (2014); Rudenko et al. (2015); Qiao et al. (2014), and intrinsic anisotropy Qiao et al. (2014); Xia et al. (2014); Wang et al. (2015) resulting in, for instance, unusual optical response Low et al. (2014); Nemilentsau et al. (2016). Compared to graphene, properties of black phosphorus are considerably less studied theoretically, which hinders the understanding of experimentally observable phenomena.

2D materials are known to be intrinsically unstable with respect to long-wavelength thermal fluctuations, resulting in the formation of a corrugated or rippled structure in accordance with the Mermin-Wagner theorem Mermin (1968). Earlier studies demonstrated that rippling is an intrinsic feature of graphene, which affects its electronic properties Meyer et al. (2007); Fasolino et al. (2007); de Juan et al. (2007); Isacsson et al. (2008); Guinea et al. (2008a, b); Cortijo and Vozmediano (2009); Costamagna et al. (2010). Other 2D structures were also shown to have a tendency to form ripples, such as, for example, in hexagonal boron nitride Meyer et al. (2009), transition metal dichalcogenides Brivio et al. (2011); Miro et al. (2013), and black phosphorus Kou et al. (2015); Zhou et al. (2016); Wang et al. (2016). Although 2D materials are usually deposited on substrates, which may suppress the formation of intrinsic rippling, surface roughness of common dielectrics like SiO2 represents by itself another source of structural corrugations Ishigami et al. (2007); Geringer et al. (2009); Lui et al. (2009).

Understanding the dynamics of charge carriers in 2D materials under realistic conditions is a problem of practical importance as it determines observable transport properties. Magnetotransport measurements offer a powerful tool to probe carrier dynamics at the quantum level. Recently, several studies have reported quantum transport measurements in few-layer black phosphorus Chen et al. (2015); Li et al. (2015); Tayari et al. (2015); Gillgren et al. (2015); Li et al. (2016); Long et al. (2016); Tran et al. (2017); Long et al. (2017). Interpretation of experimental observations is usually carried out on a phenomenological level without explicit consideration of their microscopic nature. On the other hand, theoretical description of quantum transport at the level of model HamiltoniansYuan et al. (2015); Zhou et al. (2015a); Tahir et al. (2015); Jiang et al. (2015); Zhou et al. (2015b); Yuan et al. (2016) has limited capability to capture essential environmental effects caused by impurities, substrates, and structural corrugations. The role of those effects in magnetotransport properties of few-layer black phosphorus is not well understood.

In this paper, we study the role of long-range structural corrugations on the Landau levels (LLs) and magnetotransport properties of monolayer black phosphorus (MBP) in the presence of a tilted magnetic field. We use a perturbative approach to obtain first-order corrections to the energy spectrum induced by a corrugation potential. We find noticeable deviations of LLs from the linear dependence on magnetic field, which are also apparent in the calculated longitudinal and Hall resistivities at not very strong fields.

The paper is organized as follows. The theory part is presented in Sec. II, where we first consider unperturbed Hamiltonian for MBP in perpendicular magnetic field (Sec. II A), and then obtain a correction to the Hamiltonian in the presence of a corrugation potential in tilted magnetic field (Sec. II B). In Sec. II C, we present the formalism of the linear response theory, which is used to calculate magnetotransport properties of MBP. The results and their discussion are presented in Sec. III. In Sec. IV, we briefly summarize our results and conclude the paper.

II Theory

II.1 Pristine MBP in perpendicular magnetic field

The energy spectrum of MBP can be described by a four-band tight-binding model Rudenko and Katsnelson (2014). However, $C_{2h}$ group invariance of the MBP lattice enables one to describe the system by a two-band model Rodin et al. (2014). An effective continuum model can be obtained by expanding the tight-binding Hamiltonian around the $\mathrm{\Gamma}$ point, yielding a good agreement with the tight-binding results in the energy range of $\sim$3.5 eV Ezawa (2014); Pereira and Katsnelson (2015). In the long-wavelength limit, the continuum Hamiltonian of MBP can be written as Pereira and Katsnelson (2015),

[TABLE]

Eq. (1) can be divided into quadratic and linear terms, respectively, as

[TABLE]

[TABLE]

Here, $\bar{\eta}_{i}=\eta_{i}/\hbar^{2}$ ($\eta_{x}=0.58$ eVÅ2 and $\eta_{y}=1.01$ eVÅ2), $\bar{\gamma}_{i}=\gamma_{i}/\hbar^{2}$ ($\gamma_{x}=3.93$ eVÅ2 and $\gamma_{y}=3.83$ eVÅ2), $\bar{\chi}=\chi/\hbar$ ($\chi=5.25$ eVÅ), $u_{0}=-0.42$ eV and $\delta=0.76$ eV Pereira and Katsnelson (2015), and $\pi_{i}$ is the 2D canonical momentum. If magnetic field is applied normal to the MBP plane, B=(0,0,B), in symmetric gauge $\pi_{x}=p_{x}-(eB/2)y$ and $\pi_{y}=p_{y}+(eB/2)x$, where $p_{i}$ is the momentum operator. One can express $p_{i}$ and $r_{i}$ in terms of the creation $b_{i}^{{\dagger}}$ and annihilation $b_{i}$ operators as

[TABLE]

where $\lambda$ is the band index taking $+$1 ($-$1) for the conduction (valence) band, and $i$ refers to $x$ or $y$. $\omega_{\lambda}=eB/\sqrt{m^{\lambda}_{x}m^{\lambda}_{y}}$ is the cyclotron frequency, which takes $\omega_{+}=2.668\omega_{e}$ ($\omega_{-}=2.195\omega_{e}$) for electrons (holes) with $\omega_{e}=eB/m_{0}$, and $m^{\lambda}_{i}$ are the effective masses: $m^{+}_{x}=\hbar^{2}/2(\eta_{x}+\gamma_{x})=0.846m_{0}$ and $m^{+}_{y}=\hbar^{2}/2(\eta_{y}+\gamma_{y}+\chi^{2}/2\delta)=0.166m_{0}$ for the conduction band, $m^{-}_{x}=\hbar^{2}/2(\gamma_{x}-\eta_{x})=1.140m_{0}$ and $m^{-}_{y}=\hbar^{2}/2(\eta_{y}-\gamma_{y}-\chi^{2}/2\delta)=0.182m_{0}$ for the valance band, with $m_{0}$ being the free electron mass. By diagonalizing $\mathcal{H}_{0}^{\mathrm{quad}}$, eigenvalues of the quadratic Hamiltonian become

[TABLE]

where

[TABLE]

In turn, diagonalization of the linear term $\mathcal{H}_{0}^{\mathrm{lin}}$ yields

[TABLE]

Due to the gauge independent degeneracy of LLs, we assume $n_{x}=n_{y}=n$, thus eigenvalues of the total Hamiltonian $\mathcal{H}_{0}$ can be written as (see Appendix A for more details)

[TABLE]

The last term can be expanded as $\delta\left[1+(\bar{\chi}^{2}m^{\lambda}_{y}/\delta^{2})\hbar\omega_{\lambda}(n+1/2)\right]^{1/2}\approx\left[\delta+(\bar{\chi}^{2}m^{\lambda}_{y}/2\delta)\hbar\omega_{\lambda}(n+1/2)\right]$. Finally, we arrive at

[TABLE]

The expression given by Eq. (8) is fully consistent with the results of previous studies Pereira and Katsnelson (2015); Zhou et al. (2015a). It is clear from Fig. 1 that for $B<10$ T, two expressions in Eqs. (7) and (8) match with each other demonstrating that the linear term ($\chi$) in the continuum Hamiltonian is less effective on LLs of MBP. From Fig. 1 one can also see that both spectra are very close to the results of tight-binding calculations performed in Appendix B. In Appendix B, we also consider the case of in-plane magnetic field, which is shown to have a negligible effect on the properties of pristine (noncorrugated) MBP.

II.2 Corrugated MBP in tilted magnetic field

We now consider a vector potential that produces a tilted magnetic field,

[TABLE]

which consists of a constant field $B_{\parallel}$ along the $xy$-plane and a constant field $B_{\perp}$ tilted with respect to the $z$-axis by angle $\theta$. Modified symmetric gauge which yields this magnetic field can be chosen as

[TABLE]

Similar gauge choices were considered before for parabolic quantum wells Haupt and Wendler (1994a, b); Hai and Peeters (1999) and other 2D materials Kandemir (2010); Mogulkoc et al. (2016). In Eq. (10), even if the tilt angle $\theta$ is set to zero, the parallel component $B_{\parallel}$ still exists which, allows us to examine the effect of $B_{\parallel}$ on the energy spectrum of MBP. In the presence of tilted magnetic field, the square of the momentum operators is given by

[TABLE]

In Eq. (11), we have introduced a corrugation potential along the $xy$-plane having the form $z(x,y)=V\cos(\textrm{K}x)\cos(\textrm{K}^{\prime}y)$ which can be considered as a small perturbation on the surface of MBP (see Fig.2). Here, $\textrm{K}=2\pi/\ell_{x}$ and $\textrm{K}^{\prime}=2\pi/\ell_{y}$, $\ell_{x}$ and $\ell_{y}$ are the length of the corrugation along the $x$ and $y$ directions, respectively. $V$ is the height (amplitude) of the corrugation, $\Theta(B,\theta)=B_{\perp}(\sin\theta+\xi)$, and $\xi=B_{\parallel}/B_{\perp}$. In what follows, the effect of corrugation potential on LLs is treated perturbatively, and assuming $B_{\perp}>B_{\parallel}$ ($\xi<1$), which preserves $C_{2h}$ group invariance of the MBP lattice. More sophisticated analysis can be, in principle, performed following the variational technique Kandemir and Mogulkoc (2010a, b). Considering the modified momentum operators in Eq. (11) and following the same procedure outlined in Sec. II A, the energy eigenvalues of the system can be written as

[TABLE]

Here, $\overline{E}_{n\lambda}^{0}$ is a modified angle-dependent version of the energy eigenvalues that appeared in Eq. (8), i.e., $\overline{E}_{n\lambda}^{0}=(u_{0}+\lambda\delta)+\lambda\hbar\overline{\omega}_{\lambda}(n+1/2)$, where $\overline{\omega}_{\lambda}=\omega_{\lambda}\cos\theta$ is the modified cyclotron frequency. $\Delta E_{n}^{x}$ and $\Delta E_{n}^{y}$ are the first-order corrections to the energy eigenvalues given by

[TABLE]

with

[TABLE]

being the spatial correlation function defined as

[TABLE]

In Eq. (15), $\Psi_{n_{x}n_{y}}(x,y)=\left\langle\mathbf{r}|n_{x}n_{y}\right\rangle$, and

[TABLE]

where $H_{n}$ are the Hermite polynomials (see Appendix A). After taking the integrals in Eq. (15) and considering the assumption $n_{x}=n_{y}=n$, we get (see Appendix for more details)

[TABLE]

where $L_{n}$ are the Laguerre polynomials. To see the oscillatory nature of the Laguerre polynomials, their asymptotic expression can be used, $e^{u/2}L_{n}(u)\approx\left(\pi^{2}nu\right)^{-1/4}\cos\left(2\sqrt{nu}-\pi/4\right)$. For large $n$, $n\rightarrow(E_{F}/\hbar\omega_{\lambda})-1/2$ Vasilopoulos and Peeters (1989); Peeters and Vasilopoulos (1993) and

[TABLE]

Here, $m_{\lambda}=(m_{x}^{\lambda}m_{y}^{\lambda})^{1/2}$ is the cyclotron mass, and $n_{c}=m_{\lambda}E_{F}/\pi\hbar^{2}$ is the carrier concentration.

The density of states (DOS) for quantized energy spectrum can be calculated as

[TABLE]

where $S$ is the area of the MBP unit cell. To calculate DOS, we use the Gaussian functions as an approximation to the Dirac function in Eq. (19), i.e., $\delta\left(E-E_{n\lambda}\right)\approx\left(1/\sigma\sqrt{\pi}\right)\exp\left[-\left(E-E_{n\lambda}\right)^{2}/\sigma^{2}\right]$, with $\sigma$ being the broadening parameter taken to be $\sigma=0.1$ meV.

II.3 Magnetotransport properties

To examine the effect of tilted magnetic field on magnetotransport properties of MBP, we make use of the linear response theory. We consider a strongly quantized regime, in which $\overline{\omega}_{\lambda}\gg\tau^{-1}$, where $\tau$ is the carrier relaxation time. In the presence of the perturbative term in Eq. (12), the carrier velocity along both $x$ and $y$ directions remains zero due to the Landau quantization. In this situation, one can distinguish between the two main contributions to the conductivity tensor, namely, transverse (Hall) $\sigma_{xy}$ and longitudinal (collisional) conductivity $\sigma_{xx}$. The Hall conductivity can be readily evaluated as

[TABLE]

which is a standard expression for conventional 2D electron gas Vasilopoulos (1985); Vasilopoulos and Peeters (1989); Peeters and Vasilopoulos (1992, 1993). Here, $g_{s}=2$ stands for the spin degrees of freedom, and $f\left(E_{n,\lambda}\right)=\left[1+\exp\beta\left(E_{n,\lambda}-E_{F}\right)\right]^{-1}$ is the Fermi-Dirac distribution function, where $E_{F}$ is the Fermi energy, $\beta=1/k_{B}T$ is the inverse temperature in energy units with $k_{B}$ being the Boltzmann constant. It is worth noting that $\sigma_{xy}$ is scattering independent.

The second contribution to the conductivity tensor, i.e., longitudinal conductivity, can be evaluated as Vasilopoulos (1985); Vasilopoulos and Peeters (1989); Peeters and Vasilopoulos (1992, 1993); Vasilopoulos and Van Vliet (1986)

[TABLE]

Here, we assume scattering on randomly distributed Coulomb impurities with density $n_{\mathrm{imp}}$. Other scattering mechanism, such as scattering on phonons can be neglected in the limit of low temperatures. In Eq. (21), $\overline{\ell}_{B}=\sqrt{\hbar/eB_{\perp}\cos\theta}$ is the magnetic length, and $U^{\lambda}=2\pi e^{2}k_{e}/\epsilon k_{s}^{\lambda}$ is the impurity Coulomb potential for small momentum transfer $q\ll k_{s}$, where $\epsilon$ is the relative dielectric permittivity, $k_{e}=1/4\pi\epsilon_{0}$ is the Coulomb constant, and $k_{s}^{\lambda}=2\pi e^{2}D^{\lambda}_{0}$ is the screening wave vector of electrons (holes) in the Thomas-Fermi approximation Katsnelson (2012) with $D^{\lambda}_{0}=m_{\lambda}/\pi\hbar^{2}$ being DOS in the absence of magnetic field. Eq. (21) represents a collisional contribution to the conductivity $\sigma_{xx}=\sigma_{xx}^{\mathrm{col}}$, which increases with impurity concentration $n_{\mathrm{imp}}$ contrary to the diffusive contribution $\sigma_{xx}^{\mathrm{dif}}\sim 1/n_{\mathrm{imp}}$ Yuan et al. (2015); Liu et al. (2016). This is because impurity scattering in the presence of a magnetic field favors electron hoppings between quantized cyclotron orbits Peeters and Vasilopoulos (1992), thus increasing the conductivity. It is also interesting to note that as long as $\sigma_{xx}^{\mathrm{col}}$ dominates, $\sigma_{xx}$ remains isotropic. We note that in contrast to earlier studies Zhou et al. (2015a), we explicitly take into account impurity- and field-induced broadening of the LLs widthVasilopoulos and Van Vliet (1986). Using Eq. (21), the Hall and longitudinal resistivity can be calculated as $\rho_{xy}=\sigma_{xy}/S$ and $\rho_{xx}=\sigma_{xx}/S$, respectively, where $S=\sigma_{xx}\sigma_{yy}-\sigma_{xy}\sigma_{yx}\approx\sigma_{xy}^{2}$ assuming $\sigma_{xy}\gg\sigma_{xx}$ in sufficiently strong fields. In the following magnetotransport calculations, we use $T=1$ K, $n_{\mathrm{imp}}=10^{12}$ cm*-2*, and $\epsilon=1$. The latter corresponds to the case of freestanding non-doped MBP Prishchenko et al. (2017).

III Results and Discussion

We first examine evolution of LLs in MBP considering tilted magnetic field in the presence of a corrugation potential. Since we use a perturbative approach to describe the effect of angle-dependent magnetic field [Eq. (12)], the product $(B_{\perp}V)^{2}$ must not be too large to ensure validity of the approach, that is to satisfy $\Delta E_{n\lambda}\ll E_{n}$. This condition holds for $B_{\perp}<10$ T and $V<5$ Å considered in this work. Unless stated otherwise, we consider fixed ratio $\xi=B_{\parallel}/B_{\perp}=0.5$, and the corrugation length along both directions $\ell_{1}=\ell_{2}=250$ Å, which is an order of intrinsic ripples length in graphene Meyer et al. (2007); Fasolino et al. (2007). The case $\theta=0$ is evaluated for $V=0$ to be consistent with the results of earlier works Pereira and Katsnelson (2015); Zhou et al. (2015a).

In the left panel of Fig. 3, we show the Landau level diagram calculated for both electron and hole states for different values of tilt angle $\theta$, and fixed amplitude of the corrugation potential $V=1$ Å. One can see that energies of LLs decrease with $\theta$, while the linearity of the curves is preserved in the regime of relatively small corrugations and not too strong magnetic fields. The effect of the tilt angle on LLs is twofold. While $B_{\perp}$ confines the motion of charge carriers in the $xy$ plane, changing the magnetic field direction increases the cyclotron radius in the $xy$ plane due to the $\mathrm{cos}\theta$ factor in $\overline{\omega}_{\lambda}$. As a result, the LL energies $E_{n}$ decrease with $\theta$, which effectively correspond to a smaller magnetic field. This effect is partially compensated by the presence of the corrugation potential, which provides an additional contribution $\Delta E_{n}$ to $E_{n}$ [Eq. (13)]. As can be inferred from Fig. 3, the main contribution to LLs comes from the first term, $\overline{E}_{n\lambda}^{0}$ in Eq. (12), which is strongly dependent on the perpendicular component of the out-of-plane magnetic field $B_{\perp}\cos\theta$. The effect of the second term, $\Delta E_{n}$ in Eq. (12) is small for corrugations as low as $V=1$ Å. The effect of tilted magnetic field on the electronic spectrum can be seen also from DOS shown for different tilt angles at $B_{\perp}=5$ T and $B_{\perp}=10$ T (right panel of Fig. 3). At $\theta\neq 0$, the energy spacing between LLs becomes smaller, which leads to more dense electronic states and larger DOS. A similar effect of tilted magnetic field on LLs of graphene was reported previouslyKandemir (2010); Krstajić (2013). For larger values of $B_{\perp}$, the energy spacing between LLs decreases, which gives rise to more pronounced oscillations in DOS shown in Fig. 3 for $B_{\perp}=10$ T.

In Fig. 4(a), the index ($n$) dependence of LLs is shown in the presence of tilted magnetic field for $V=1$ Å. LLs splitting of electron and hole states is different because of the electron-hole asymmetry and unequal effective masses. In Fig. 4(b), the magnetic field ($B_{\perp}$) dependence of LLs is shown for different corrugation heights at $\theta=30$ Å. One can see pronounced deviations from the linear behavior, which become especially clear for $V=5$ Å and $B_{\perp}=10$ T. In the chosen range of parameters, these deviations do not exceed $\hbar\overline{\omega}_{\lambda}(n+1/2)$, demonstrating the validity of the perturbative approach. The observed nonlinearity is a manifestation of long-range structural corrugations. Although at relatively weak fields, first-order correction to the LLs energy is quadratic in $V$ [Eq. (13)], the dependence at large fields may be different. Particularly, we do not exclude oscillatory behavior in this regime.

Fermi energy as a function of magnetic field is shown in Fig. 5 for different electron concentrations $n_{c}$. For fixed Fermi energy, carrier concentration can be calculated by the formula, $n_{c}=\int_{0}^{E_{F}}D\left(E\right)dE$. Here, we see the magnetic field dependence of fixed Fermi energies for different carrier concentrations.

To gain insight into the role of other model parameters on the LL spectrum, we analyze $\Delta E_{n}$ in more detail. In Fig. 6(a), $\Delta E_{n}$ is shown both for electrons and holes as a function of the carrier concentration $n_{c}$ calculated for different tilt angles $\theta$ at $B_{\perp}=10$ T. It can be seen that $\Delta E_{n}$ exhibits oscillations with $n_{c}$, and its amplitude increases for larger $\theta$. This behavior is attributed to the $\sin\theta$ factor in $\Delta E_{n}$ [Eq. (13)]. The hole states turn out to be less affected by the magnetic field direction, which is due to the higher cyclotron mass. In Fig. 6(b), we show $\Delta E_{n}$ as a function of $\theta$ calculated for different $V$. According to Eq. (13), $\Delta E_{n}\sim V^{2}(\sin\theta+\xi)^{2}$, meaning that the nonlinear effects in the spectrum of LLs increase both with $V$ and $\theta$. At small $\theta$, $\Delta E$ raises linearly, whereas at larger $\theta$, $\Delta E_{n}$ demonstrates a quadratic behavior.

In Fig. 7(a), we show the effect of a parallel magnetic field $B_{\parallel}$ by calculating the dependence of $\Delta E_{n}$ on the dimensionless parameter $\xi=B_{\parallel}/B_{\perp}$. One can see the expected from Eq. (13) $\Delta E_{n}\sim\xi^{2}$ behavior, suggesting that at large $\xi$ the in-plane field might play a role in the energy spectrum of corrugated MBP samples. The absolute effect of $B_{\parallel}$ is, however, not large and can hardly be detected experimentally under realistic field strengths and corrugation amplitudes. A similar effect of $B_{\parallel}$ on LLs was also reported previously in the context of bilayer graphene Hyun et al. (2012). Fig. 7(b) shows the effect of the corrugation length $\ell_{1}$ as well as its anisotropy $\beta=\ell_{2}/\ell_{1}$ on $\Delta E_{n}$. In this case, $\Delta E_{n}$ exhibits a complicated oscillatory behavior. Keeping in mind anisotropic ripple formation typical to MBPKou et al. (2015), we also examine $\Delta E_{n}$ as a function of $n_{c}$ for anisotropic corrugation patterns [see inset of Fig. 7(b)]. Depending on the corrugation direction, the behavior of LLs is significantly different. One can see, however, that $\Delta E_{n}$ remains weakly affected by a particular corrugation pattern as well as by the corrugation length.

We now turn to the results of our magnetotransport calculations to see whether weak effects induced by the corrugation could be observed experimentally. The Hall ($\sigma_{xy}$) and longitudinal ($\sigma_{xx}$) conductivities are shown for different tilt angles in Fig. 8. For Fermi energies in the gap region between the valence and conduction states, dc conductivity is obviously zero due to the absence of charge carriers. Beyond the gap region, $\sigma_{xy}$ exhibits distinct plateaus, arising from the discrete nature of the LL spectrum [Fig. 8(a)]. The Hall conductivity increases by $2e^{2}/h$ for each level forming the integer Hall plateaus indexed as $0,\pm 2,\pm 4,\pm 6...$. It can be seen that $\sigma_{xy}$ increases with tilt angle, which is attributed to larger DOS caused by more dense LLs (cf. Fig. 3). For the same reason, $\sigma_{xy}$ becomes smaller in stronger fields [Fig. 8(b)]. At sufficiently small $V$, $\sigma_{xy}\sim(B_{\perp}\cos\theta)^{-1}$. The longitudinal conductivity $\sigma_{xx}$ exhibits oscillatory behavior typical to the Shubnikov-de Haas (SdH) oscillations, as shown in Figs. 8(c) and (d). $\sigma_{xx}$ also increases with $\theta$ yet more slowly than $\sigma_{xy}$, because in this case $\sigma_{xx}\sim(B_{\perp}\cos\theta)^{-1/2}$ due to a factor $\overline{\ell}_{B}$ in the denominator of Eq. (21). We note that the effect of in-plane magnetic field ($B_{\parallel}$) and the corrugation lengths ($\ell_{1}$ and $\ell_{2}$) are negligible in the context of magnetotransport properties of MBP and, therefore, not presented here. The role of the corrugation amplitude is analyzed below.

In Fig. 9(a), the Hall $\rho_{xy}$ and longitudinal $\rho_{xx}$ resistivity are shown as a function of magnetic field calculated at different $\theta$ for the case of electron doping $n_{c}\approx 3\times 10^{16}$ m*-2* and $V=1$ Å. The behavior of $\rho_{xy}$ is closely related to $\sigma_{xy}$ shown in Fig. 8. As expected, $\rho_{xy}$ increases linearly with $B_{\perp}$, while larger $\theta$ correspond to effectively weaker fields. At large fields, $\rho_{xy}$ becomes quantized increasing by the unit of $\rho_{0}=h/2e^{2}$. For a given magnetic field, the Hall plateaus observed at different $\theta$ correspond to different filling factors $j$. The filling factors increase with $\theta$, meaning that the effect of the tilt angle is opposite to that of the Fermi energy, i.e., larger $\theta$ corresponds to shifting the position of the Hall plateaus toward smaller $E_{F}$ and vice versa. Similar effect of tilted magnetic field on resistivities was reported for grapheneKrstajić (2013). As can be seen from Fig. 9(a), the behavior of $\rho_{xx}$ is also similar to $\sigma_{xx}$, exhibiting pronounced SdH oscillations as well as a $\theta$-dependence of the peak amplitudes. One of the most interesting result is presented in Fig. 9(b), which shows the effect of the corrugation amplitudes $V$ on $\rho_{xy}$ and $\rho_{xx}$ calculated for a fixed $\theta$. Although the effect of $V$ is less pronounced compared to $\theta$, it becomes clearly seen at fields $B_{\perp}>5$ T. Larger $V$ shift the Hall plateaus in $\rho_{xy}$ as well as SdH oscillation peaks in $\rho_{xx}$ toward weaker magnetic fields. For sufficiently strong fields, the difference in corrugation amplitudes of a few Å results in a notable contraction of the $\rho(B_{\perp})$ spectrum along the field axis reaching 0.5 T at $B_{\perp}\sim 8$ T. Given that a realistic corrugation pattern would be represented by a superposition of different corrugation amplitudes, we expect a broadening of the SdH peaks increasing with $B_{\perp}$ under experimental conditions. Although such a behavior is typical to experimentally measured longitudinal and Hall resistivity in few-layer BPTayari et al. (2015); Li et al. (2015, 2016); Long et al. (2016, 2017), it is usually attributed to the Zeeman spin-splitting, well described by the standard Lifshitz-Kosevich formula for 2D resistivity. To reveal the effect of corrugations in the resistivity measurements, the broadening must increase with tilt angle, which is apparently not observed in known experiments on few-layer BP. The resolution of the available experimental spectra also does not allow us to observe nonlinear effects in the SdH oscillations. We note, however, that the existing magnetotransport measurements has been performed on a few-layer BP, which should be significantly less affected by the structural corrugations compared to single-layer samples.

IV Conclusions

In summary, we studied the effects of tilted magnetic field and long-range structural corrugations on LLs and magnetotransport properties of MBP. We considered an analytical model and obtained first-order corrections to the LL energies induced by corrugations in the long-wavelength limit. The energies of LLs are found to be expectedly decreasing with the tilt angle due to the $\cos\theta$ factor in the modified cyclotron frequency. At sufficiently strong fields, however, the corrugation potential induces nonlinear deviations in the dependence of LL energies on magnetic field. We find that these deviations are predominantly affected by the corrugation amplitude, whereas the corrugation length and its specific direction are less relevant. We also examined the magnetotransport properties of MBP in the presence of corrugations under tilted magnetic field within the scheme of linear response theory. Overall, the tilt angle modifies the resistivity spectra considerably, effectively reducing the magnetic field strength. In the presence of long-range corrugations, we find that both Hall and longitudinal resistivity spectra display: (i) a shift toward weaker magnetic fields, and (ii) additional broadening of the SdH peaks increasing with magnetic field, not related to the Zeeman splitting. The obtained effects are noticeable even at moderate ($B<10$ T) fields, which allows us to expect that they might be observable experimentally for MBP samples deposited on sufficiently corrugated (e.g., SiO2) substrates.

Acknowledgements.

A. Mogulkoc would like to thank Professor B. S. Kandemir for fruitful discussions. A.N.R. acknowledges support from the Ministry of Education and Science of the Russian Federation, Project No. 3.7372.2017/BP.