Magnetic Effects on Topological Chiral Channels in Bilayer Graphene

TL;DR

This paper investigates how magnetic fields influence topological chiral channels in bilayer graphene, revealing unexpected shifts and layer imbalances, with analytical and numerical insights into these effects.

Contribution

It provides new understanding of magnetic effects on topological edge states in bilayer graphene, including analytic solutions and a combined approximation approach.

Findings

Magnetic field shifts chiral channels away from electric interfaces.

Magnetic field causes layer occupation imbalance in channels.

System behavior can be approximated as a weighted sum of electric and magnetic dominant limits.

Abstract

We study the effect of a magnetic field on topological chiral channels of bilayer graphene at electric domain walls. The persistence of chiral edge states is attributed to the difference in valley Chern number in the regions of opposite electric field. We explore the regime of large electric and magnetic fields perpendicular to the lattice. The magnetic field shifts the channel away from our electric interface in a way that is inconsistent with the semiclassical expectation from the Lorentz force. Moreover, the magnetic field causes an imbalanced layer occupation preference to the chiral channels. These behaviors admit analytic solutions in the limits that either the electric or the magnetic field dominates. We numerically show in the general case that the system can be well-approximated as a weighted sum of the two limits.

Equations

$$H=$$

$$-\frac{V}{2}\sum _{\mathbf{r},s}(-1{)}^s(a_{\mathbf{r},s}^{†}{a}_{\mathbf{r},s}+b_{\mathbf{r},s}^{†}{b}_{\mathbf{r},s})$$

$$=$$

$$H(\mathbf{k})=\left(\begin{array}{cccc}V/2& tf(\mathbf{k})& 0& u\\ t{f}^{\ast }(\mathbf{k})& V/2& 0& 0\\ 0& 0& -V/2& tf(\mathbf{k})\\ u& 0& t{f}^{\ast }(\mathbf{k})& -V/2\end{array}\right)$$

