text:
\section{Introduction}\label{sec:intro}
In recent years, there has been considerable interest in the electronic properties of bilayer graphene~\cite{McCann07,NetoGuineaPeresNovoselovGeim09,McCann13,ROZHKOV20161} including a quadratic dispersion with Berry phase of $2\pi$ and an unusual quantum Hall effect~\cite{Kovoselov06,FalKo08}. The inequivalence of the $K$ and $K'$ 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~\cite{RycerzTworzydloBeenakker07,Behnia12,AngYangZhangZhongshuiAng17}. New techniques have also made it possible to produce large flakes of Bernal (AB-stacked) graphene on the order of $10^{-4}$ m,~\cite{Ji17} 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~\cite{Min07} and experimentally by applying an electric field perpendicular to the layers~\cite{Zhang09}. In pure samples where inter-valley mixing is suppressed~\cite{MorozovNovoselovKatsnelsonSchedinPonomarenkoJiangGeim06,McCannKechedzhiFalkoSuzuuraAndoAltshuler06,MorpurgoGuinea06}, a non-local valley symmetry emerges and the insulating electronic system carries a non-trivial valley Chern number~\cite{ZhangMacDonaldMele13}. This symmetry-protected topological nature~\cite{VaeziLiangNgaiYangKim13} of the ground state is revealed by a one-dimensional highway of electrons along an electric domain wall~\cite{MartinBlanterMorpurgo08,QiaoJungNiuMacDonald11,JungZhangQiaoMacDonald11,ZareniaPereiraFariasPeeters11}, 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~\footnote{Except along commensurate ones, such as the armchair edge, where the valley momentum difference is projected out, ${\bf a}_{\parallel}\cdot(K-K')\in2\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~\cite{Ju15,Li16,Lee17,Li18}. 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.~\onlinecite{Lee17} using bilayer graphene encapsulated between atomically clean hexagonal boron nitride single crystals, which has been implemented in similar systems to reduce intervalley scatterings~\cite{DeanHone10,KimLee16}. 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.~\onlinecite{Ju15} reported a domain wall mean free path of $l_0\sim400$ nm, and Ref.~\onlinecite{Li16} inferred a mean free path of $l_0\sim200$ 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~\cite{Datta95}. 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.~\onlinecite{Li16} at around $B=8$ T for $L=0.4$ $\mu$m, and similar ballistic transport has been reported in Ref.~\onlinecite{Li18} under similar field strength and length scale. The realization of differential conductance approaching the ballistic limit without a magnetic field in Ref.~\onlinecite{Lee17} 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~\cite{ZareniaPereiraFariasPeeters11,WangRenDengYangJungQiao17}. In particular, we focus on the asymmetric lateral and layer distribution~\cite{Li16} 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~\cite{Li18}, and introduces new complexities to the critical transport behaviors of bilayer graphene domain wall quantum point contacts~\cite{WiederZhangKane15}. 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 \ref{sec:tightbinding}, 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 \ref{sec:landau} 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 \ref{sec:chiral}. 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 \ref{sec:top} including an argument for why the valley Chern number is a topological invariant of our system. In Section \ref{sec:strong}, 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 \ref{sec:con}.

\begin{figure}[b!] \label{fig:bilayerdiagram}
\centering
\includegraphics[width=0.45\textwidth]{bilayerlattice.pdf}
\caption{AB-stacked bilayer graphene with electric domain wall at $y = 0$.} \label{blg_diagram}
\end{figure}

\section{Tight Binding Hamiltonian} \label{sec:tightbinding}

We review a tight-binding model Hamiltonian~\cite{McCann13} 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 $\vb{a_{1}} = a(1,0)$ and $\vb{a_{2}} = a(-1/2, \sqrt{3}/2)$ with lattice constant $a = 2.46$ \AA\cite{Saito98}, 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 interaction\cite{Yao07}, 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$ eV\cite{McCann13}, interlayer hopping parameter $u = 0.381$ eV between the $A1$ and $B2$ sites that lie directly on top of each other\cite{McCann13}, and electric potential energy difference $V$ between the layers, the tight-binding Hamiltonian can be written as
\begin{align}
    H = &- t\sum_{\ev{{\bf r},{\bf r}'}, s}(a_{{\bf r},s}^{\dagger}b_{{\bf r}',s} + \text{h.c.}) - u\sum_{{\bf r}}({a}_{{\bf r},1}^{\dagger}b_{{\bf r},2} + \text{h.c.}) \nonumber\\
    &-\frac{V}{2}\sum_{{\bf r},s}(-1)^s(a_{{\bf r},s}^{\dagger}a_{{\bf r},s} + b_{{\bf r},s}^{\dagger}b_{{\bf r},s}) \label{TBH}\\
		=&\int\frac{d^2{\bf k}}{(2\pi)^2}{\bf c}_{\bf k}^\dagger H({\bf k}){\bf c}_{\bf k}\nonumber
\end{align}
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
\begin{equation} \label{eq:bloch}
H(\vec{k})=\begin{pmatrix}
V/2 & tf(\vec{k}) & 0 & u\\
tf^{*}(\vec{k}) & V/2 & 0 & 0\\
0 & 0 & -V/2 & tf(\vec{k})\\
u & 0 & tf^{*}(\vec{k}) & -V/2
\end{pmatrix}
\end{equation}
where $f(\vb{k}) = 1 + e^{-i\vb{k}\vdot\vb{a_{2}}} + e^{-i\vb{k}\vdot(\vb{a_{1}} + \vb{a_{2}})}$ are the interlayer hopping terms in momentum space describing the three nearest-neighbor points. Diagonalizing (\ref{eq:bloch}) gives a quadratic dispersion with zero energy gap when $V=0$ at the inequivalent $K$ and $K'$ points, on the edge of the Brillouin zone (Figure \ref{fig:mkgamma}). An energy gap $\Delta$ is introduced with a non-vanishing potential
\begin{equation}
    \Delta = \frac{V}{\sqrt{1 + V^{2}/u^{2}}} \approx \begin{cases}V, & V \ll u \\ u, & V \gg u\end{cases}.
\end{equation}

\begin{figure}
    \centering
    \raisebox{-.5\height}{\includegraphics[width=.12\textwidth]{BZ.png}}
    \raisebox{-.5\height}{\includegraphics[width=.33\textwidth]{gmk.png}}
    \caption{Band structure (right) of pristine bilayer graphene along a parameterized path through the $M$, $K$, and $\Gamma$ points in the Brillouin zone (left). The inset shows a gap at the $K$ point for $V = 0.1$ eV.}
    \label{fig:mkgamma}
\end{figure}

If we are only concerned with low energy perturbations, we may expand $H$ about the $K$ and $K'$ points
\begin{align}
H_{\mathrm{k}\cdot\mathrm{p}} &= [H_0({\bf k})+\nabla_{\bf k}H({\bf k})\vdot \delta {\bf k}]\big|_{{\bf k}=K,K'} \nonumber\\
&= \begin{pmatrix}
V/2 & \hbar v_F\pi & 0 & u\\
\hbar v_F\pi^{\dagger} & V/2 & 0 & 0\\
0 & 0 & -V/2 & \hbar v_F\pi \\
u & 0 & \hbar v_F\pi^{\dagger} & -V/2
\end{pmatrix}\label{eq:kp}\\
&=\hbar v_F(\nu\delta k_x\sigma_x-\delta k_y\sigma_y)\otimes\openone_2\nonumber\\
&\;\;\;+\frac{V}{2}\openone_2\otimes\tau_z+\frac{u}{2}(\sigma_x\otimes\tau_x-\sigma_y\otimes\tau_y)\nonumber
\end{align}
 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'$. Here, $\delta k_x$ and $\delta k_y$ are small momenta deviations away from the Dirac points, and the $2\times2$ Pauli matrices $\sigma$ ($\tau$) act on the $AB$-sublattice (layer) degree of freedom.

\section{Landau Levels} \label{sec:landau}
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 $\vb{B} = B\vu{z}$. The momentum in the tight-binding Hamiltonian is replaced by the canonical momentum
\begin{equation}
    {\bf k}\longrightarrow {\bf k}-\frac{e{\bf A}}{\hbar}.
\end{equation}
The $\mathrm{k}\cdot\mathrm{p}$ Hamiltonian in \eqref{eq:kp} has the symmetry $[H_{\mathrm{k}\cdot\mathrm{p}}(K),Z]=[H_{\mathrm{k}\cdot\mathrm{p}}(K'),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=\pm1$. Focusing on one of the two eigenspaces, the energy spectrum near the Fermi level may be found by the effective two-band Hamiltonian\cite{McCann06,Pereira07}
\begin{equation}
    H_{\mathrm{eff}} = -\frac{v_{F}^{2}}{u}\begin{pmatrix}0 & \pi^{2} \\ \qty(\pi^{\dagger})^{2} & 0\end{pmatrix}
\end{equation}
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
\begin{equation} \label{eq:landau}
    E_{n} = \operatorname{sgn}(n)\frac{2\hbar^{2}v_{F}^{2}}{ul_{B}^{2}}\sqrt{\abs{n}(\abs{n}-1)}
\end{equation}
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 $\vb{A} = -By\vu{x}$, our tight-binding Hamiltonian about the Dirac points is then given by (\ref{eq:kp}) 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_xx}$) are
\begin{align}
    \psi_{B,1}(y) &= \frac{1}{\sqrt{\mathcal{N}}}e^{-(y+l_B^2 \delta k_x)^2/2l_{B}^{2}}
    \begin{pmatrix}0 \\ \sqrt{l_B^2/2+(\hbar v_F)^2/u^2} \\ 0 \\ 0 \end{pmatrix} \nonumber\\
    \psi_{B,2}(y) &= \frac{1}{\sqrt{\mathcal{N}}}e^{-(y+l_B^2 \delta k_x)^2/2l_{B}^{2}}
    \begin{pmatrix}0 \\ y+l_B^2 \delta k_x \\ 0 \\ -(\hbar v_F)/u \end{pmatrix}
\end{align}
where the normalization constant is
\begin{equation}
\begin{aligned}

		\mathcal{N} &= \frac{\sqrt{\pi} l_B^3}{2}\left[1+2\left(\frac{\hbar v_F}{u l_B}\right)^2\right]
\end{aligned}
\end{equation}
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 $\ev{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.
\begin{equation}
   \ev{y} = -l_{B}^{2} \delta k_{x}+\frac{l_B^2 u\Re(\beta_{1}^{*}\beta_{2})}{\sqrt{(\hbar v_F)^2+l_B^2 u^2/2}}
   \label{eq:yB}
\end{equation}
This result for $K'$ is similar in that $\ev{y}\propto -l_B^2 \delta k_x$.

The localization length of $\psi_B$, $\sigma_y = \sqrt{\ev{ y^2} -\ev{y}^2}$, has no dependence on $\delta k_x$.
\begin{equation}
    \sigma_y^2 = \frac{l_B^2 (\hbar v_F)^2+l_B^4 u^2[1+\Re(c_1^{*}c_2)-2\Im(c_1^{*}c_2)^2]}{2 (\hbar v_F)^2+l_B^2 u^2}
\end{equation}

\begin{figure}
    \centering
    \includegraphics[width=.45\textwidth]{sigma_B.png}
    \caption{Localization of magnetic states. For large fields, the localization approaches $l_B/\sqrt{2}.$}
    \label{fig:stddev_purelyB}
\end{figure}

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}$
\begin{equation}
    P_{\text{bottom}} = \frac{1+|\beta_{2}|^2-|\beta_{1}|^2}{2+l_B^2 u^2/(\hbar v_F)^2}.
    \label{eq:pbot}
\end{equation}
$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}}$
\begin{equation}
    0\leq P_{\text{bottom}} \leq \frac{2}{2+l_B^2 u^2/(\hbar v_F)^{2}}.
\end{equation}
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~\cite{Li16}.
\begin{equation}
P_{\text{bottom}} = \frac{e \hbar v_F^2}{u^2} \qty(1+|\beta_{2}|^{2}-|\beta_{1}|^{2}) B.
\end{equation}
At the $K'$ point, the opposite layer preference was found.

\begin{figure}
    \centering
    \includegraphics[width=.5\textwidth]{Ptop_B.png}
    \caption{Probability for an electron to occupy the top layer. Any magnetic state that is a linear combination of these two states will have a similar curve that lies in between these two curves. For realistic fields, the layer preference is approximately linear in B.}
    \label{fig:prob_purelyB}
\end{figure}

\section{Chiral Modes} \label{sec:chiral}
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 (\ref{eq:kp}) 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 \ref{blg_diagram}. 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
\begin{equation}
    \kappa = \frac{\abs{V}}{2\sqrt{2}\hbar v_f}\sqrt{-1+\sqrt{1+8\qty(\frac{u}{V})^2}}
\end{equation}
The value of $k$ is to be interpreted as the momentum where our chiral channels are located in reciprocal space.
\begin{equation}
    \label{eqn:dkx}
    k = \delta k_{x}=\frac{\abs{V}}{4\hbar v_F}\sqrt{1+\sqrt{1+8\qty(\frac{u}{V})^{2}}}.
\end{equation}
The solution takes the following form:
\begin{equation}
\psi_{V} (y) =
\begin{cases}
\psi_{V}^{-}(y), & y<0\\
\psi_{V}^{+}(y), & y\geq 0\\
\end{cases}
\end{equation}
\begin{equation}
\psi_{V}^{-}(y)=
\text{\footnotesize $
c_{V} \begin{pmatrix}
(V/2\hbar v_f) e^{i k y}+i (V/2\hbar v_f) e^{-i k y} \\ ((1-i)k-\kappa)e^{i k y}+((-1+i)k-i\kappa)e^{-i k y} \\ ((1-i)k-i \kappa)e^{i k y}+((-1+i)k-\kappa)e^{-i k y} \\ i (V/2\hbar v_f)e^{i k y}+(V/2\hbar v_f)e^{-i k y}
\end{pmatrix} e^{\kappa y}$}
\end{equation}
\begin{equation}
\psi_{V}^{+}(y)=
\text{\footnotesize $
c_{V} \begin{pmatrix}
(V/2\hbar v_f)e^{i k y}+i (V/2\hbar v_f) e^{-i k y} \\ ((-1+i)k-\kappa)e^{i k y}+((1-i)k-i \kappa) e^{-i k y}  \\ ((-1+i)k-i \kappa)e^{i k y}+((1-i)k-\kappa) e^{-i k y} \\ i (V/2\hbar v_f)e^{i k y}+(V/2\hbar v_f) e^{-i k y}
\end{pmatrix} e^{-\kappa y}$}
\end{equation}
where the normalization constant is
\begin{equation}
c_{V}=\frac{\sqrt{\kappa}}{4 k}.
\label{eq:psiv}
\end{equation}
This describes a total of four zero-energy modes in the Brillouin zone, at $K\pm \delta k_x$ and $K' \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'$ 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$ \AA $\ $in the $V \gg 1$ eV limit.
\begin{equation}
\sigma_y = \sqrt{\frac{1}{2\kappa^2}+\frac{\kappa^2}{2(k^2+\kappa^2)^2}}
\end{equation}

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.

\begin{figure
    \centering
    \includegraphics[width=.45\textwidth]{potentials.png}
    \caption{Three potential profiles used for calculations. The plots show the potential on the top layer; the bottom layer is given by the negative of this value. The linear and sinusoidal profiles are described by a width $w$, the characteristic gradient length scale.} \label{fig:potprofs}
    \includegraphics[trim={3cm 0 0 0},clip,width=.5\textwidth]{StdDev_w.png}
    \caption{Localization of an interface state for different potential profiles as a function of $w$. Parameters used in the calculation are $E_{f}=0$ and $V=0.1$ eV. The dashed line indicates the step potential considered above.}
\end{figure}

We considered three different potential profiles: step-function, linear, and sinusoidal, as seen in Figure \ref{fig:potprofs}, 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'$ 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.

\begin{figure}
    \centering
    \includegraphics[width=.45\textwidth]{StdDev_V.png}
    \caption{Localization of an interface state as a function of V. Plane wave solutions are recovered when the electric interface vanishes ($V=0$).}
    \label{fig:std_vs_V}
\end{figure}

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~\ref{blg_diagram}). 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_xa}$. Therefore the screw eigenvalues are $\pm e^{ik_xa/2}$. Figure~\ref{fig:evals_screw} 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_xa/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 \eqref{eqn:dkx} 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_xa/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~\ref{fig:evals_screw}, and are well-approximated by the $\mathrm{k}\cdot\mathrm{p}$ approximation represented by the continuous curve.