$$\Delta =\frac{V}{\sqrt{1+V^2/u^2}}\approx \{\begin{array}{ll}V,& V\ll u\\ u,& V\gg u\end{array}.$$

$$H_{\mathrm{k}\cdot \mathrm{p}}$$

$$=\left(\begin{array}{cccc}V/2& \hslash v_F\pi & 0& u\\ \hslash v_F{\pi }^{†}& V/2& 0& 0\\ 0& 0& -V/2& \hslash v_F\pi \\ u& 0& \hslash v_F{\pi }^{†}& -V/2\end{array}\right)$$

$$=\hslash v_F(\nu \delta k_x{\sigma }_x-\delta k_y{\sigma }_y)\otimes {\text{openone}}_2$$

$$+\frac{V}{2}{\text{openone}}_2\otimes {\tau }_z+\frac{u}{2}({\sigma }_x\otimes {\tau }_x-{\sigma }_y\otimes {\tau }_y)$$

$$\mathbf{k}⟶\mathbf{k}-\frac{e\mathbf{A}}{\hslash }.$$

$$H_{\mathrm{eff}}=-\frac{v_F^2}{u}\left(\begin{array}{cc}0& {\pi }^2\\ \text{quantity}({\pi }^{†}{)}^2& 0\end{array}\right)$$

$$E_n=\mathrm{sgn}(n)\frac{2{\hslash }^2{v}_F^2}{u{l}_B^2}\sqrt{\text{absolutevalue}n(\text{absolutevalue}n-1)}$$

$${\psi }_{B,1}(y)$$

$${\psi }_{B,2}(y)$$

$$\mathcal{N}$$

$$\text{expectationvalue}y=-l_B^2\delta k_x+\frac{l_B^{2}u\Re ({\beta }_1^{\ast }{\beta }_2)}{\sqrt{(\hslash v_F{)}^2+l_B^2{u}^2/2}}$$

$${\sigma }_y^2=\frac{l_B^2(\hslash v_F{)}^2+l_B^4{u}^2[1+\Re (c_1^{\ast }{c}_2)-2\text{imaginary}(c_1^{\ast }{c}_2{)}^2]}{2(\hslash v_F{)}^2+l_B^2{u}^2}$$

$$P_{\text{bottom}}=\frac{1+|{\beta }_2{|}^2-|{\beta }_1{|}^2}{2+l_B^2{u}^2/(\hslash v_F{)}^2}.$$

$$0\le P_{\text{bottom}}\le \frac{2}{2+l_B^2{u}^2/(\hslash v_F{)}^2}.$$

$$P_{\text{bottom}}=\frac{e\hslash v_F^2}{u^2}\text{quantity}(1+|{\beta }_2{|}^2-|{\beta }_1{|}^2)B.$$

$$\kappa =\frac{\text{absolutevalue}V}{2\sqrt{2}\hslash v_f}\sqrt{-1+\sqrt{1+8\text{quantity}(\frac{u}{V}{)}^2}}$$

$$k=\delta k_x=\frac{\text{absolutevalue}V}{4\hslash v_F}\sqrt{1+\sqrt{1+8\text{quantity}(\frac{u}{V}{)}^2}}.$$

$${\psi }_V(y)=\{\begin{array}{ll}{\psi }_V^{-}(y),& y<0\\ {\psi }_V^{+}(y),& y\ge 0\end{array}$$

$${\psi }_V^{-}(y)=\text{cV(V/2ℏvf)eiky+i(V/2ℏvf)e−iky((1−i)k−κ)eiky+((−1+i)k−iκ)e−iky((1−i)k−iκ)eiky+((−1+i)k−κ)e−ikyi(V/2ℏvf)eiky+(V/2ℏvf)e−ikyeκy}$$

$${\psi }_V^{+}(y)=\text{cV(V/2ℏvf)eiky+i(V/2ℏvf)e−iky((−1+i)k−κ)eiky+((1−i)k−iκ)e−iky((−1+i)k−iκ)eiky+((1−i)k−κ)e−ikyi(V/2ℏvf)eiky+(V/2ℏvf)e−ikye−κy}$$

$$c_V=\frac{\sqrt{\kappa }}{4k}.$$

$${\sigma }_y=\sqrt{\frac{1}{2{\kappa }^2}+\frac{{\kappa }^2}{2(k^2+{\kappa }^2{)}^2}}$$

$${\mathrm{Ch}}_1=\frac{1}{2\pi }\sum _n\int {\text{differential}}^{2}k({\mathbf{\Omega }}_n{)}_z,$$

$${\mathbf{\Omega }}_n=i⟨\frac{\partial u_n}{\partial \mathbf{k}}|\text{crossproduct}|\frac{\partial u_n}{\partial \mathbf{k}}⟩$$

$$\mathrm{Ch}_{1}=\frac{1}{24\pi^{2}}\int\differential{\tilde{E}}\differential^{2}k\,\mathrm{Tr}\bigg{[}\epsilon^{\mu\nu\rho}G\partial_{\mu}G^{-1}G\partial_{\nu}G^{-1}G\partial_{\rho}G^{-1}\bigg{]}$$

$${\mathrm{Ch}}_1=\frac{1}{8{\pi }^2}\int \text{differential}\tilde{E}{\text{differential}}^{2}k\mathrm{Tr}\text{quantity}[G[G{\partial }_1{G}^{-1},G{\partial }_2{G}^{-1}]].$$

$${\mathrm{Ch}}_1=-\nu \mathrm{sgn}V.$$

$$3m{1}^{\prime }={\mathbb{Z}}_2^T\times C_{3v}.$$

Full text

I Introduction

In recent years, there has been considerable interest in the electronic properties of bilayer graphene McCann, E. et al. (2007); Castro Neto et al. (2009); McCann and Koshino (2013); Rozhkov et al. (2016) including a quadratic dispersion with Berry phase of $2\pi$ and an unusual quantum Hall effect Novoselov et al. (2006); Fal’Ko (2008). The inequivalence of the $K$ and $K^{\prime}$ Dirac points in the Brillouin zone leads to a valley degree of freedom, which can play a similar role as electron spins and may be manipulated in “valleytronic” devices Rycerz et al. (2007); Behnia (2012); Ang et al. (2017). New techniques have also made it possible to produce large flakes of Bernal (AB-stacked) graphene on the order of $10^{-4}$ m, Ji et al. (2017) motivating the experimental study of low-energy electron transport.

Notably, a tunable bandgap has been realized both theoretically as a result of layer-asymmetry between the on-site energies Min et al. (2007) and experimentally by applying an electric field perpendicular to the layers Zhang et al. (2009). In pure samples where inter-valley mixing is suppressed Morozov et al. (2006); McCann et al. (2006); Morpurgo and Guinea (2006), a non-local valley symmetry emerges and the insulating electronic system carries a non-trivial valley Chern number Zhang et al. (2013). This symmetry-protected topological nature Vaezi et al. (2013) of the ground state is revealed by a one-dimensional highway of electrons along an electric domain wall Martin et al. (2008); Qiao et al. (2011); Jung et al. (2011); Zarenia et al. (2011), across which the interlayer potential difference, and consequently the valley Chern number, change sign. The electronic highway consists of four counter-propagating pairs of chiral channels, where electrons are laterally confined but delocalized along the line interface. The forward and backward propagating directions are locked with the two valley index, so that all electron modes of the same valley species propagate in the same direction. Inter-valley scatterings are forbidden by the valley symmetry and lattice momentum conservation along generic line interfaces 111Except along commensurate ones, such as the armchair edge, where the valley momentum difference is projected out, ${\bf a}_{\parallel}\cdot(K-K^{\prime})\in 2\pi\mathbb{Z}$ for the primitive parallel lattice vector ${\bf a}_{\parallel}$.. Moreover, they are chiral anomalous and associate a non-conservative valley current under an interface-parallel electric field. The anomaly is resolved by connecting the high-energy modes to a higher dimensional topological bulk, in this case sandwiched between two valley Chern insulators with opposite topological indices, or allowing the switching between valleys in high-energy.

These ballistic topological chiral channels have been realized and experimentally observed in bilayer graphene electric domain walls Ju et al. (2015); Li et al. (2016); Lee et al. (2017); Li et al. (2018). The valley-symmetry-protected channels carry a quantized differential conductance $\sigma_{0}=dI/dV=4e^{2}/h$ at zero bias, where the factor of 4 associates to the four counter-propagating pairs of Dirac modes along the interface. This value has been reported in Ref. Lee et al., 2017 using bilayer graphene encapsulated between atomically clean hexagonal boron nitride single crystals, which has been implemented in similar systems to reduce intervalley scatterings Dean et al. (2010); Kim et al. (2016). On the other hand, intervalley scatterings, which may be induced by local disorder from uneven substrates and gates among other factors, are non-negligible in other setups. Ref. Ju et al., 2015 reported a domain wall mean free path of $l_{0}\sim 400$ nm, and Ref. Li et al., 2016 inferred a mean free path of $l_{0}\sim 200$ nm. Both are shorter than or comparable to a typical domain wall length $L$, which can range between $200$ nm to 1 $\mu$m. This leads to a significant deviation of the differential conductance from the theoretical quantized value, $\sigma=\sigma_{0}/(1+L/l_{0})$, according to the Landauer–Büttiker formula Datta (1995). Intervalley scattering can be suppressed in these systems by a uniform perpendicular magnetic field $B$. A $4e^{2}/h$ approaching differential conductance has been observed in ref. Li et al., 2016 at around $B=8$ T for $L=0.4$ $\mu$m, and similar ballistic transport has been reported in Ref. Li et al., 2018 under similar field strength and length scale. The realization of differential conductance approaching the ballistic limit without a magnetic field in Ref. Lee et al., 2017 suggests the observed conductance cannot be attributed to the quantum Hall effect.

In addition to enhancing the robustness of the topological chiral channels, the magnetic field provides an external parameter that controls the microscopic properties of the electron modes Zarenia et al. (2011); Wang et al. (2017). In particular, we focus on the asymmetric lateral and layer distribution Li et al. (2016) of the electronic wavefunctions due to explicit symmetry breaking by the magnetic field. The wavefunction deviation from the domain wall center position accounts for some of the ballistic transport phenomena observed in a valley valve and electron beam splitter Li et al. (2018), and introduces new complexities to the critical transport behaviors of bilayer graphene domain wall quantum point contacts Wieder et al. (2015). In this paper, we concentrate on the single-body non-interacting understanding and quantitative description of a single domain wall under a large electric potential kink and a strong magnetic field. With potential future applications in micro- or nano-circuit devices in mind, we are interested in the dependence of microscopic quantities, for example the average channel position and lateral localization, on external control parameters, such as Fermi energy, electric and magnetic field. We show that a semiclassical description of these topological channels is inaccurate. For instance, the spatial deviation of the chiral channels under a magnetic field is inconsistent with the prediction using the classical Lorentz force, and can only be explained using Landau level physics.

This paper is organized as follows. In Section II, we review the tight-binding model applied to bilayer graphene in the presence of a uniform external electric field. The Bloch Hamiltonian for bilayer graphene is calculated and linearized about the Dirac points using the k $\cdot$ p approximation. Section III discusses how the Hamiltonian is modified in the presence of a strong magnetic field. Analytic solutions for the energy levels, wavefunctions, average position, localization, and layer preference are found. The emergence of chiral modes at electric domain walls is derived in Section IV. We calculate the lattice momenta, wavefunction, and localization of the electrons on the zero-energy chiral channels. A topological explanation of the chiral modes is given in Section V including an argument for why the valley Chern number is a topological invariant of our system. In Section VI, we consider the regime where both the electric and magnetic fields are important. We numerically find how the energy bands, average position, localization, and layer preference depend on the magnetic and electric fields as well as the Fermi energy. Here, we introduce an approximation scheme to understand the prominent features in the case of large electric and magnetic fields. Concluding remarks are given in Section VII.

II Tight Binding Hamiltonian

We review a tight-binding model Hamiltonian McCann and Koshino (2013) of bilayer graphene in the Bernal $AB$-stacking under a perpendicular electric field. Bilayer graphene may be described using a four atom unit cell consisting of the points $(A1,B1,A2,B2)$, describing the inequivalent $A$ and $B$ sites on each layer. The translations of this unit cell form a hexagonal Bravais lattice with the primitive vectors $\mathbf{a_{1}}=a(1,0)$ and $\mathbf{a_{2}}=a(-1/2,\sqrt{3}/2)$ with lattice constant $a=2.46$ ÅSaito et al. (1998), as in monolayer graphene. AB-stacked bilayer graphene is formed by stacking two monolayer sheets with the two layers displaced by a distance $a/\sqrt{3}$. In our tight-binding approximation, we consider interactions between nearest neighbors. We ignore the spin degree of freedom due to weak spin-orbit interactionYao et al. (2007), the direct $B1$ to $A2$ hopping, as well as next-nearest-neighbor and other weaker tunnelling processes. Using an intralayer nearest-neighbor hopping parameter $t=3.16$ eVMcCann and Koshino (2013), interlayer hopping parameter $u=0.381$ eV between the $A1$ and $B2$ sites that lie directly on top of each otherMcCann and Koshino (2013), and electric potential energy difference $V$ between the layers, the tight-binding Hamiltonian can be written as

[TABLE]

where ${\bf c}_{\bf k}=(a_{{\bf k},1},b_{{\bf k},1},a_{{\bf k},2},b_{{\bf k},2})^{T}$, $s=\{1,2\}$ sums over the layers, and ${\bf r}$ runs over Bravais lattice vectors $m_{1}{\bf a}_{1}+m_{2}{\bf a}_{2}$. The Bloch Hamiltonian is

[TABLE]

where $f(\mathbf{k})=1+e^{-i\mathbf{k}\dotproduct\mathbf{a_{2}}}+e^{-i\mathbf{k}\dotproduct(\mathbf{a_{1}}+\mathbf{a_{2}})}$ are the interlayer hopping terms in momentum space describing the three nearest-neighbor points. Diagonalizing (2) gives a quadratic dispersion with zero energy gap when $V=0$ at the inequivalent $K$ and $K^{\prime}$ points, on the edge of the Brillouin zone (Figure 2). An energy gap $\Delta$ is introduced with a non-vanishing potential

[TABLE]

If we are only concerned with low energy perturbations, we may expand $H$ about the $K$ and $K^{\prime}$ points

[TABLE]

where $\pi=\nu\delta k_{x}+i\delta k_{y}$, $v_{F}=\sqrt{3}at/2\hbar$ is the Fermi velocity for monolayer graphene, and $\nu$ is the valley index, $+1$ for $K$ and $-1$ for $K^{\prime}$. Here, $\delta k_{x}$ and $\delta k_{y}$ are small momenta deviations away from the Dirac points, and the $2\times 2$ Pauli matrices $\sigma$ ($\tau$) act on the $AB$-sublattice (layer) degree of freedom.

III Landau Levels

In this section, we discuss the other limit where the layer potential energy difference $V$ is absent and the bilayer system is under a uniform magnetic field $\mathbf{B}=B\hat{\mathbf{z}}$. The momentum in the tight-binding Hamiltonian is replaced by the canonical momentum

[TABLE]

The $\mathrm{k}\cdot\mathrm{p}$ Hamiltonian in (4) has the symmetry $[H_{\mathrm{k}\cdot\mathrm{p}}(K),Z]=[H_{\mathrm{k}\cdot\mathrm{p}}(K^{\prime}),Z]=0$ for $Z=\sigma_{z}\otimes\tau_{z}$ at the Dirac points where $\delta{\bf k}=0$. This allows the partition of our basis into the dimer $(A2,B1)$ and $(A1,B2)$ sites that correspond to the two eigenspaces $Z=\pm 1$. Focusing on one of the two eigenspaces, the energy spectrum near the Fermi level may be found by the effective two-band HamiltonianMcCann and Fal’ko (2006); Pereira et al. (2007)

[TABLE]

which is valid for $E\ll u,t$. We assume that $\psi_{A2}$ and $\psi_{B1}$ are shifted harmonic oscillator eigenstates. Through appropriate use of ladder operators, we may determine the energy spectrum

[TABLE]

where $l_{B}=\sqrt{\hbar/eB}$ is the magnetic length.

The energy spectrum near the Dirac points is independent of momentum, corresponding to flat bands. The low-energy electrons in the bulk of the material are bound to cyclotron orbits. Choosing the Coulomb gauge $\mathbf{A}=-By\hat{\mathbf{x}}$, our tight-binding Hamiltonian about the Dirac points is then given by (4) with $V=0$, $\pi=\nu(\delta k_{x}+eBy/\hbar)+\partial_{y}$ and $\pi^{\dagger}=\nu(\delta k_{x}+eBy/\hbar)-\partial_{y}$, where we replace $k_{y}\to-i\partial_{y}$ because our Hamiltonian is now $y$-dependent. We seek zero-energy states in an infinite lattice; this requires two of the four components to be zero (which components are zero depends on the value of $\nu$). For concreteness, we consider the case $\nu=+1$. Solving for the two components introduces two arbitrary constants. One is fixed by enforcing normality, while the other is freely chosen to simplify calculations. We choose the remaining constant such that one of the solutions is purely localized to the top layer. The wavefunctions (suppressing the plane wave factor $e^{ik_{x}x}$) are

[TABLE]

where the normalization constant is

[TABLE]

A general state at the lowest Landau level will then be a linear combination $\psi_{B}=\beta_{1}\psi_{B,1}+\beta_{2}\psi_{B,2}$. We observe that $\expectationvalue{y}=\bra{\psi_{B}}y\ket{\psi_{B}}$ depends linearly on the momentum $\delta k_{x}$ and is inversely proportional to the magnetic field strength.

[TABLE]

This result for $K^{\prime}$ is similar in that $\expectationvalue{y}\propto-l_{B}^{2}\delta k_{x}$.

The localization length of $\psi_{B}$, $\sigma_{y}=\sqrt{\expectationvalue{y^{2}}-\expectationvalue{y}^{2}}$, has no dependence on $\delta k_{x}$.

[TABLE]

The asymmetry in the two nonzero components of $\psi_{B}$ suggests that in a magnetic field, electrons have a greater probability of occupying one of the two layers. Since we have two independent states, the layer preference will depend on how we take a linear combination of these two states. So, we can simply calculate the probability for an electron to occupy the bottom layer, $P_{\text{bottom}}$, by integrating the squared magnitude of the fourth component of $\psi_{B}$, which is just the fourth component of $\psi_{B,2}$

[TABLE]

$P_{\text{bottom}}$ is independent of $\delta k$, so for a fixed magnetic field and only small deviations away from the Dirac points, one can predict how the electrons localize to a certain layer. Even without knowledge of $\beta_{1}$ and $\beta_{2}$, one can put a bound on $P_{\text{bottom}}$

[TABLE]

Note that the lower and upper bounds are obtained using the states $\psi_{B,1}$ and $\psi_{B,2}$, respectively. For magnetic fields of up to 21 T, $l_{B}^{2}u^{2}/(\hbar v_{F})^{2}>10$, so the layer preference is approximately linear in $B$. This is the magnetic regime used in experimental setups Li et al. (2016).

[TABLE]

At the $K^{\prime}$ point, the opposite layer preference was found.

IV Chiral Modes

Here, we summarize our findings from applying a first order $\mathrm{k}\cdot\mathrm{p}$ expansion about a Dirac point to bilayer graphene with a nonuniform electric field. Equation (4) is modified by making the replacement $V\to V\operatorname{sgn}y$. This corresponds to introducing an electric field perpendicular to the lattice that switches direction at $y=0$ along the screw axis as shown in Figure 1. We must again replace $\delta k_{y}\to-i\partial_{y}$ owing to translation symmetry breaking in $y$. Finding the eigenstates of this matrix requires solving a system of four coupled differential equations with appropriate boundary conditions at the interface. We explicitly found solutions in the zero-energy case. The solutions are proportional to $e^{i(k+i\kappa)y}$ where $\kappa$ is given by

[TABLE]

The value of $k$ is to be interpreted as the momentum where our chiral channels are located in reciprocal space.

[TABLE]

The solution takes the following form:

[TABLE]

[TABLE]

[TABLE]

where the normalization constant is

[TABLE]

This describes a total of four zero-energy modes in the Brillouin zone, at $K\pm\delta k_{x}$ and $K^{\prime}\pm\delta k_{x}$. The band structure indicates that two are right-moving modes near the $K$ point and the other two are left-moving near the $K^{\prime}$ point that connect the valence and conducting bands. These chiral modes are exponentially localized directly on the electric interface, and the localization asymptotically approaches $\hbar v_{F}/u\approx 17.7$ Å  in the $V\gg 1$ eV limit.

[TABLE]

The probability distribution falls off as $e^{-y}$, unlike the magnetic states which fall off as $e^{-y^{2}}$. We provided further confirmation using numerical results from the tight-binding model.

We considered three different potential profiles: step-function, linear, and sinusoidal, as seen in Figure 6, the latter two having a characteristic gradient length scale $w$. All of these profiles produce two chiral channels near the $K$ point and two chiral channels near the $K^{\prime}$ point. This is a result of the bulk topology, as will be described in the next section. The primary difference between the three profiles is the wavefunction localization; for a fixed $V$, the step potential leads to the most localized states. For simplicity, only the step function will be considered for the rest of the paper, though it is important to note that the localization can also be experimentally tuned by modulating the electric field. In the limit of small $w$ for a smoothly varying $V(y)$, one would then observe states with a localization approaching those in a step-function profile.

When a magnetic field is absent, bilayer graphene with the potential interface recovers a screw symmetry along the potential interface on the $x$-axis (see figure 1). Explicitly, we may write the symmetry transformation as $H(x,y,z,k_{x})\to H(x+a/2,-y,-z,k_{x})$, where $z\to-z$ describes a layer interchange. The screw symmetry ensures equal probability distribution between the two layers. The screw symmetry is represented by a unitary matrix $S_{2}$ which squares to a unit translation in $x$, $S_{2}^{2}=e^{i{\bf k}\cdot{\bf a}_{1}}=e^{ik_{x}a}$. Therefore the screw eigenvalues are $\pm e^{ik_{x}a/2}$. Figure 8 shows the effect of the potential on the argument of the screw eigenvalues for the zero-energy modes. In the $V=0$ limit, they begin at $\pm e^{i\pi/3}=\pm e^{iK\cdot{\bf a}_{1}/2}=\pm e^{iK_{x}a/2}$. Increasing or decreasing the potential shifts the Fermi momenta of the zero-modes to $K_{x}\pm\delta k_{x}$, where $\delta k_{x}$ takes the analytic form in (16) in the $\mathrm{k}\cdot\mathrm{p}$ approximation. The screw eigenvalues of the two zero modes become $\pm e^{i(\pi/3\pm\delta k_{x}a/2)}=\pm e^{i(K_{x}\pm\delta k_{x})a/2}$. The evaluation of the argument of the screw eigenvalue $\langle\psi_{0}|S_{2}|\psi_{0}\rangle$ therefore provides a method to calculate the Fermi momentum shift $\delta k_{x}$ of the zero-mode $\psi_{0}$ in the discrete lattice tight-binding model. These values are shown as discrete points in Figure 8, and are well-approximated by the $\mathrm{k}\cdot\mathrm{p}$ approximation represented by the continuous curve.

V Bulk topology and symmetry considerations

The existence of chiral modes is a consequence of the change of bulk topology across the electric interface. By the bulk-boundary correspondence Volovik (2003); Nakahara (2003); Chiu et al. (2016), the number of gapless (valley) Dirac edge modes is identical to the change in (valley) Chern numberThouless et al. (1982); Zhang et al. (2013) across the edge. The Chern number $\mathrm{Ch}_{1}$ is a topological property of the bulk, so by studying the topology of bilayer graphene with a perpendicular electric field, we can determine the behavior at the interface of these two configurations. $\mathrm{Ch}_{1}$ is typically calculated by integrating the Berry curvature $\mathbf{\Omega}$ over the Brillouin zoneTse et al. (2011)

[TABLE]

where the sum is taken over the two occupied bands $n$ below the Fermi level at 0 energy and

[TABLE]

For bilayer graphene, the integral of the Berry curvature over the entire Brillouin zone vanishes as a result of time-reversal symmetry which requires the Berry curvature to be an odd function.Xiao et al. (2010) However, the Berry curvature is localized near the $K/K^{\prime}$ points, so one may consider the integral of the Berry curvature near these points. Using k$\cdot$p perturbation, we found the Berry curvature to be an odd function of $V$ and $\nu$. Integrating the local Berry curvature over all space gives the valley Chern number, though this is difficult to calculate using this formulation.

In the non-interacting limit, there is an equivalent formulation of the Chern number given in terms of the Green’s function.Martin et al. (2008); Wang and Zhang (2012); Tse et al. (2011); Volovik (2003) Here, the Green’s function takes on the form $G=(\tilde{E}-H(\mathbf{k}))^{-1}$

[TABLE]

where $\tilde{E}$ represents a complex energy and the energy integral is taken over a contour around the valence bands. The indices run from 0 to 2, where $\partial_{0}$ is an energy derivative and $\partial_{1}$ and $\partial_{2}$ are $k_{x}$ and $k_{y}$ derivatives, respectively. This expression is simpler to evaluate since it only involves derivatives of $\tilde{E}-H(\mathbf{k})$ rather than of the gauge-dependent normalized eigenstates as in (23). Since $\partial_{0}G^{-1}=1$, we have that

[TABLE]

The integrand only depends on the magnitude of the momentum, so the momentum integral can be evaluated in polar coordinates. The energy integral is found by evaluating the two residues at negative real energy to obtain

[TABLE]

At the interface where $V$ changes sign, the difference in the Chern number is $-2/+2$ at the $K/K^{\prime}$ point, in agreement with our previous observations of two left/right moving modes at the $K/K^{\prime}$ point as obtained from the band structure.

V.1 Symmetry protected topology

The spinless bilayer graphene model (1) under a uniform electric field preserves time-reversal symmetry and charge $U(1)$ conservation, and thus belongs to class AI according to the tenfold classification Altland and Zirnbauer (1997) of band theory. In addition to these local symmetries, the model is also symmetric under the non-centrosymmetric space group (wallpaper group) P3m1, which is broken from the centrosymmetric P$\bar{3}$m1 by the layer-asymmetric electric potential. In particular, we focus on the point group symmetry $C_{3v}=\{1,r,r^{2},\mu_{x},r\mu_{x},r^{2}\mu_{x}\}$, which is generated by a threefold rotations $r$ about the $z$-axis and a vertical mirror plane $\mu_{x}$ perpendicular to the $x$-axis. The twofold screw rotation described previously (in section IV and figure 8) only applies along the domain wall potential interface and is absent in the bulk where the electric field is uniform. Together with the time-reversal symmetry group $\mathbb{Z}_{2}^{T}$, which is generated by the anti-unitary time-reversal operator $T$ that squares to $T^{2}=+1$, they form the magnetic point group

[TABLE]

In this subsection, we discuss the symmetry-protected topology.

Mirror symmetry requires the Berry curvature to obey the asymmetry $\Omega_{n}({\bf k})=-\Omega_{n}(\mu_{x}{\bf k})$, where the minus sign is associated to the orientation reversing nature of reflection. Similar asymmetry relations hold for the other two mirror planes, $\Omega_{n}({\bf k})=-\Omega_{n}(r\mu_{x}{\bf k})=-\Omega_{n}(r^{2}\mu_{x}{\bf k})$. Therefore, as curvatures cancel between mirror opposite momentum points, the Chern number (22), when integrated over the entire Brillouin zone, must be identically zero. In particular, along the three mirror-symmetric lines $L_{1}=\Gamma M_{1}\Gamma$, $L_{2}=\Gamma M_{2}\Gamma$ and $L_{3}=\Gamma M_{3}\Gamma$ where $\mu_{x}{\bf k}={\bf k}$, $r\mu_{x}{\bf k}={\bf k}$ and $r^{2}\mu_{x}{\bf k}={\bf k}$ respectively (see the dashed lines in figure 10), the Berry curvature must vanish, $\Omega_{n}|_{L_{j}}=0$. Here, $M_{1},M_{2},M_{3}$ are the three inequivalent time-reversal symmetric momenta apart from the Brillouin zone origin $\Gamma$. The three mirror-symmetric lines enclose a triangular region $\mathrm{BZ}^{1/2}_{+}$ (represented by the blue shaded region in figure 10) that contains a single $K$ point and traces out half of the Brillouin zone. Mirror (or time-reversal) maps $\mathrm{BZ}^{1/2}_{+}$ to another triangular region $\mathrm{BZ}^{1/2}_{-}$ that contains $K^{\prime}$. Together, they generate the entire Brillouin zone $\mathrm{BZ}=\mathrm{BZ}^{1/2}_{+}\cup\mathrm{BZ}^{1/2}_{-}$. These two mirror-conjugated half-regions are almost mutually disjoint as they only intersects along the mirror-symmetric lines, $\mathrm{BZ}^{1/2}_{+}\cap\mathrm{BZ}^{1/2}_{-}=L_{1}\cup L_{2}\cup L_{3}$.

The valley Chern number is defined to be the integral

[TABLE]

over the half-Brillouin zone containing only a single valley. We now show that the point group symmetry $C_{3v}$ guarantees that the valley Chern number can only take discrete integral values. Consequently, $v\mathrm{Ch}_{1}$ is stable against any perturbation that preserves the $C_{3v}$ symmetries and excitation energy gap, and defines a symmetry-protected topological invariant. The Berry curvature is $d$-exact on the half-Brillouin zone, and is identical to the differential $\mathbf{\Omega}_{n}=\nabla_{\bf k}\times\boldsymbol{\alpha}_{n}$, where the Berry connection is

[TABLE]

From the Stokes’ theorem, the valley Chern number is

[TABLE]

where the boundary of the half-Brillouin zone consists of the three mirror-symmetric lines, $\partial\mathrm{BZ}^{1/2}_{+}=L_{1}\cup L_{2}\cup L_{3}$. Without symmetries, the holonomy (also known as polarization)

[TABLE]

can take any real value. Time-reversal symmetry requires them to take integral values Fu and Kane (2006); Qi et al. (2008); Chiu et al. (2016). Each one of the mirror-symmetric lines $L_{j}$ is closed under time-reversal symmetry in the sense that if ${\bf k}$ belongs in $L_{j}$, so is its time-reversal conjugate $T{\bf k}=-{\bf k}$. Moreover, each line is a closed loop and is topologically equivalent to the 1D Brillouin zone. The band Hamiltonian $H({\bf k})|_{L_{j}}$ restricted on each of these lines is hence identical to a 1D time-reversal symmetric band insulator, which is known to have integral electric polarization. Combining the polarizations, the valley Chern number (30) therefore has integral value. Lastly, since the Berry curvature is gauge invariant, so is the valley Chern number (28) and the sum of polarizations along $L_{1}$, $L_{2}$ and $L_{3}$ in (30).

We notice that the symmetry-protected topology actually only relies on the combination of mirror and time-reversal symmetry, rather than the individuals. These combinations form a symmetry subgroup

[TABLE]

inside the full magnetic point group $3m1^{\prime}$ in (27). The valley Chern number must still take integral values based on this subset of magnetic point group symmetries. This is because the mirror-time-reversal combinations $T\mu_{x},Tr\mu_{x},Tr^{2}\mu_{x}$ take the same role as the local time-reversal $T$ along the mirror-symmetric lines $L_{1},L_{2},L_{3}$, and they still enforce the integrality of the polarizations $P_{j}$. In other words, the valley Chern number continues to provide a topological characterization of the system even when time-reversal symmetry is broken as long as the combined symmetries in (32) are preserved. This applies, for instance, in the presence of a magnetic field on the bilayer graphene. This is because magnetic field is a pseudo-vector and is flipped under any improper rotation such as mirror. At the same time, it is also odd under time-reversal. Consequently, a magnetic field ${\bf B}=B_{z}\hat{\bf z}$ along the perpendicular direction is invariant under the magnetic point group $3m^{\prime}$ in (32). Although some of the magnetic-mirror symmetries may be broken by the electric domain wall potential interface or a particular gauge choice ${\bf B}=\nabla\times{\bf A}$, we speculate that the robustness of the topological chiral channel along an electric domain wall in the presence of a magnetic field may be a consequence of the lingering non-trivial bulk topology protected by a magnetic symmetry. The numerical calculations presented in the following section will provide results consistent with this claim.

VI Strong Electric and Magnetic Field

Tight-binding approximations were used to study bilayer graphene in the presence of strong electric and magnetic fields perpendicular to the bilayer in the lattice limit. The breaking of translation symmetry in $y$ with the introduction of an electric interface requires us to extend our unit cell to include the entire vertical dimension. To preserve translational symmetry in $x$, the Coulomb gauge $\mathbf{A}=-By\hat{\mathbf{x}}$ is chosen with $y=0$ set at the interface.

Introducing a magnetic field requires the Peierls substitution in the hopping parameter

[TABLE]

where the integral is evaluated along the path between nearest-neighbor sites, and $\phi_{0}=h/e$ is the magnetic flux quantum. It will be useful to know when the magnetic energy scale is comparable to the electric energy scale set by $V$. The magnetic energy scale is set by the spacing between Landau levels. From equation (7), $E_{2}-E_{1}=2(\hbar v_{F})^{2}eB\sqrt{2}/\hbar u\approx(0.005$ eV/T$)B$, so for $V\approx 0.1$ eV, a comparable magnetic energy scale requires a field of about $B\approx 20$ T. For many of the plots shown in this section, we use these values for $V$ and $B$.

We numerically study this comparable limit in the tight binding model by considering lattice commensurate magnetic filling fractions, where the magnetic field takes on the form

[TABLE]

where the rational number $p/q$ is the amount of magnetic flux (in unit of the flux quantum $\phi_{0}$) through a single hexagon plaquette. To study the topological chiral channels along the electric interface, we choose an open geometry along the vertical $y$ direction. Although a closed cylindrical geometry could avoid the irrelevant edge modes, the vector potential would become discontinuous because the geometry would enclose magnetic monopoles. The Peierls substitution (33) would in general be discontinuous as well unless the system circumference $L_{y}$ is also in some commensurate length. The discontinuity would correspond to an unphysical edge where the magnetic field diverges and additional edge modes arise. On the other hand, the Coulomb gauge vector potential is continuous in an open geometry. In a large system where the boundary edges at $y=L_{y}/2$ and $-L_{y}/2$ are far separated from the electric interface at $y=0$, the edge modes do not mix with interface states due to their exponentially localized wavefunction thanks to the bulk energy gap. For all of our calculations, we use at least $L_{y}=2000$ atoms in an open edge geometry. Although the edge modes are still visible in the electronic band structure, they are not of interest in our study. These irrelevant modes can be discarded by focusing only on quantum states that localized along the electric interface. Numerical analysis will be given for the relevant chiral modes near the $K$ point, and an explanation will be given on how to relate these results to the chiral modes near the $K^{\prime}$ point.

The band structure for bilayer graphene with an electric domain wall and no magnetic field indicates two right-moving modes near the $K$ point and two left-moving modes near the $K^{\prime}$ point, as explained in section IV. Plotting the wavefunctions corresponding to the chiral modes confirms that these states are exponentially localized at the domain wall as given by equation 20. These wavefunctions have no dependence on the Fermi energy. The localization $\sigma_{y}$ in the numerical lattice tight binding model agrees with our continuum calculations, as shown in Figure 7.

When both $B$ and $V$ are non-negligible, both flat bands (Landau levels) and interface modes can be observed in the band structure (Figure 11). Although classically one would expect electrons in a magnetic field to form cyclotron orbits leading to Landau levels at zero energy, the fixed chirality of the interface states prevents electrons from completing cyclotron orbits Zarenia et al. (2011). The chiral modes are still localized near the interface, but their mean position $\expectationvalue{y}$ is not located exactly at the interface due to the formation of a secondary peak in the wavefunction (Figure 12). Electronic transport in this configuration would occur parallel to, but shifted away from, the electric domain wall. A nonzero Fermi energy also contributes to an overall shifting of the wavefunction peak (Figure 13). This shifting only occurs when the magnetic field is nonzero.

Classically, one might expect the shifting of the edge modes from the electric domain wall to be related to how moving charges are bent in a magnetic field. By the Lorentz force, electrons on two different channels near $K$ moving in the same direction (or equivalently, electrons in the same valley) should be bent in the same direction, but we have found the co-propagating modes shift in opposite directions as the magnetic field is increased. The sign of the Fermi velocity is the same for the two modes at a fixed valley, but the force exerted on the electrons is opposite. Therefore, the Lorentz force does not provide an explanation of this phenomenon. We will later show that the shifting can be understood by the overlap of the electric and magnetic states at a fixed momentum.

At a fixed Fermi energy, the deviation from the electric interface is directly proportional to the magnetic field strength with direction determined by the sign of B (Figure 14). Specifically, $\expectationvalue{y}(B)=-\expectationvalue{y}(-B)$. This shifting of $\expectationvalue{y}$ away from the interface is a result of the secondary peak in the wavefunction previously mentioned. Here, the electric field keeps the states localized to the interface while the magnetic field tends to shift the states away from the interface. Increasing the Fermi energy exaggerates the shifting for $B>0$ (Figure 15). As one might expect, a small change in the Fermi energy does not have a significant effect on the shifting for larger $V$ because the Fermi level remains close to the center of the bulk band gap determined by $V$ (Figure 16).

The shifting can be better understood by partitioning the wavefunction into components residing on the top layer and those residing on the bottom layer. We use $\psi_{-}^{t}$ to denote the component of $\psi(K-\delta k_{x})$ that lives on the top layer and $\psi_{-}^{b}$ for the bottom component of $\psi(K-\delta k_{x})$. Similar notation is used for the top and bottom components of $\psi(K+\delta k_{x})$. Note that while $\psi_{+}$ and $\psi_{-}$ shift in opposite directions, as $V$ becomes much larger than $B$, both wavefunctions approach the interface (Figure 17). Figure 18 describes the shift and localization of the wavefunction as a function of $B$ for fixed $E_{f}$ and $V$. The peak in the figure on the right showing localization length $\sigma_{y}$ as a function of magnetic field is a result of two competing effects. For small magnetic fields, the magnetic contributions to the wavefunction are negligible, so $\sigma_{y}$ approaches that of the electric state. For very large magnetic fields, the magnetic states sit at the interface with a greater localization, consistent with equation 10. The peak occurs when the magnetic field strength is intermediate between these two scenarios. Then both the electric and magnetic states have non-negligible contributions, and the net effect is that the overall wavefunction broadens. The localization is identical for the four different chiral modes.

As the Fermi energy is increased, the two chiral modes in a given valley shift in the same direction (Figure 19). The dependence of $\expectationvalue{y}$ on the Fermi energy is to be expected because $\expectationvalue{y}$ depends linearly on the Fermi momentum by equation 10. This differs from the $B=0$ case, where $\expectationvalue{y}=0$ regardless of the Fermi energy.

As discussed before, bilayer graphene with $B=0$, $V\neq 0$ obeys a screw symmetry along the axis of the interface. Since the vector potential is $y$-dependent, the symmetry is broken and a layer preference is allowed to form (Figure 20). The layer preference depends on the valley index but not the sign of $\delta k_{x}$, consistent with the $V=0$ case. Furthermore, the layer preference has no dependence on the Fermi energy, since changing the Fermi energy only changes the Fermi momentum, and the layer preference does not depend on $\delta k_{x}$ as given in (12).

In order to understand why a magnetic field causes a peculiar shifting of the wavefunctions, we define $\psi_{BV}$ to be the linear combination of the electric and magnetic states

[TABLE]

where $k_{x}=k_{x}(B,V)$ takes on the value of one of our four Fermi momenta and the coefficients $\alpha$, $\beta_{1}$, and $\beta_{2}$ are given by inverting the following matrix equation.

[TABLE]

While the magnetic states are orthogonal to each other, the inner product between a magnetic state and the electric state generally does not vanish. Qualitatively, one might hope that $\psi_{BV}$ describes the actual wavefunction well. Our numerical results confirm that this is a good approximation for experimentally realizable fields (Figure 22). The magnitudes of $\alpha$ and $\beta$ give an indication of in what regimes $\psi_{V}$ dominates.

Employing our approximation leads to a more intuitive understanding of the microscopic properties of the composite system. From our study of the magnetic states, we expect $\psi_{B}$ to shift in $y$ in proportion to $-\nu\delta k_{x}$. However, $\psi_{V}$ is localized at $y=0$ with momentum fixed by (16). For this momentum, $\psi_{B}$ is localized at nonzero $y$. Therefore, the combination of these two separate wavefunctions with fixed momentum causes two peaks to form, the larger primary one directly on the interface is the contribution from $\psi_{V}$ and the smaller secondary hump shifted slightly from the interface is the contribution from $\psi_{B}$. It is important to stress that the hump approaches the interface with increasing $B$ consistent with equation 10, but this causes a greater overlap with the electric and magnetic states, increasing the magnitude of the hump, and leading to an overall probability distribution shifted away from the interface for experimentally realizable fields. In fact, for a small electric field $V\approx.1$ eV, the states will shift away from the interface with increasing $B$ for $B\lessapprox 50$ T. Figure 12 displays this behavior. This explains the counter-intuitive notion that increasing the magnetic field shifts the magnetic wavefunction toward the interface, yet it shifts the overall wavefunction away from the interface for values of $V$ and $B$ used in experiment. In our approximation, this is equivalent to the more intuitive idea that strengthening $B$ increases the weight of our magnetic states in $\psi_{BV}$, as given by $|\beta|^{2}$.

VII Discussion and Conclusion

It is well-known that the band gap of bilayer graphene can be tuned by changing the interlayer potential difference, such as by applying a perpendicular electric field. A perpendicular magnetic field provides an additional means of controlling the system. The magnetic field causes peculiar yet predictable alterations to the electronic wavefunction at electric domain walls. The lateral shifting and asymmetric layer distribution could partially explain the suppressed backscattering in the presence of a magnetic field as detected in experiments Li et al. (2016, 2018). With no magnetic field, the counter-propagating modes at $K\pm\delta k$ and $K^{\prime}\mp\delta k$ are located at the same position in real space. An introduction of a magnetic field separates these wavefunctions, decreasing their spatial overlap. For example, at $V=0.1$ eV and $E_{f}=0$, $|\innerproduct{\psi^{*}(K-\delta k)}{\psi(K^{\prime}+\delta k)}|^{2}$ equals 1 when $B=0$ but equals 0.73 when $B=12$ T. While the counter-propagating modes at $K\pm\delta k$ and $K^{\prime}\pm\delta k$ still lie at the same y-value, they are localized to opposite layers with a nonzero magnetic field. This layer distribution is a new observation that could be detected using scanning tunneling microscopy.

The applications for valleytronic devices and transport measurements motivate a theoretical study for bilayer graphene with a spatially varying potential difference and a magnetic field. Attempting to diagonalize this Hamiltonian in the continuum limit results in a system of differential equations that admits no elementary solutions, its direct treatment being restricted to numerical methods. We present results showing that the actual solution can be approximated as a linear combination of the electric and magnetic states, all of which possess analytic forms. This combination of numerical and analytical results allows for high accuracy as well as a qualitative understanding of the microscopic properties. Unlike results obtained perturbatively, these results hold even in the regime that both $V$ and $B$ are large. Treating $\alpha$, $\beta_{1}$, and $\beta_{2}$ as pure numbers that can be found numerically by equation 36, one may find equations for the average position, localization, and layer preference of the combined electric and magnetic state by using the results for the purely electric and purely magnetic states.

Acknowledgements.

We thank Prof. Jun Zhu for insightful and inspiring discussion that initiated this work. We thank Matt Walker and Juan Velasco for helping with calculations and figures. We thank the Mead Endowment for supporting the undergraduate research program. JCYT is supported by the National Science Foundation under Grant No. DMR-1653535.