\begin{figure}
    \centering
    \includegraphics[width=.45\textwidth]{screwevals.png}
    \caption{Argument of eigenvalues of screw operator at $K-\delta k_x$.}
    \label{fig:evals_screw}
\end{figure}

\section{Bulk topology and symmetry considerations} \label{sec:top}
The existence of chiral modes is a consequence of the change of bulk topology across the electric interface. By the bulk-boundary correspondence~\cite{Volovikbook,Nakaharabook,RMP}, the number of gapless (valley) Dirac edge modes is identical to the change in (valley) Chern number\cite{TKNN,ZhangMacDonaldMele13} 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 $\vec{\Omega}$ over the Brillouin zone\cite{Tse11}
\begin{equation} \label{chernintegral}
    \mathrm{Ch}_1=\frac{1}{2\pi}\sum_n \int \dd^2 k (\vb{\Omega}_n)_z,
\end{equation}
where the sum is taken over the two occupied bands $n$ below the Fermi level at 0 energy and
\begin{equation} \label{berrycurvature}
   \vb{\Omega}_n=i \bra{\frac{\partial u_n}{\partial \vec{k}}} \cross \ket{\frac{\partial u_n}{\partial \vec{k}}}
\end{equation}
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.\cite{Xiao10} However, the Berry curvature is localized near the $K/K'$ 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.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.47\textwidth]{berryberry.png}
    \caption{Berry curvature around the $K/K'$ point (left/right). At the $K'$ point, the sign of the Berry curvature is reversed, but the magnitude is identical.}
    \label{fig:berry_curvature}
\end{figure}

In the non-interacting limit, there is an equivalent formulation of the Chern number given in terms of the Green's function.\cite{MartinBlanterMorpurgo08, Wang12, Tse11, Volovikbook} Here, the Green's function takes on the form $G = (\tilde{E} - H(\vb{k}))^{-1}$
\begin{equation}
    \mathrm{Ch}_{1}=\frac{1}{24 \pi^2}\int\dd{\tilde{E}}\dd^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]
\end{equation}
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(\vec{k})$ rather than of the gauge-dependent normalized eigenstates as in \eqref{berrycurvature}. Since $\partial_0 G^{-1}=1$, we have that
\begin{equation}
    \mathrm{Ch}_{1}=\frac{1}{8 \pi^2}\int \dd{\tilde{E}} \dd^2 k\, \mathrm{Tr}\qty[G[G\partial_{1}G^{-1},G\partial_{2}G^{-1}]].
\end{equation}
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
\begin{equation}
    \mathrm{Ch}_1=-\nu\operatorname{sgn} V.
\end{equation}
At the interface where $V$ changes sign, the difference in the Chern number is $-2/+2$ at the $K/K'$ point, in agreement with our previous observations of two left/right moving modes at the $K/K'$ point as obtained from the band structure.

\subsection{Symmetry protected topology}
The spinless bilayer graphene model \eqref{TBH} 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~\cite{AltlandZirnbauer97} 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~\ref{sec:chiral} and figure~\ref{fig:evals_screw}) 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 \begin{align}3m1'=\mathbb{Z}_2^T\times C_{3v}.\label{mpg}\end{align} In this subsection, we discuss the symmetry-protected topology.

\begin{figure}[htbp]
\centering\includegraphics[width=0.3\textwidth]{halfBZ.pdf}
\caption{The shaded triangular region represents the half Brillouin zone $\mathrm{BZ}^{1/2}_+$ that is used as the integration domain of the valley Chern number around $K$. The dashed lines represent the mirror-symmetric lines $L_1=\Gamma M_1\Gamma$, $L_2=\Gamma M_2\Gamma$ and $L_3=\Gamma M_3\Gamma$.}\label{fig:halfBZ}
\end{figure}

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 \eqref{chernintegral}, 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~\ref{fig:halfBZ}), 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~\ref{fig:halfBZ}) 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'$. 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 \begin{align}v\mathrm{Ch}_1=\frac{1}{2\pi}\sum_{E_n({\bf k})<0}\int_{\mathrm{BZ}^{1/2}_+}d^2k\Omega_n({\bf k})\label{vChern}\end{align} 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 $\vb{\Omega}_n=\nabla_{\bf k}\times\boldsymbol{\alpha}_n$, where the Berry connection is \begin{align}\boldsymbol{\alpha}_n=i\left\langle u_n|\nabla_{\bf k}u_n\right\rangle.\end{align} From the Stokes' theorem, the valley Chern number is \begin{align}v\mathrm{Ch}_1=\frac{1}{2\pi}\sum_{E_n({\bf k})<0}\oint_{\partial\mathrm{BZ}^{1/2}_+}d{\bf k}\cdot\boldsymbol{\alpha}_n({\bf k}),\label{vChern2}\end{align} 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) \begin{align}P_j=\frac{1}{2\pi}\oint_{L_j}d{\bf k}\cdot\boldsymbol{\alpha}_n({\bf k})\end{align} can take any real value.
Time-reversal symmetry requires them to take integral values~\cite{FuKane06,QiHughesZhang08,RMP}. 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 \eqref{vChern2} therefore has integral value. Lastly, since the Berry curvature is gauge invariant, so is the valley Chern number \eqref{vChern} and the sum of polarizations along $L_1$, $L_2$ and $L_3$ in \eqref{vChern2}.

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 \begin{align}3m'=\{1,r,r^2,T\mu_x,Tr\mu_x,Tr^2\mu_x\}\label{mpg0}\end{align} inside the full magnetic point group $3m1'$ in \eqref{mpg}. 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 \eqref{mpg0} 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'$ in \eqref{mpg0}. 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.

\section{Strong Electric and Magnetic Field} \label{sec:strong}
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 $\vb{A} = -By\vu{x}$ is chosen with $y = 0$ set at the interface.

Introducing a magnetic field requires the Peierls substitution in the hopping parameter
\begin{equation}
t\rightarrow t e^{i\theta}, \quad \theta\equiv\frac{2\pi}{\phi_0}\int_{r}^{r'}\vec{A}\vdot\dd{\vb{l}}
\label{eq:Peierls}
\end{equation}
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 (\ref{eq:landau}), $E_2-E_1=2(\hbar v_F)^2 e B \sqrt{2}/\hbar u\approx (0.005$ eV/T$)B$, so for $V\approx0.1$ eV, a comparable magnetic energy scale requires a field of about $B\approx20$ 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
\begin{equation}
    B = \frac{\phi_{0}}{\text{Area}}\frac{p}{q}
\end{equation}
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 \eqref{eq:Peierls} 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'$ point.

\begin{figure}
    \centering
    \includegraphics[width=0.47\textwidth]{bandplots2.png}
    \caption{Electronic band structure with $B = 0$ T, $V = 0.1$ eV (left); $B = 20$ T, $V = 0.1$ eV (right), with $E_{f} = 0$ eV for both figures. The pairs of lines connecting the valence and conducting bands near the $K$ point are the chiral modes living on the interface (circled in the right figure), with a Fermi velocity given by the slope of the line. Additional lines crossing the Fermi level are localized states on the edge of the graphene strip, a result of using an open geometry in our calculations. Landau levels are observable on the right when $B \neq 0$.}
    \label{fig:bands}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=.5\textwidth]{WavefnRangeofB.png}\hfill
    \includegraphics[width=.5\textwidth]{WavefnLargeB.png}
    \caption{The shape of the wavefunction changes with the introduction of a magnetic field. The presence of a hump far from the interface becomes more prominent with a larger magnetic field. The side the hump forms on is determined by the sign of $B$. Here, $V=0.1$ eV and $E_f = 0$ eV.} \label{fig:ebstates}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=.45\textwidth]{WavefnEf.png}
    \caption{The peak of the wavefunction shifts with a nonzero Fermi energy when $B\neq0$. Here, $V=0.1$ eV and $B = 20$ T.} \label{fig:ebstates2}
\end{figure}

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'$ point, as explained in section \ref{sec:chiral}. Plotting the wavefunctions corresponding to the chiral modes confirms that these states are exponentially localized at the domain wall as given by equation \ref{eq:psiv}. 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 \ref{fig:std_vs_V}.

When both $B$ and $V$ are non-negligible, both flat bands (Landau levels) and interface modes can be observed in the band structure (Figure \ref{fig:bands}). 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 \cite{ZareniaPereiraFariasPeeters11}. The chiral modes are still localized near the interface, but their mean position $\ev{y}$ is not located exactly at the interface due to the formation of a secondary peak in the wavefunction (Figure \ref{fig:ebstates}). 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 \ref{fig:ebstates2}). 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 \ref{fig:vsBandVfixedEf}). Specifically, $\ev{y}(B)=-\ev{y}(-B)$. This shifting of $\ev{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 \ref{fig:vsBandEffixedV}). 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 \ref{fig:vsVandEffixedB}).

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{AvgY_B_V.png}\hfill
    \caption{Overall shift of the wavefunction from the potential interface at $E_f = 0$.}
    \label{fig:vsBandVfixedEf}
\end{figure}

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{AvgY_B_Ef.png}\hfill
    \caption{Overall shift of the wavefunction from the potential interface at $V=0.1$ eV.}
    \label{fig:vsBandEffixedV}
\end{figure}

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{AvgY_V_Ef.png}\hfill
    \caption{Overall shift of the wavefunction from the potential interface at $B=20$ T.}
    \label{fig:vsVandEffixedB}
\end{figure}

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 \ref{fig:vsVfixedBandEf}). Figure \ref{fig:vsBfixedVandEf} 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 \ref{eq:yB}. 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 \ref{fig:vsEffixedVandB}). The dependence of $\ev{y}$ on the Fermi energy is to be expected because $\ev{y}$ depends linearly on the Fermi momentum by equation \ref{eq:yB}. This differs from the $B=0$ case, where $\ev{y}=0$ regardless of the Fermi energy.

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{AvgY_V.png}\hfill
    \caption{Overall shift of the wavefunction from the potential interface for $B=20$ T and $E_f=0$ eV. }
    \label{fig:vsVfixedBandEf}
\end{figure}

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{AvgY_B.png}\hfill
    \caption{Overall shift and localization length of the wavefunction from the potential interface and standard deviation at $K$ for $E_{f} = 0$ eV and $V = 0.1$ eV.}
    \label{fig:vsBfixedVandEf}
\end{figure}

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{AvgY_Ef.png}\hfill
    \caption{Layer shifting for the two states near the $K$ point. $B = 20$ T, $V = 0.1$ eV.}
    \label{fig:vsEffixedVandB}
\end{figure}

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 \ref{fig:layer_pref}). 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 (\ref{eq:pbot}).

\begin{figure}
    \centering
    \includegraphics[width=.5\textwidth]{Ptop2.png}
    \caption{Probability for an electron to occupy the top layer.}
    \label{fig:layer_pref}
\end{figure}

\begin{figure}
    \centering

    \includegraphics[width=.475\textwidth]{DiagramsKandKp.png}
    \caption{Electron propagation with $B=0,E_f=0$ (top), $B>0,E_f=0$ (middle), and $B>0,E_f>0$ (bottom) for electrons near the $K$ point. The size of the arrow is proportional to the probability of an electron to occupy that state.}
    \label{fig:electronpropagation}
\end{figure}

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
\begin{equation}
    \psi_{BV}(k_x)\approx \alpha \psi_{V}+\beta_1 \psi_{B,1}+\beta_2 \psi_{B,2}
\end{equation}
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.
\small{
\begin{equation}
\begin{pmatrix}
1 & \bra{\psi_V}\ket{\psi_{B,1}} & \bra{\psi_V}\ket{\psi_{B,2}}\\
\bra{\psi_{B,1}}\ket{\psi_{V}} & 1 & 0\\
\bra{\psi_{B,2}}\ket{\psi_{V}} & 0 & 1
\end{pmatrix}
\begin{pmatrix}
\alpha \\ \beta_{1} \\ \beta_{2}
\end{pmatrix}
=
\begin{pmatrix}
\braket{\psi_V}{\psi} \\ \braket{\psi_{B,1}}{\psi} \\ \braket{\psi_{B,2}}{\psi}
\end{pmatrix}
\label{eq:bvapprox}
\end{equation}}

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 \ref{fig:approx}). 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 (\ref{eqn:dkx}). 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 \textit{approaches} the interface with increasing $B$ consistent with equation \ref{eq:yB}, 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 \textit{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 \ref{fig:ebstates} 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$.

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{MagnitudeOfProjection.png}\hfill
    \caption{Modulus of the projection of our approximated state $\psi_{BV}$ onto the actual state $\psi$. This is a measure of the quality of our approximation.} \label{fig:approx}
\end{figure}

\begin{figure
    \centering
    \includegraphics[width=\linewidth]{alphabeta.png}\hfill
    \caption{Modulus of the electric component $\alpha$ and magnetic component $\beta$ of $\psi_{BV}$, where $\beta$ is defined by $|\beta|^2=|\beta_1|^2+|\beta_2|^2$. The components change smoothly with varying $B$ and $V$.}
\end{figure}

\section{Discussion and Conclusion} \label{sec:con}
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~\cite{Li16,Li18}. With no magnetic field, the counter-propagating modes at $K\pm \delta k$ and $K'\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$, $|\braket{\psi^{*}(K-\delta k)}{\psi(K'+\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'\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 \ref{eq:bvapprox}, 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.

\acknowledgments
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.


text:
\section{Introduction}

Following Hart and Nisan work in \cite{hart2012approximate} we are looking at the scenario of a single seller selling multiple items to a single buyer. This paper is limited to the case of two items with independently distributed values and an additive valuation. The distributions (given by a cumulative distribution $F$) are known to the seller but not the chosen valuations.

For one item the problem is completely solved by Myerson's classic result \cite{myerson1981optimal}. If the seller offer to sell it for a price $p$ then the probability that the buyer will buy is $1 - F(p)$, and the revenue will be $p\cdot(1 - F(p))$. The seller will choose a price $p^*$ that maximizes this expression. Myerson's characterization of optimal auctions concludes that the take-it-or-leave-it offer at the above price $p^*$ yields the optimal revenue among all mechanisms. Myerson's result also applies when there are multiple buyers, in which case $p^*$ would be the reserve price in a second price auction.

Hart and Nisan showed that selling each item separately using the optimal Myerson's price, gains at least half of the revenue attainable by optimal auction for two items.

In addition, they showed that when the items valuations are drown from two i.i.d distribution the ratio between the revenues is tighter and stand on at least $73\%$.

We show that on the other hand, when the two distributions of the values correspond to two different revenues (when selling each item separately) the revenue ratio is again grater then half. Thus, covering another set of distributions.

\subsection*{Notations and Preliminaries}
\subsubsection*{Mechanisms}
\begin{itemize}
\item $x=(x_1,...,x_k)$ - buyer valuation getting each item.
\item $q(x)=(q_1,...,q_k)$ - denotes the probability
that item i is allocated to the buyer when his valuation is $x$
\item $s(x)=(s_1,...,s_k)$ - denotes the expected payment that the seller receives from the buyer when the buyer's valuation is $x$
\item $b(x)=(b_1,...,b_k)$ - denotes the utility of the buyer when his valuation is $x$,
i.e., $b(x) = x\cdot q(x) - s(x)$
\item IR - Individually Rational - $b(x)\geq 0$ for all $x$
\item IC - Incentive Compatible - for all $x,x'$: $x\cdot q(x) - s(x)\geq x\cdot q(x') - s(x')$
\end{itemize}

\subsubsection*{Revenue}

For a cumulative distribution $F$ on \(\mathbb{R}^{k}_{+}\) (for \(k \geq 1\)), we consider the optimal revenue
obtainable from selling $k$ items to a (single, additive) buyer whose valuation for the $k$
items is jointly distributed according to $F$:
\begin{itemize}
\item \(Rev(F)\) is the maximal revenue obtainable by any incentive compatible
and individually rational mechanism.
\item \(SRev(F)\) is the maximal revenue obtainable by selling each item separately.
\end{itemize}

\section{Result}
\begin{theorem}For every one-dimensional distributions $F_1$ and $F_2$,
$$\left(2-\frac{\alpha - 1 - ln(\alpha)}{1+\alpha}\right)SRev(F_1\times F_2) \geq Rev(F_1 \times F_2)$$
where $\alpha=\max{\{\frac{Rev(F_1)}{Rev(F_2)},\frac{Rev(F_2)}{Rev(F_1)}\}}$
\end{theorem}

From now on we denote $R_i=Rev(F_i)$. w.l.o.g. we assume that $R_1\geq R_2$ and denote $\frac{R_1}{R_2}$ as $\ratio$ ($\geq 1$).

\begin{lemma}\label{r1_r2_min_is_grater_then_rev} Let $A,\ B$ be two items with value distributions $X$ and $Y$ associated with $SRev$ of $R_1,\ R_2$ respectively. $R_1 + R_2 + \mathbb{E}(\min\{X,Y\})\geq \mathbb{E}s(X,Y) $
\end{lemma}
\begin{proof}
we use a similar technique to the one used in appendix A in \cite{hart2012approximate}. Take any IC and IR mechanism $(q,s)$. We will split its expected revenue into two parts, according to which one of $X$ and $Y$ is maximal: $\mathbb{E}(s(X,Y)) \leq \mathbb{E}(s(X,Y)\mathbbm{1}_{X\geq Y} ) + \mathbb{E}(s(X,Y )\mathbbm{1}_{Y\geq X})$

for every fixed value $y$ of $Y$ define a mechanism $(\tilde{q}^y, \tilde{s}^y)$ for $X$ by $\tilde{q}^y(x) := q_1(x,y)$ and $\tilde{s}^y(x) :=s(x,y)-yq_2(x,y)$ for every $x$ (so the buyer's payoff remains the same: $\tilde{b}^y(x) = b(x,y)$).

The mechanism $(\tilde{q}^y, \tilde{s}^y)$ is IC and IR for $X$, since $(q,s)$ was IC and IR for $(X,Y)$ (only the IC constraints with $y$ fixed, i.e., $(x′,y)$ vs. $(x,y)$, matter). We do the same for $(\tilde{q}^x, \tilde{s}^x)$.
Therefore,
\begin{equation} \label{item1}
R_1 \geq \mathbb{E}_{x\sim X}(S_1^y(x)) \geq \mathbb{E}(S_1^y(x)\mathbbm{1}_{x\geq y}) \geq \mathbb{E}((s(x,y)-y)\mathbbm{1}_{X\geq y})
\end{equation}
\begin{equation} \label{item2}
R_2 \geq \mathbb{E}_{y\sim Y}(S_2^x(y)) \geq \mathbb{E}(S_2^x(y)\mathbbm{1}_{y\geq x}) \geq \mathbb{E}((s(x,y)-x)\mathbbm{1}_{Y\geq x})
\end{equation}

summing equations \ref{item1}, \ref{item2} and taking expectation over $X, Y$ we get:
\begin{equation}
R_1 + R_2 \geq \mathbb{E}s(X,Y) - \mathbb{E}(Y\mathbbm{1}_{X\geq Y} + X\mathbbm{1}_{Y\geq X}) =  \mathbb{E}s(X,Y) - \mathbb{E}(\min\{X,Y\})
\end{equation}

\end{proof}

\begin{lemma}\label{r1_much_bigger_r2}
$\mathbb{E}(\min\{X,Y\})\leq (2+\ln(\ratio))R_2$
\end{lemma}
\begin{proof}

\[\mathbb{E}(\min\{X,Y\}) = \int_{0}^{\infty}{\bar{F_1}(u)\bar{F_2}(u)du}\]
\[ \leq \int_{0}^{R_2}{1\cdot du} + \int_{R_2}^{R_1}{\frac{R_2}{u} du} + \int_{R_1}^{\infty}{\frac{R_1R_2}{u^2} du}\ = R_2 + R_2\ln\left(\frac{R_1}{R_2}\right) + \frac{R_1R_2}{R_1} = R_2(2+\ln\left(\frac{R_1}{R_2}\right))\]
\end{proof}

combining the lemmas we get that
\[\mathbb{E}s(X,Y)\leq R_1 + (3+\ln(\frac{R_1}{R_2}))R_2\]

\[ = \left(2-\frac{\alpha - 1 - ln(\alpha)}{1+\alpha}\right)(R_1+R_2)\]

thus, finishing the proof of Theorem 1.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
The lightest chemical element in which beta-delayed proton ($\beta$p)
emission has been observed so far is carbon, in the decay of
$^9$C. This paper explores the possibility of seeing the decay mode
also in beryllium and boron, more specifically in the decays of
$^{11}$Be and $^8$B. These nuclei are the lightest (not counting the
deuteron) one-neutron and one-proton halo nuclei, respectively, and
the possible $\beta$p decays are intimately connected to this
structure. There are two reasons for this connection, the first being
the factorization of the wavefunction into a halo and a core part that
suggests ``independent decays'' of the two parts \cite{Nil00} which
naturally leads to final states with a continuum proton; we shall
discuss this model in the next section. The second connection is
specific for $\beta^-$ decays where one can derive the general
expression \cite{Jon01}
\begin{equation}
 Q_{\beta p}(^AZ) = (m_n-m_H)c^2 - S_n(^AZ) =\mathrm{782\, keV} - S_n(^AZ)
\end{equation}
from which it is seen that the $\beta$p decay only occurs for
nuclei with very low neutron separation energy. Similar relations
exist for beta-delayed deuteron and triton emission that are known to
take place close to the neutron dripline \cite{Nil00,Pfu11}. So
far $\beta$p decays have only been observed following $\beta^+$ or EC
decays. General recent reviews of these processes can be found in
\cite{Pfu11,Bla08}.

The following section discusses the characteristics of the proposed
decays in more detail including the estimated order of magnitude of
the decay rate. The expected branching ratios are quite low and some
general considerations are made on how they could be measured.
Sections 3 and 4 give details on first measurements of
the decays in $^{11}$Be and $^8$B from which upper limits on
the branching ratios can be given. The final section gives our
conclusions and the prospects of identifying the two decays.

\section{Expected rates}
A crucial feature of halo states is their intrinsic clustering so that
(a major part of) the wavefunction will factorize into a halo part and
a core part \cite{Rii00,Jen04}. Formally, the beta-decay can then be
thought of in terms of separate decays of the core (c) and halo (h)
parts \cite{Nil00}
\begin{equation}  \label{eq:fac}
  \mathcal{O}_{\beta} (|c\rangle |h\rangle) = (\mathcal{O}_{\beta}
  |c\rangle) |h\rangle + |c\rangle (\mathcal{O}_{\beta} |h\rangle)
\end{equation}
but the actual final states will of course in general be
superpositions of the two terms on the right side e.g.\ to ensure they
have proper isospin. If the states on the right hand side are
close to being eigenstates for the final system one would expect a
large resemblance between the decays of the core nucleus and the halo
nucleus. Experimental confirmation of this idea may be found in the
decays of $^{12}$Be and $^{14}$Be, which are quite analogous
\cite{Jon04}.

\subsection{$^{11}$Be($\beta$p)}
The most recent mass tables \cite{mas11} give $Q_{\beta}= 11509.2 \pm
0.5$ keV for the decay of $^{11}$Be, which leaves 281 keV as the
energy window for beta-delayed proton emission.  Different assumptions
on the decay mechanism can give widely different rates \cite{Ang98}.
We shall therefore present several simplified models that lead to
rather straightforward expressions for the decay rate, and start by
considering sequential decays as described in R-matrix theory that can
handle broad intermediate states. However, the classical R-matrix
expressions may not give a sufficiently accurate description of the
decay mechanism. At the other extreme are models that consider decays
to proceed directly to continuum states, as has been assumed for most
calculations of the related $\beta$d decays \cite{Nil00,Pfu11}.

Allowed beta decay of $^{11}$Be will populate $1/2^+$ and
$3/2^+$ levels in $^{11}$B, but since these would decay to $^{10}$Be+p
with s-wave and d-wave protons, respectively, only the $1/2^+$ states
can be expected to give a sizeable contribution. Out of the known
states \cite{Ajz90,Yam11} in the excitation energy range 10.0--12.5
MeV the only one that may have spin and parity $1/2^+$ is a state at
11.444 MeV with width 103 keV that has been seen to decay by emission
of $\alpha$ particles, but where the spin is so far undetermined.
The R-matrix expressions for $\beta$-decay have been derived for
several cases, see \cite{Bar89} and references therein, but we shall
here only consider a crude approximation based on the single-channel
single-level case in which the energy spectrum is:
\begin{equation}
  \frac{w_p}{w_{tot}} = T_{1/2} \frac{g_A^2}{K} f_{\beta}(Q-E)
  B_{GT}  \frac{\Gamma/(2\pi)}{(E-E_0)^2 + \Gamma^2/4}
  \frac{P(E-E_t)}{P(E_0-E_t)} \;,
\end{equation}
where $w$ denote the decay rates for the proton branch and the total
beta decay, $T_{1/2}$ the experimental halflife, $g_A$ and $K$ the
standard weak decay constants \cite{Pfu11,Bla08}, $f_{\beta}$ the beta
decay phase space factor, $B_{GT}$ the reduced Gamow-Teller matrix
element squared, $E_0$ and $\Gamma$ the level energy and width, and $P(E-E_t)$
the proton penetrability ($E_t$ being the proton threshold
energy). Several energy dependent terms in the denominator are here
approximated as constants. To take into account that the level can
decay both by emission of $\alpha$ particles and protons one should
multiply by $\Gamma_p/\Gamma_{tot}$ to get the intensity of the proton
spectrum. This ratio is unknown, but may in the best case (if the
proton width is close to a single-particle width) reach $1/2$. The GT
strength $B_{GT}$ can at most be 3 (the value for a pure neutron to
proton decay). Inserting these maximum values the total branching
ratio becomes $1.0 \cdot 10^{-6}$, but a realistic value could be
lower by an order of magnitude or more. The predicted spectrum is
given in figure \ref{fig:e11be} for the case where the intensity (the
product $B_{GT}\Gamma_p/\Gamma_{tot}$) is a factor 50 below the
maximum, i.e.\ where the branching ratio is $2.0 \cdot 10^{-8}$.

As an intermediate step we consider next the suggested modification of
the R-matrix formalism containing the opposite time ordering
\cite{Bor93}, i.e.\ where the neutron is emitted before beta decaying
into a proton. As argued by Barker \cite{Bar94}, who presents a more
detailed derivation analogous to the way photon emission is treated,
this may be a simple way to effectively include the decays to
continuum states into the R-matrix framework.  In the limit where only
the ``emission before decay'' process is included the decay rate
becomes
\begin{equation}
  \frac{w_p}{w_{tot}} = T_{1/2} \frac{g_A^2}{K} f_{\beta}(Q-E)
  B_{GT}  \frac{1}{\pi} \frac{P(E-E_t)\gamma^2}{(S_n+E)^2} \;,
\end{equation}
where $S_n$ is the neutron separation energy in $^{11}$Be and the
parameter $\gamma^2$ quantifies the neutron emission probability. For
a decaying level it is related to the level width through $\Gamma = 2
P \gamma^2$ and it must be less than the single particle width (the
Wigner limit, here 5.8 MeV). The total branching ratio is proportional
to $\gamma^2$ and is at most $2.5 \cdot 10^{-8}$. The predicted energy
distribution is also plotted in figure \ref{fig:e11be}.

\begin{figure}
 \hspace{2cm}
 \resizebox{0.7\textwidth}{!}{\includegraphics{fig1_be11t.eps}}
\caption{Energy spectra of beta-delayed protons from the decay of
   $^{11}$Be predicted in three different models presented in the text.
   Sequential decay through the 11.444 MeV state in $^{11}$B
   (dashed line), emission before decay (doted line),
   decay to continuum (full line). \label{fig:e11be}}
\end{figure}

Finally, the simple model for decays directly to continuum states
\cite{Rii90} gives
\begin{equation}
  \frac{w_p}{w_{tot}} = T_{1/2} \frac{g_A^2}{K} f_{\beta}(Q-E)
  B_{GT}(E) \frac{P(E-E_t)}{R} \frac{mc^2}{2\pi^2 (\hbar c)^2} \;,
\end{equation}
where $m$ is the reduced mass of the outgoing proton and $R$ is the
channel radius used to evaluate the penetrability. The Gamow-Teller
strength is now explicitly energy dependent since the final state
wavefunction is a continuum state. The expression above results if one
uses a plane wave approximation for the final state when calculating
$B_{GT}(E)$.  If one instead uses a Coulomb wave the effect of the
transmission through the Coulomb potential is included in $B_{GT}(E)$
and one must replace the ratio $P/R$ by the
wavenumber $k$. The total branching ratio is in this latter case $2.3
\cdot 10^{-8}$, the corresponding energy distribution is given in
figure \ref{fig:e11be}.  More realistic calculations will use
distorted waves so that the final state strong interactions are also
included. A recent more sophisticated calculation \cite{Bay11} in a
two-body potential model gives a broad energy spectrum peaking at
0.1--0.2 MeV and a branching ratio of $3.0 \cdot 10^{-8}$. (This
specific calculation includes both Fermi and Gamow-Teller
contributions to the decay rate.)

To summarize, the different models predict somewhat different shapes of
the energy distribution and give a branching ratio that typically is a
few times $10^{-8}$, but could range up to $10^{-6}$. To see the decay
mode an experiment therefore needs to be sensitive down to the
$10^{-8}$ level, while considerably more intensity will be needed to
separate the different models.

\subsection{$^8$B($\beta$p)}
Turning now to $^8$B two facts make the decay directly to continuum
states less likely here: the smaller spatial extent of the halo in this
nucleus and the fact that there is an excited state in $^8$Be within
the $\beta$-decay window (actually the EC-decay window)
that the decay may pass through, namely the
$1^+$ $T=1$ state at 17.640 keV. It is situated 385 keV above the
proton threshold and is know to decay mainly by proton emission \cite{Til04}. If
equation (\ref{eq:fac}) is appropriate, we can estimate the matrix element of
the transition to be the same as for the ground state decay of $^7$Be
into the ground state of $^7$Li (the larger amplitude of the
wavefunction in $^8$B increases the rate by a factor 2.1 \cite{Bam77}
but this is compensated by a factor 1/2 from the isospin
Clebsch-Gordan coefficient squared). Both transitions are electron
capture decays for which the phase space goes as the Q-value squared,
so by scaling the halflife of $^7$Be we obtain an estimate of the
branching ratio of the $\beta$p transition of $2.3 \cdot 10^{-8}$.

Direct support for this estimate can be found from the three-cluster
calculations in \cite{Gri00} that predict a $B_{GT}$ to the 17.64 MeV
state of 1.366 and 1.997 for two different potentials employed. This
is in the same range as the $B_{GT}$ of 1.83 for the $^7$Be ground
state transition.

\subsection{Experimental considerations}
Direct detection of these $\beta$p decays is challenging since one has
to look for particles of low energy, with a low branching ratio, and
where the decay branch has to be identified e.g.\ through particle
identification (PID) techniques. As an alternative one may try to
detect the presence of the daughter nucleus or, more ambitiously, its
growing-in. In the cases in question here the final nucleus is stable
or longlived so one cannot make use of its decay, but must resort to
other methods. One possibility is accelerator mass spectroscopy (AMS)
as explained in the next section.

A low branching ratio by itself does not prohibit direct detection.
As an example, branching ratios well below $10^{-10}$ have been
detected for cluster radioactivity \cite{Pri89}. The challenge is
rather to identify the decay branch uniquely. For proton kinetic
energies below 1 MeV one could employ gas telescopes and use the
energy loss $\Delta E$ versus full energy $E$ to identify the
particle. However, the energy loss for a proton peaks just below 100
keV energy so the $\Delta E$-$E$ curves for protons, deuterons and
tritons will cross in this region, making separation e.g.\ between
protons and tritons (relevant for $^{11}$Be) impossible.
Identification via detection of $E$ and time-of-flight is an
alternative, but will typically have much less efficiency. Separating
protons from the background of energy loss signals from beta particles
could also be a problem.

If the proton can be recorded in coincidence with the recoiling
nucleus this may be used to discriminate against other channels. This
``ratio cut'' method is useful to separate true events from response
tails of events at higher energy, as done in the $\beta\alpha$ decay
of $^{16}$N \cite{Azu94} and needed here for the case of $^8$B. It can
also be used to single out one decay channel among many others, as
done for the $^{11}$Li beta-decay to t+$^8$Li \cite{Mad09} and needed
here for $^{11}$Be. The main challenge with this method is to obtain
efficient detection of the recoil with sufficient accuracy in
energy, in the examples given the energy was higher than for our
$\beta$p decays.

\section{The $^{11}$Be experiment}
Four beta-delayed particle branches are energetically open for
$^{11}$Be, their Q-values are listed in table \ref{tab:QS}. We may
neglect the $\beta$n branch since it would proceed to the $^{10}$B
$3^+$ ground state and require a d-wave neutron with very low
energy. The $\beta$t branch leads to $^8$Be and would therefore give a
three-body final state, this could also be a quite interesting decay
mode. The main background for the $\beta$p branch will be
$\beta\alpha$ whose branching ratio was measured to be $3.1\pm0.4$\%
\cite{Mil82} (a recent determination gave $3.47\pm0.12$\%
\cite{Bus10}).

\begin{table}
 \caption{Q-values and separation energies in $^{11}$Be
 beta-delayed particle decays. \label{tab:QS}}
 \begin{tabular}{lcccc}
 \br
 $x$ & p & n & $\alpha$ & t \\ \mr
 S$_x$($^{11}$B) (keV)$^a$ & $11228.5 \pm 0.4$ & $11454.12 \pm 0.16$ &
 $8664.1 \pm 0.4$ & $11223.5 \pm 0.4$  \\
 $Q_{\beta x}$ (keV)$^a$ & $280.7 \pm 0.3$ & $55.1 \pm 0.5$ &
 $2845.2 \pm 0.2$ & $285.7 \pm 0.2$ \\
 \br
 \end{tabular}
$^a$Mass values from ref. \cite{mas11}
\end{table}

\subsection{The set-up and source production}
Due to the challenges involved in direct detection of the $\beta$p
decays we decided to explore the possibility of indirect
detection. The $\beta$p-daughter nucleus, $^{10}$Be, is radioactive
with a halflife of $1.39 \cdot 10^6$ years and is only present on earth
in minute quantities due to production by cosmic rays.
It is used as a tracer in earth sciences
and procedures have been developed for detecting very small samples of
$^{10}$Be via accelerator mass spectroscopy (AMS), see e.g.\ \cite{Elm87,Mul08}.
Stated briefly, our experiment then consists in collecting a large
number of $^{11}$Be atoms, measure their number by detecting
$\gamma$-radiation from the $^{11}$Be decay and later measure the
number of produced $^{10}$Be atoms via AMS.

The $^{11}$Be source was collected in September 2001 at the ISOLDE
facility at CERN. The $^{11}$Be atoms were produced
by proton bombardement in a Ta target, ionized by laser ionization
\cite{Fed03}, accelerated and transported to the measuring station
where the ion beam passed two sets of collimators before being
implanted into a Be foil (Goodfellow, LS226536 L C, thickness 0.01mm,
purity 99.8+\%).
The average intensity of the collected beam was $3.6\cdot 10^6$
ions/s, the collection took place for close to 40 hours.
The distance from the first collimator to the collection point was 177
mm.
The transmission was optimized with a stable $^{23}$Na beam and a
radioactive $^{27}$Mg beam, the fraction of the activity that
deposited on the collimator was at most 12\%.

\subsection{Source strength}
Most of the decays of $^{11}$Be proceed through states in $^{11}$B that
decay by gamma emission, the most prominent $\gamma$-ray being the
2124.47 keV line with a total branching ratio, including feeding from
higher levels, of $b_{\gamma} = 0.355 \pm 0.018$ \cite{Mil82}.
The gamma detection was done with a Ge-detector
that was placed about 40 cm
from the collection point in the opposite direction of
the collimator. It was further shielded by 3 mm Pb to reduce the total
count rate. Still, the multi channel analyzer used for data taking had a
deadtime of $14\pm 1$\% determined from the ratio of live time and
real time of the data taking system.

The Ge detector was energy and efficiency calibrated \textit{in situ} with
absolutely calibrated sources of $^{60}$Co, $^{137}$Cs and $^{228}$Th
that bracket the energies from 239 keV to 2615 keV.
The total gamma spectrum from $^{11}$Be is shown in figure \ref{fig:g11be}.
The spectrum is dominated by gamma lines from the decay of $^{11}$Be
with a small background mainly from $^{40}$K and the $^{222}$Rn decay
chain.
The 2124 keV line is clearly separated from background and gives,
after correction for detector efficiency and deadtime, a deduced amount
of $^{11}$Be atoms collected of $(5.2 \pm 0.3) \cdot 10^{11}$. As a
cross-check of this number we shall look at two other $\gamma$-lines
recorded from the $^{11}$Be decay. First the 2895 keV line that again
is  nicely separated in our spectra; we deduce an intensity
ratio with respect to the 2124 keV line of $0.241 \pm 0.011$\% in
agreement with the literature value of $0.227 \pm 0.008$\% \cite{Mil82}.
Secondly, we consider the 478 keV line that follows $\beta\alpha$
decays of $^{11}$Be to the first excited state in $^7$Li. The line is
therefore recoil broadened and care must be taken when extracting its
intensity since it furthermore is situated on top of a large background. Our
deduced intensity relative to the 2124 keV transition is $0.75 \pm 0.06$\%.

\begin{figure}
 \hspace{2cm}
 \resizebox{0.8\textwidth}{!}{\includegraphics{fig2_begam.eps}}
\caption{Gamma energy spectra recorded from the decay of
   $^{11}$Be, the number of counts per 0.4 keV is shown versus gamma-ray
   energy. The star marks the 478 keV line. \label{fig:g11be}}
\end{figure}

The literature value for this last ratio can be derived from the decay
scheme in \cite{Alb81} to be 1.13\% (it was not observed directly
since gamma detection was done with a NaI detector). However, as we
shall argue now, this value is most likely wrong.
Alternative determinations of this ratio
measured with Ge detectors can be extracted from experiments on the
decay of $^{11}$Li in which $^{11}$Be is fed. From two such recent
experiments accurate results can be obtained (although not published
previously): an experiment we did at ISOLDE \cite{Fyn04} gave
$0.67 \pm 0.11$\% and an experiment at TRIUMF \cite{Mat09,Sar09} gave
$0.69 \pm 0.07$\%. An unpublished recent $\beta\alpha$ experiment \cite{Bus10}
has revised the branching ratio of alpha particles to the 478 keV
excited state in $^7$Li to be $8.0 \pm 0.4$\% of all alpha decays, a result
based on a detailed analysis of the alpha spectrum but in disagreement
with the literature value of $12.6 \pm 1.2$\% \cite{Alb81}. Using the new
value one derives a 478/2125 ratio of $0.78 \pm 0.06$\%. Since
the three new determinations all are consistent with our
value for the ratio, we conclude that our spectra are internally
consistent which adds confidence to our derived source intensity.

\subsection{Source purity}
The collection foil must of course be free of previous $^{10}$Be
activity, but we also need to worry about contamination coming along
with the $^{11}$Be beam. The most worrisome contaminant is $^{11}$Li
that can be ionized by surface ionization in the hot cavity where
laser ionization of Be takes place. The primary production yield is
lower by a factor at least 2000, but the main decay
branch of $^{11}$Li is beta-delayed one-neutron emission that leads to
$^{10}$Be and suppression of $^{11}$Li is therefore needed.

The mass difference between $^{11}$Be and $^{11}$Li, $\Delta M =
Q(^{11}Li) = 20.551$ MeV \cite{mas11}, is sufficiently big that the
resolving power of ISOLDE's high resolution separator HRS, which is up
to $M/\Delta M = 5000$, can easily separate the two nuclei
that have $M/\Delta M = 500$. A tail of $^{11}$Li activity may
nevertheless be present at the $^{11}$Be mass position, but we have
estimated its magnitude by
looking for the corresponding tails of the nuclei $^9$Li and $^8$Li a
factor $M(^{11}Be)/M(^{11}Li)$ below their nominal setting. In both
cases background prevented a direct identification --- there are no
$\gamma$-rays in either decay, so only the (non-unique)
$\beta$-detection is possible --- but a lower limit of the suppression
factor of $10^3$ was extracted in both cases.
A second source was collected with the separator deliberately set
off-mass by two thirds of the $^{11}$Li-$^{11}$Be mass difference (below
the $^{11}$Be position). The intensity of the $^{11}$Be gamma lines was
here reduced by a factor 500, consistent with the results at the lower masses.

A further reduction can be obtained by blocking collection for the
first 150 ms after the radioactivities are produced by the proton
pulse. Most $^{11}$Li, of halflife 8.5 ms, will have decayed by
then. However, at the end of our collection it was discovered that the
onset of the blocking by mistake was delayed by a few ms allowing up
to 15\% of the $^{11}$Li to leak through. Thus up to $2 \cdot 10^4$
$^{10}$Be could have contaminated the sample through decay of collected $^{11}$Li.

A second possible contaminant could be the BeH molecule $^1$H$^{10}$Be
that would be almost on top of the $^{11}$Be mass since $M/\Delta M =
36600$ for this case. It is quite unlikely that this molecule will be
formed and ionized in sufficient quantities to be a problem, but the
possibility should be kept in mind in future experiments. As will be
seen we do not have indications for any of these contaminants in the
present sample.

\subsection{AMS measurement}

The $^{10}$Be accelerator mass spectrometry (AMS) measurements were
performed at the Tandem Laboratory, Uppsala University, in Sweden
\cite{Ber09}. A dedicated setup based on a NEC 5 MV tandem pelletron
accelerator devoted to high precision and low background detection of
$^{10}$Be was employed. The measurements were carried out according to
the principle presented by Middletone and Klein \cite{Mid86,Mid93}
where ion currents were normalized by $^{17}$O from the $^9$Be$^{17}$O
molecule. The absolute transmission in the AMS-system was determined
with the NIST SRM 4325 standard ($^{10}$Be/$^9$Be= $3.06 \cdot
10^{-11}$ \cite{Mid86}; NIST certificates the ratio to be R$_{SRM} = 2.68 \cdot
10^{-11}$).

The Be-foil used in the investigation was divided into three $10 \times
10$ mm$^2$ samples where one subfoil A was just chemically prepared
and used as a blank value representative for the whole AMS
procedure. Subfoil B and C were both mounted in two aluminium frames
with a circular aperture (diameter 8 mm). Both were installed in the
ISOLDE facility. There was no implantation in subfoil B while
subfoil C was used as the catcher for $^{11}$Be.

All three subfoils were chemically prepared by dissolution in HCl
followed by precipitation to Be(OH)$_2$ by adding NH$_3$ which in turn
was converted to BeO. Finally the BeO was mixed with Nb as a binder
before the AMS-analysis in the accelerator.

\begin{table}
 \caption{Results of AMS measurements. \label{tab:AMS}}
 \begin{tabular}{cccc}
 \br
 Subfoil & Mass (mg) & R=R$_{foil}$/R$_{SRM}$ & $^{10}$Be ($10^6$ atoms/mg)  \\ \mr
 A (blank) & 3.04 & 0.00149(12) & 2.67(22) \\
 B (background) & 3.65 & 0.00152(11) & 2.73(20) \\
 C (catcher) & 4.22 & 0.00169(13) & 3.03(23) \\
 \br
\end{tabular}
\end{table}

A summary of the results is given in table \ref{tab:AMS}.
No significant difference is observed for the blank and background
foils which indicates no additional $^{10}$Be contamination in the
handling in ISOLDE and during transportation. It is, however,
noteworthy to compare the background values for the foil used with
what is obtained from our normal Be carrier which for 1 mg sample
gives a R $= 0.0005(1)$ corresponding to a factor two lower $^{10}$Be content.

If we use subfoil B as a background to subfoil C the $^{10}$Be signal
can be estimated to be $(1.28 \pm 1.29) \cdot10^6$ atoms. With a total
number of $^{11}$Be nuclei in the sample of $(5.2 \pm 0.3) \cdot
10^{11}$ this leads to a deduced branching ratio of $(2.5 \pm 2.5)
\cdot 10^{-6}$ where the uncertainty is on the one sigma level. There
is therefore no positive evidence for the $^{11}$Be $\beta$p decay,
only an upper limit of the order of the most optimistic theoretical
branching ratios.

\section{The $^8$B experiment}
The $^8$B experiment was carried out at the IGISOL (Ion Guide Isotope
Separator On-Line) separator at the Accelerator Laboratory of
University of Jyv\"{a}skyl\"{a}. The main aim of this experiment was
to determine precisely the neutrino spectrum from the decay of $^8$B
through analysis of the beta-delayed alpha spectrum. The result of
this part of the experiment along with a detailed description of the
experimental procedures is described in \cite{Kir10,Kir11}; the set-up
consisted of four DSSSD detectors each backed by a thick Si detector
in close geometry around the collection point. As a byproduct
the high statistics of the experiment allowed to identify for the
first time an electron capture branch in the decay of $^8$B, namely to
the $2^+$ 16.922 MeV state.  We are here concerned with the possible
electron capture branch to the even higher-lying 17.640 MeV state.

Since the $^8$B experiment recorded a total of 16 million decays, the
expected number of proton detections is 0.36 and we can only report an
upper limit.  The 17.640 keV state will decay by emitting a 337 keV
proton and a 48 keV recoiling $^7$Li ion that could not be detected. The
main background from the decay will be $\alpha$-$\alpha$ coincidences
at low energy or positrons that happen to deposit around 340 keV in
the detector. The latter background can be suppressed by requiring
anti-coincidence with the backing detectors. The former background
needs to be addressed indirectly since the set-up did not include
particle identification. One can suppress the $\alpha$ particles
by doing anti-coincidence with the opposite DSSSD, but to do this
efficiently one needs to restrict the solid angle and only consider
events close to the center of each DSSSD since beta-recoil can cause
the $\alpha$ coincidences to deviate from a linear geometry. This
reduces the effective solid angle by a factor about 9. A small
background nevertheless remains, see figure \ref{fig:sp8b},
and the 95\% confidence level upper limit of the number of counts in a
337 keV peak turns out to be 18. As seen in the figure there are
indications for a peak at 349 keV; however, the uncertainty on the
energy scale is only 7 keV since there were accurately known
calibration lines from the decay of $^{23}$Al close by \cite{Kir11a}.
The 95\% confidence level upper limit on the number of counts in
the 349 keV peak is 45. Combing all factors gives an
upper limit of $2.6 \cdot 10^{-5}$ on the $\beta$p branching ratio at
95\% confidence level. To reach the range of predicted intensity the
experiment needs to improve by a factor 1000. One may gain a factor of
ten by including particle identification, which will reduce the
background in the interesting energy range. Another factor of ten may
be obtainable by increasing the effective solid angle, but one still would
need a factor of ten higher overall production of $^8$B ions.

\begin{figure}[ht]
 \hspace{2cm}
 \resizebox{0.7\textwidth}{!}{\includegraphics{fig3_b8.eps}}
\caption{The low-energy spectrum from the decay of $^8$B after
   background suppression. The red curve shows the two sigma
   upper limit from of a fit  where the peak position was allowed
   to vary within uncertainties.
\label{fig:sp8b}}
\end{figure}

\section{Outlook}
The theoretical estimates in section 2 showed that beta-delayed proton
emission in the decays of $^{11}$Be and $^8$B is likely to take place
with a branching ratio that can be expected to be a few times
$10^{-8}$. The energies of the emitted proton and the recoiling
nucleus will be small. These two factors make detection of the
processes challenging.

We have carried out first experiments to look for these decay branches
and have obtained upper limits on the branching ratios. These do not
yet reach the level of the theoretical predictions, improvements are
needed by factors of at least $10^2$ and $10^3$ for the two cases.
The continuing improvements of radioactive beam facilities may give
increases in the yields that will allow experiments to reach the
needed level of sensitivity, but other possible improvements should
also be considered. As explained above for $^8$B changes to a more
dedicated set-up should be considered. For the case of $^{11}$Be
the AMS measurement can be considerably improved if the amount of
stable Be in the sample is reduced.  This can be achieved by using a different
catcher foil (e.g.\ Cu, Au or Nb) and after adding a small amount of Be
carrier ($< 0.1$ mg with low $^{10}$Be content) do a refined chemical
separation \cite{Ber09a}. An enhanced sensitivity of two orders of magnitude can
be expected from such improved sample handling.

Once the decay modes have been established and the branching ratio
securely determined, the next important step will be to determine the
energy distribution of the emitted protons. This will give a more
sensitive test of the theoretical calculations and may in particular
show whether our hypothesis of a direct relation between the $\beta$p
decay mode and the halo structure in these two nuclei is correct.

\bigskip

\ack
We would like to thank I. Mukha and L. Grigorenko for
discussions and F. Sarazin for communicating non-published results on
the $^{11}$Li decay.
This work was partly financed by the Spanish Research funding agency under
project CICYT FPA2007-62170 and FPA2009-07387 and by MINECO through
projects FPA2010-17142 and CPAN CSD-2007-00042.

\section*{References}