Quantum Anomalous Parity Hall Effect in Magnetically Disordered Topological Insulator Films

TL;DR

This paper demonstrates a novel quantum anomalous parity Hall phase in magnetically disordered topological insulator films, enabling electrical control of Majorana modes without external magnetic fields.

Contribution

It introduces the quantum anomalous parity Hall phase arising from magnetic disorder and shows how to manipulate Majorana zero modes electrically in this system.

Findings

Disorder-induced helical edge modes protected by reflection symmetry.

Electrical creation and manipulation of Majorana zero modes.

Potential for zero-field topological quantum computing.

Abstract

In magnetically doped thin-film topological insulators, aligning the magnetic moments generates a quantum anomalous Hall phase supporting a single chiral edge state. We show that as the system de-magnetizes, disorder from randomly oriented magnetic moments can produce a `quantum anomalous parity Hall' phase with \emph{helical} edge modes protected by a unitary reflection symmetry. We further show that introducing superconductivity, combined with selective breaking of reflection symmetry by a gate, allows for creation and manipulation of Majorana zero modes via purely electrical means and at zero applied magnetic field.

Equations

$$\mathcal{H}=v(k_x{\sigma }_x+k_y{\sigma }_y){\tau }_z+t(k){\tau }_x+m(\mathbf{r}){\sigma }_z.$$

$$⟨m(\mathbf{r})⟩=\bar{m},⟨m(\mathbf{r})m({\mathbf{r}}^{\prime })⟩-{\bar{m}}^2=\delta m^{2}K(\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$$h_{\pm }=-v(k_x{\sigma }_x+k_y{\sigma }_y)+[m(\mathbf{r})\pm t(k)]{\sigma }_z.$$

$$C_{\pm }=\left[\mathrm{sgn}(\bar{m}\pm t_0)\mp \mathrm{sgn}(t_2)\right]/2.$$

$$\begin{array}{rl}{\Sigma }_{\pm }(\omega ,\mathbf{k})=\delta m^2\int & \frac{{\mathrm{d}}^{2}q}{(2\pi {)}^2}\tilde{K}(\mathbf{q}){\sigma }_z[\omega -{\bar{h}}_{\pm }(\mathbf{k}-\mathbf{q})\\ & {-{\Sigma }_{\pm }(\omega ,\mathbf{k}-\mathbf{q})]}^{-1}{\sigma }_z.\end{array}$$

$${\Sigma }_{\pm }(0,\mathbf{k})={\Sigma }_{\pm }^{(1)}(k_x{\sigma }_x+k_y{\sigma }_y)+({\Sigma }_{\pm }^{(0)}+{\Sigma }_{\pm }^{(2)}{\mathbf{k}}^2){\sigma }_z,$$

$$\delta m_{\pm }^{\mathrm{c}}={\left[\frac{t_{2,\pm }^{\prime }\left(t_0\pm \bar{m}\right)}{{\xi }^2\ln \left\{1+{\left[t_{2,\pm }^{\prime }/\left(v_{\pm }^{\prime }\xi \right)\right]}^2\right\}}\right]}^{1/2},$$

$$\begin{split}H&=\sum_{n_{x}=1}^{L_{x}}\sum_{n_{y}=1}^{L_{y}}\Big{\{}c_{\boldsymbol{n}}^{\dagger}\left(t_{0}^{\prime}\tau_{x}+m_{\boldsymbol{n}}\sigma_{z}\right)c_{\boldsymbol{n}}\\ +&\sum_{\boldsymbol{d}=\hat{x},\hat{y}}\left[c_{\boldsymbol{n}}^{\dagger}\left(i\frac{v}{2}\boldsymbol{\sigma}\cdot\boldsymbol{d}\tau_{z}-t_{2}\tau_{x}\right)c_{\boldsymbol{n}+\boldsymbol{d}}+h.c.\right]\Big{\}},\end{split}$$

$$H=\sum _{n_x=1}^{L_x}\left\{{\vec{\psi }}_{n_x}^{†}{h}_{n_x}{\vec{\psi }}_{n_x}+\left[{\vec{\psi }}_{n_x}^{†}V{\vec{\psi }}_{n_x+1}+h.c.\right]\right\},$$

$$r(\omega )=1-2\pi i{\rho }_{\mathrm{R}}{V}^{†}[G_{L_x}^{-1}(\omega )+i\pi {\rho }_{\mathrm{R}}V{V}^{†}{]}^{-1}V,$$

$$G_{n_x}(\omega )=[\omega -h_{n_x}-V^{†}{G}_{n_x-1}V{]}^{-1}.$$

$$G\equiv \frac{\mathrm{d}I}{\mathrm{d}V}=\frac{e^2}{h}\mathrm{Tr}[t^{†}t]=\frac{e^2}{h}\mathrm{Tr}[1-r^{†}r],$$

$$H=\frac{1}{2}\sum _{n_x=1}^{L_x}\left\{{\vec{\Psi }}_{n_x}^{†}{\tilde{h}}_{n_x}{\vec{\Psi }}_{n_x}+\left[{\vec{\Psi }}_{n_x}^{†}\tilde{V}{\vec{\Psi }}_{n_x+1}+h.c.\right]\right\},$$

$${\Sigma }_{+}(k,\omega )={\Sigma }_0(k,\omega )+{\Sigma }_x(k,\omega ){\sigma }_x+{\Sigma }_y(k,\omega ){\sigma }_y+{\Sigma }_z(k,\omega ){\sigma }_z.$$

$$h_{+}^{\mathrm{eff}}(\mathbf{k})=-v{k}_x{\sigma }_x-v{k}_y{\sigma }_y+\left[\bar{m}+t_0+t_2{\mathbf{k}}^2\right]{\sigma }_z+{\Sigma }_{+}(\mathbf{k},0).$$

$${\Sigma }_0(\mathbf{k})=\frac{\delta m^2}{(2\pi {)}^2}\int \frac{{\mathrm{d}}^{2}q\tilde{K}(\mathbf{q}){\Sigma }_0(\mathbf{k}-\mathbf{q})}{{\sum }_{i=x,y}{\left[v(k_i-q_i)-{\Sigma }_i(\mathbf{k}-\mathbf{q})\right]}^2+{\left[\bar{m}+t_0+t_2(\mathbf{k}-\mathbf{q}{)}^2+{\Sigma }_z(\mathbf{k}-\mathbf{q})\right]}^2-{\Sigma }_0^2(\mathbf{k}-\mathbf{q})},$$

$${\Sigma }_{x,y}(\mathbf{k})=\frac{\delta m^2}{(2\pi {)}^2}\int \frac{{\mathrm{d}}^{2}q\tilde{K}(\mathbf{q})\left[-v(k_{x,y}-q_{x,y})+{\Sigma }_{x,y}(\mathbf{k}-\mathbf{q})\right]}{{\sum }_{i=x,y}{\left[v(k_i-q_i)-{\Sigma }_i(\mathbf{k}-\mathbf{q})\right]}^2+{\left[\bar{m}+t_0+t_2(\mathbf{k}-\mathbf{q}{)}^2+{\Sigma }_z(\mathbf{k}-\mathbf{q})\right]}^2-{\Sigma }_0^2(\mathbf{k}-\mathbf{q})},$$

$${\Sigma }_z(\mathbf{k})=-\frac{\delta m^2}{(2\pi {)}^2}\int \frac{{\mathrm{d}}^{2}q\tilde{K}(\mathbf{q})\left[\bar{m}+t_0+t_2(\mathbf{k}-\mathbf{q}{)}^2+{\Sigma }_z(\mathbf{k}-\mathbf{q})\right]}{{\sum }_{i=x,y}{\left[v(k_i-q_i)-{\Sigma }_i(\mathbf{k}-\mathbf{q})\right]}^2+{\left[\bar{m}+t_0+t_2(\mathbf{k}-\mathbf{q}{)}^2+{\Sigma }_z(\mathbf{k}-\mathbf{q})\right]}^2-{\Sigma }_0^2(\mathbf{k}-\mathbf{q})}.$$

$${\Sigma }_{x,y}={\Sigma }^{(1)}{k}_{x,y};{\Sigma }_z={\Sigma }^{(0)}+{\Sigma }^{(2)}{\mathbf{k}}^2.$$

$${\Sigma }^{(1)}\mathbf{k}$$

$${\Sigma }^{(0)}+{\Sigma }^{(2)}{\mathbf{k}}^2$$

$${\Sigma }^{(1)}\mathbf{k}$$

$$-\bar{m}-t_0+{\Sigma }^{(2)}{\mathbf{k}}^2$$

$$\delta m_{\mathrm{c}}^2=\frac{\left(\bar{m}+t_0\right)t_2^{\prime }}{{\xi }^2\ln \left\{1+{\left[t_2^{\prime }/(v^{\prime }\xi )\right]}^2\right\}},$$

$$v^{\prime }$$

$$t_2^{\prime }$$

$$h_{+}=-v(k_x{\sigma }_x+k_y{\sigma }_y)+[m(x)+t_2{\mathbf{k}}^2]{\sigma }_z.$$

$$iv{\sigma }_x{\partial }_x|\varphi (x)⟩+\left[m(x)-t_2{\partial }_x^2\right]{\sigma }_z|\varphi (x)⟩=0,$$

$$|\varphi (x)⟩=e^{{\gamma }_{0}x}|{\downarrow }_y⟩\{\begin{array}{ll}{e}^{-{\gamma }_{+}x},& x\ge 0\\ \cosh ({\gamma }_{-}x)-\frac{{\gamma }_{+}}{{\gamma }_{-}}\sinh ({\gamma }_{-}x),& x<0,\end{array}$$

$${\lambda }_{\pm }=\underset{x\to \pm \infty }{\lim }\left|x/\ln [|\varphi (x){|}^2]\right|=\frac{1}{2|\mathrm{Re}({\gamma }_0-{\gamma }_{\pm })|}.$$

$$\Delta x=\frac{t_2}{v}-\frac{v\bar{m}}{\delta m^2-{\bar{m}}^2}.$$

$$\delta m_{\mathrm{c}}=v\sqrt{\frac{\bar{m}}{t_2}}$$

Full text

Introduction**.** Thin films of magnetically doped topological insulators (TIs) provide an experimental realization of the quantum anomalous Hall (QAH) effect (Yu et al., 2010; Jiang et al., 2012; Wang et al., 2013a, b, 2014, 2015; Chang et al., 2013; Checkelsky et al., 2014; Kou et al., 2014; Bestwick et al., 2015; Chang et al., 2015; Kandala et al., 2015; Chang et al., 2016; Mani and Benjamin, 2018a, b), wherein a quantized Hall response emerges in the absence of an external magnetic field. For a ‘pure’ TI thin film (without magnetic dopants), the top and bottom surfaces host Dirac cones (Fu et al., 2007; Zhang et al., 2009; Hsieh et al., 2008; Xia et al., 2009; Chen et al., 2009; Moore and Balents, 2007) that can gap out via hybridization through the narrow bulk—naturally yielding a trivial insulator. When the Zeeman energy from polarized magnetic moments overwhelms the inter-surface hybridization, the TI film instead enters a QAH phase that exhibits a nontrivial Chern number together with a single chiral edge state that underlies conductance quantization. Studies of the magnetic structure (Lachman et al., 2015, 2017; Grauer et al., 2015) suggest that the magnetic dopants form weakly interacting, nanometer-scale islands and interact via easy-axis ferromagnetic coupling within each island. In typical experiments, these islands—which generically exhibit different coercive fields—are polarized by an external magnetic field, though ferromagnetic interactions allow the sample to remain magnetized even as the field is eliminated.

In this Letter we examine the de-magnetization process for the TI film, focusing on the regime in which the net magnetization vanishes. While one might expect that a trivial insulator supplants the QAH phase here, we show that a more interesting scenario can quite naturally emerge. In particular, a TI film with zero net magnetization experiences strong magnetic disorder and features locally polarized magnetic domains that cancel only on average. We show that this magnetic disorder can drive the system into a quantum-spin-Hall-like phase (first described in Ref. Hattori, 2015 in a different context) supporting helical edge modes that can be detected via standard transport measurements. Unlike the canonical quantum-spin-Hall state that is protected by time-reversal symmetry Kane and Mele (2005); Bernevig et al. (2006); König et al. (2007), these modes are protected by a unitary reflection symmetry that interchanges the top and bottom surfaces. The local magnetization in a thin TI film is not expected to vary appreciably along the perpendicular direction [see Fig. 1(a)]; the magnetically disordered film can then at least approximately preserve this reflection symmetry under appropriate gating conditions.

The physics we uncover can be viewed as a disorder-driven, zero-field counterpart to the very recently reported ‘quantum parity Hall effect’ in trilayer graphene, where edge channels protected by mirror symmetry arise Stepanov et al. (see also Ref. Barlas, 2018). We therefore refer to our helical phase as a ‘quantum anomalous parity Hall’ (QAPH) state. As an appealing application, we argue that the helical edge channels in a magnetically disordered TI film provide an ideal venue for pursuing Majorana zero modes (MZMs) Alicea (2010); Beenakker (2013); Aguado (2017); Lutchyn et al. (2018). In a usual quantum-spin-Hall state, Majorana zero modes bind to domain walls separating regions of the edge gapped by proximity-induced superconductivity and by time-reversal symmetry breaking Fu and Kane (2009). Our QAPH phase requires only breaking of reflection symmetry—thereby eschewing applied magnetic fields altogether and enabling dynamical manipulation of Majorana zero modes using purely electrical means.

Model**.** We consider a magnetically doped thin TI film (Zhang et al., 2010; Li et al., 2010; Shan et al., 2010; Hattori, 2015; Zhang et al., 2013) described at low energies by

[TABLE]

Here $\boldsymbol{r}=(x,y)$ is a coordinate along the film, $k_{x,y}=-i\partial_{x,y}$ is the momentum, $\sigma_{x,y,z}$ are Pauli matrices acting in spin space, and $\tau_{x,y,z}$ are Pauli matrices in the basis of states belonging to the top and bottom surfaces. The first term encodes the Dirac spectrum for each surface ($v$ is the velocity). The second hybridizes the two surfaces with tunneling matrix element $t(k)=t_{0}+t_{2}k^{2}$. In the last term, $m(\boldsymbol{r})$ is a Zeeman field induced by easy-axis magnetic dopants; note that the Zeeman field depends on $\boldsymbol{r}$ but is identical for the top and bottom surfaces. Upon disorder averaging we assume

[TABLE]

where $K(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ decays with correlation length $\xi$ and is normalized so that $K(0)=1$. With this normalization the disorder strength is set by $\delta m$.

Equation (1) commutes with $\mathcal{M}=\tau_{x}\sigma_{z}$, which implements a reflection about the $(x,y)$ plane. In the basis that diagonalizes $\mathcal{M}$, the Hamiltonian therefore acquires a block diagonal form, ${\mathcal{H}}=h_{+}\oplus h_{-}$, with

[TABLE]

Each block describes a single Dirac cone with a disordered mass term. As an illuminating primer, let us examine the clean limit where $m(\boldsymbol{r})=\bar{m}$. We assume a regularization of Eq. (3) such that the Chern number for block $h_{\pm}$ in this case is given by (Bernevig et al., 2006; Qi and Zhang, 2011)

[TABLE]

For $|\bar{m}|>|t_{0}|$, only one of the blocks has zero Chern number, and the overall Chern number is $C={\rm sgn}(\bar{m})$. This regime corresponds to the QAH phase that hosts a single chiral edge mode. For $|\bar{m}|<\left|t_{0}\right|$, the total Chern number necessarily vanishes. When $t_{0}t_{2}>0$ a trivial phase with $C_{\pm}=0$ arises. However, if $t_{0}t_{2}<0$ then the two blocks have nonzero and opposite Chern number: $C_{\pm}=\pm{\rm sgn}(t_{0})$. Here the system realizes a pristine QAPH phase supporting helical edge modes, with a right-mover coming from one block and a left-mover from the other. These edge modes are protected from gapping only when reflection symmetry is maintained. Henceforth we will assume $t_{0}t_{2}>0$—which precludes the QAPH phase in the clean limit. Below we show that introducing magnetic disorder through a spatially varying $m(\boldsymbol{r})$ nevertheless stabilizes the QAPH phase in a mechanism akin to that of the ‘topological Anderson insulator’ (Li et al., 2009; Jiang et al., 2009; Groth et al., 2009; Prodan, 2011; Yamakage et al., 2011; Xing et al., 2011).

Analysis of magnetic disorder**.** We now restore spatially non-uniform $m(\boldsymbol{r})$ in Eq. (1). The phase boundaries separating the QAH, trivial, and QAPH states highlighted above can be analytically estimated using the self-consistent Born approximation, wherein disorder effects are captured by a self-energy term $\Sigma_{\pm}(\omega,\boldsymbol{k})$ associated with block $h_{\pm}$. The self-energy follows from the self-consistent equation

[TABLE]

Here $\tilde{K}(\boldsymbol{q})$ is the Fourier transform of $K(\boldsymbol{r})$, and $\bar{h}_{\pm}$ is defined as $h_{\pm}$ evaluated with $m(\boldsymbol{r})\rightarrow\bar{m}$. Hereafter we set $\omega=0$, which allows us to extract the Chern numbers for the disordered system from an effective Hamiltonian $h_{\pm}^{{\rm eff}}(\boldsymbol{k})=\bar{h}_{\pm}(\boldsymbol{k})+\Sigma_{\pm}(\omega=0,\boldsymbol{k})$.

To facilitate analytical progress, we choose the function describing disorder correlations to be $K(\boldsymbol{r})=2(\xi/r)J_{1}(r/\xi)$, where $J_{1}$ is a Bessel function. The Fourier transform then takes a particularly simple form: $\tilde{K}(\boldsymbol{q})=4\pi\xi^{2}\Theta(1-\xi q)$ with $\Theta(x)$ the Heaviside step function. The low-momentum expansion of the self-energy takes the form SM

[TABLE]

which allows us to obtain corrections to $v$, $t_{2}$, and $\bar{m}\pm t_{0}$. We are after the critical disorder strength, $\delta m_{\rm c}^{\pm}$, at which $h^{\rm eff}_{\pm}$ changes Chern number. This transition occurs when $\Sigma_{\pm}^{(0)}=-(\bar{m}\pm t_{0})$. Using Eq. (6) and expanding the right-hand side of Eq. (5) to $\mathcal{O}(k^{2})$ yields SM

[TABLE]

where $v^{\prime}_{\pm},t_{2,\pm}^{\prime}$ are the disorder-renormalized $v$ and $t_{2}$ parameters SM . Note that $v^{\prime}_{\pm}\approx v$ and $t_{2,\pm}^{\prime}\approx t_{2}$ when either $(\bar{m}\pm t_{0})t_{2}/v^{2}$ or $v\xi/t_{2}$ are sufficiently small.

Dashed lines in Fig. 1(b) sketch the corresponding phase boundaries (see caption for parameters). Most interestingly, when $\delta{m}>\delta{m}_{+}^{{\rm c}},\delta{m}_{-}^{{\rm c}}$, the two blocks have nontrivial and opposite Chern numbers, and the system realizes the QAPH phase as advertised. If instead $\delta{m}<\delta{m}_{+}^{{\rm c}},\delta{m}_{-}^{{\rm c}}$ a trivial insulator emerges; otherwise the QAH state appears.

The physical picture underlying Eq. (7) can be understood from the limit of long disorder-correlation length, $\xi\gg v/\delta m,|t_{2}|/v$. Suppose first that $t_{0}=\bar{m}=0$ and $t_{2}=0^{+}$. In this limit, the $h_{+}$ block in Eq. (3) describes domains of characteristic size $\xi$ with either magnetization $\delta m$ (yielding trivial Chern number $C_{+}=0$) or $-\delta m$ (yielding $C_{+}=-1$); a chiral edge state propagates along each domain wall, reflecting the change in Chern number. Since the typical sizes for trivial and topological domains are equal here, the $h_{+}$ block overall is critical, in accordance with the percolation picture of the Chalker-Coddington network model (Chalker and Coddington, 1988).

We can then examine the effect of a finite $t_{2}>0$ on the position of the boundary mode between two domains, described by the Schrödinger equation $iv\sigma_{x}\partial_{x}|\phi(x)\rangle+\left[\delta m\cdot{\rm sgn}(x)-t_{2}\partial_{x}^{2}\right]\sigma_{z}|\phi(x)\rangle=0$ where $x=0$ is taken as the boundary (see the Supplemental Material for details (SM, )). The position of the edge mode with respect to the boundary can be quantified by the difference between the decay lengths towards either side of the boundary, $\Delta x=\lambda_{+}-\lambda_{-}$, where $\lambda_{\pm}=\lim_{x\to\pm\infty}\left|x/\ln[|\phi(x)|^{2}]\right|$. One obtains that to first order in $t_{2}\delta m/v^{2}$, the edge state shifts into the trivial domain by a distance $\Delta x=t_{2}/v$ [see Fig. 2(a)], thereby enlarging the topological region and pushing the block into the topological phase. Alternatively, introducing small but finite average magnetization and inter-surface tunneling, $\bar{m},t_{0}\neq 0$, instead shifts the edge state by $\Delta x=-v(t_{0}+\bar{m})/\delta{m}^{2}$. The phase transition therefore occurs when these shifts cancel, corresponding to $\delta{m}^{c}_{+}=v[(t_{0}+\bar{m})/t_{2}]^{1/2}$, which indeed agrees with Eq. (7) in the limit of $\xi\gg t_{2}/v$ and $(\bar{m}+t_{0})t_{2}/v^{2}\ll 1$. Similarly, for the $h_{-}$ block [Fig. 2(b)] one finds a transition at $\delta m_{-}^{c}=v[(t_{0}-\bar{m})/t_{2}]^{1/2}$, again in agreement with Eq. (7).

The above analytic results can be corroborated numerically. To this end we discretize Eq. (1) on an $L_{x}\times L_{y}$ square lattice, resulting in

[TABLE]

where $c_{\boldsymbol{n},\rho,s}^{\dagger}$ creates an electron with spin $s$ on site $\boldsymbol{n}$ of surface $\rho$ (we suppress $\rho,s$ indices above); lengths are measured in units of the lattice constant; and $t_{0}^{\prime}=t_{0}+4t_{2}$. The Zeeman field $m_{\boldsymbol{n}}$ is now taken to be normally distributed, with average $\left\langle m_{\boldsymbol{n}}\right\rangle=\bar{m}$ and correlations $\left\langle m_{\boldsymbol{n}}m_{\boldsymbol{n}^{\prime}}\right\rangle-\bar{m}^{2}=\delta m^{2}\exp\left\{-\left[\left(\boldsymbol{n}-\boldsymbol{n}^{\prime}\right)/2\xi\right]^{2}\right\}$. We connect the system to two leads, as depicted in Fig. 1(a), and compute the scattering matrix for incoming and outgoing electrons (SM, ) using a recursive procedure that involves gradually increasing the system’s length in the $x$ direction Lee and Fisher (1981). The two-terminal conductance obtained from the Landauer-Büttikker formalism is $G=\left(e^{2}/h\right){\rm Tr}(t^{\dagger}t),$ where $t$ is the transmission matrix between the leads for electrons at the Fermi energy, $E_{\rm F}$.

The color map in Fig. 1(b) displays the conductance $G$ versus $\bar{m}$ and $\delta m$. Each edge mode contributes $e^{2}/h$ to the conductance. The trivial, QAH, and QAPH phases are thus readily diagnosed by quantized conductances $G=0$, $e^{2}/h$, and $2e^{2}/h$, respectively. In the clean limit ($\delta m=0$) one obtains the familiar scenario where the system passes from the trivial to the QAH phase phase when $\left|\bar{m}\right|>t_{0}$. Magnetic disorder instead drives the system into the QAPH state, as found analytically; note the good agreement between the analytical and numerical phase boundaries, despite the different disorder correlations used.

Reflection symmetry breaking**.** To study the effects of breaking the reflection symmetry $\mathcal{M}=\tau_{x}\sigma_{z}$ that protects the QAPH edge states, we include an electric potential near the sample boundary that is opposite for the top and bottom surfaces; experimentally such a term can be controllably generated via asymmetric gating of the TI film. We specifically perturb Eq. (1) [or its lattice counterpart, Eq. (8)] with $\mathcal{H}^{\prime}=V_{{\rm A}}(\textbf{r})\tau_{z}$, where $V_{\rm A}(\boldsymbol{r})=V_{\rm A}^{0}$ within a distance $W_{\rm edge}$ from the edges and $V_{\rm A}(\boldsymbol{r})=0$ otherwise. Figure 3(a) presents the two-terminal conductance versus $V_{{\rm A}}$ for different linear system sizes $L_{x}=L_{y}\equiv L$, assuming system parameters corresponding to the QAPH phase. For $V_{{\rm A}}^{0}=0$, the conductance reaches the quantized value of $G=2e^{2}/h$ as expected, independent of system size. Increasing $V_{{\rm A}}^{0}$ generates backscattering among the helical edge states and thus reduces the conductance. The system-size dependence is further explored in Fig. 3(b), which plots the conductance versus $L$ for three values of $V_{{\rm A}}^{0}$. For a given $V_{{\rm A}}^{0}\neq 0$, the probability of edge-mode backscattering increases with system size, thereby decreasing the conductance—albeit rather slowly. Interestingly, since the counter-propagating modes are not related by symmetry, they generally do not overlap in space (see also (SM, )). This property can suppress backscattering by reflection-symmetry-breaking terms such as $V_{{\rm A}}^{0}$.

Superconductivity and Majorana modes**.** The helical edge states in the QAPH phase serve as a natural platform for realizing MZMs upon coupling the edge to a conventional superconductor (Fu and Kane, 2009). In this setup a MZM localizes to the boundary between a section of the edge gapped by superconductivity and a section gapped due to a reflection-symmetry-breaking potential $V_{{\rm A}}$ (see above). The latter gapped regions can be accessed by making $V_{\rm A}$ arbitrarily large without deleteriously impacting the parent superconductor—contrary to applied magnetic fields which alternative approaches typically require Fu and Kane (2009); Oreg et al. (2010); Lutchyn et al. (2010); Pientka et al. (2013). Furthermore, locally controlling $V_{\rm A}$ through gates enables all-electrical manipulation of MZMs.

To demonstrate the realization of MZMs, we simulate superconductivity in the setup from Fig. 4(a) by adding a pairing term, $\Delta_{\boldsymbol{n}}c_{\boldsymbol{n}\uparrow}^{\dagger}c_{\boldsymbol{n}\downarrow}^{\dagger}+{\rm h.\,c.}$, to Eq. (8). Proximitizing the superconductor generally also induces an asymmetry in the chemical potential, which is simulated by a term $V_{{\rm A},\boldsymbol{n}}\boldsymbol{c_{n}}^{\dagger}\tau_{z}\boldsymbol{c_{n}}$. The potentials $\Delta_{\boldsymbol{n}}$ and $V_{{\rm A},\boldsymbol{n}}$ assume the values $\Delta_{0}$ and $V_{\rm A}^{0}$ beneath the superconductor but otherwise vanish. We then recalculate the scattering matrix, which now includes a block $r_{\rm he}$ describing Andreev reflection (SM, ). Figure 4(b) presents the total Andreev reflection, $R_{{\rm he}}={\rm Tr}[r_{{\rm he}}^{\dagger}r_{{\rm he}}^{\phantom{\dagger}}]$, for incident electrons at zero energy, for different values of $V_{{\rm A},\boldsymbol{n}}$ (see caption for parameters). For $\left|\bar{m}\right|\lesssim 0.8$, the system forms a QAPH phase whose helical edge states are gapped by superconductivity, and accordingly Andreev reflection occurs with near-unit probability. The perfect Andreev reflection at zero energy signals the emergence of a MZM at the boundary of the superconducting section of the edge (Law et al., 2009; Fidkowski et al., 2012). For $\bar{m}\gtrsim 0.8$, the QAH phase appears, accompanied by a precipitous suppression of $R_{\rm he}$ due to the inability of a chiral mode to be reflected (either through normal or Andreev processes).

Discussion**.** We have shown that experimentally motivated disorder originating from randomly oriented magnetic islands in TI thin films can stabilize a QAPH state. This phase harbors helical edge states that are protected by a reflection symmetry that interchanges the top and bottom surfaces, and can be straightforwardly detected: In a two-terminal measurement, the QAPH state is characterized by a quantized $2e^{2}/h$ conductance that is enhanced compared to that of the proximate QAH phases; recall Fig. 1(b). In a Hall-bar measurement, the QAPH should appear as a $\sigma_{xy}=0$ plateau Feng et al. (2015), together with $\sigma_{xx}=2e^{2}/h$.

Breaking reflection symmetry suppresses the conductance below $2e^{2}/h$, though edge conduction should still be observable in a finite system (see Fig. 3). Tuning in and out of the reflection-symmetric regime, while keeping the chemical potential fixed, can be achieved by employing both bottom and top gates. We argued that the ability to electrically control the helical edge modes in this manner renders the proximitized QAPH system an ideal venue for exploring MZMs. Furthermore, having local control over the breaking of reflection symmetry (e.g. using gates) can allow for binding fractional charges at domain walls between regions where the reflection-symmetry-breaking term switches sign, analogous to the bound states discussed in Ref. Qi et al., 2008.

We close by providing a complementary perspective on our findings. In a clean integer quantum Hall system, populating a spin-degenerate Landau level changes the Hall conductivity $\sigma_{xy}$ from 0 to $2e^{2}/h$. Including disorder that breaks spin conservation generically splits this plateau transition and opens an intervening quantum Hall phase with $\sigma_{xy}=e^{2}/h$ Lee and Chalker (1994); see Fig. 5(a). In a clean magnetic TI film with no inter-surface tunneling, reversing the magnetization similarly changes $\sigma_{xy}$ from $-e^{2}/h$ to $+e^{2}/h$. By analogy with the plateau transition, it is natural to expect that disorder can open up an intervening phase with $\sigma_{xy}=0$ Wang et al. (2014) as sketched in Fig. 5(b). Interestingly, we have shown that this disorder-induced phase need not be trivial—even though the Hall conductivity vanishes—when reflection symmetry is present. This viewpoint may prove useful for discovering other disorder-induced symmetry-protected phases, e.g., in bands with higher Chern numbers.

Acknowledgments**.** We thank A. Thomson for her insightful comments. We acknowledge support from the Army Research Office under Grant Award W911NF-17-1-0323 (J.A.); the NSF through grants DMR-1723367 (J.A.); the Caltech Institute for Quantum Information and Matter, an NSF Physics Frontiers Center with support of the Gordon and Betty Moore Foundation through Grant GBMF1250 (J.A.); and the Walter Burke Institute for Theoretical Physics at Caltech (J.A. and A.H.).

Supplemental Material

I Numerical simulation

We begin by rewriting the lattice Hamiltonian, Eq. (8) of the main text, in the following form

[TABLE]

where $\{h_{n_{x}}\}_{n_{x}=1}^{L_{x}}$ and $V$ are $4L_{y}\times 4L_{y}$ matrices, $\vec{\psi}_{n_{x}}^{\dagger}=(c_{(n_{x},1),{\rm b},\uparrow}^{\dagger},c_{(n_{x},1),{\rm b},\downarrow}^{\dagger},c_{(n_{x},1),{\rm t},\uparrow}^{\dagger},c_{(n_{x},1),{\rm t},\downarrow}^{\dagger},\dots,c_{(n_{x},L_{y}),{\rm b},\uparrow}^{\dagger},c_{(n_{x},L_{y}),{\rm b},\downarrow}^{\dagger}$, $c_{(n_{x},L_{y}),{\rm t},\uparrow}^{\dagger},c_{(n_{x},L_{y}),{\rm t},\downarrow}^{\dagger})$ is a $1\times 4L_{y}$ vector of creation and annihilation operators, and ‘b’ and ‘t’ indicate bottom and top surfaces, respectively.

We place two normal-metal leads at $x=0$ and $x=L_{x}$. The reflection matrix for electrons incident from the right can be calculated from the Green function, $G_{n_{x}}(\omega)$, using  Fisher and Lee (1981); Iida et al. (1990)

[TABLE]

where $\rho_{{\rm R}}$ is the density of states in the right lead, and $G_{L_{x}}$ is the Green function matrix at the right-most sites of the system, obtained through the recursive relation Lee and Fisher (1981)

[TABLE]

For every $n_{x}=1,\dots,L_{x}$, the Green function $G_{n_{x}}(\omega)$ is a $4N_{y}\times 4N_{y}$ matrix (indices running over ‘b’/‘t’, spin, and $n_{y}$), and $G_{0}=i\pi\rho_{{\rm L}}$, where $\rho_{{\rm L}}$ is the density of states in the left lead. The linear conductance is then given by

[TABLE]

where $t$ is the transmission matrix evaluated at the Fermi energy $\omega=E_{{\rm F}}$.

When simulating the proximity-coupled superconductor, we amend the Hamiltonian via $H\to H+\sum_{\boldsymbol{n},\rho}(\Delta_{\boldsymbol{n}}c_{\boldsymbol{n}\rho\uparrow}^{\dagger}c_{\boldsymbol{n}\rho\downarrow}^{\dagger}+{\rm H.c.})+\sum_{\boldsymbol{n},s}V_{{\rm A},\boldsymbol{n}}(c_{\boldsymbol{n},{\rm b},s}^{\dagger}c_{\boldsymbol{n},{\rm b},s}-c_{\boldsymbol{n},{\rm t},s}^{\dagger}c_{\boldsymbol{n},{\rm t},s})$ and rewrite it in a Bogoliubov-de Gennes form,

[TABLE]

where $\vec{\psi}_{n_{x}}^{\dagger}=(\vec{\psi}_{n_{x}}^{\dagger},\vec{\psi}_{n_{x}})$, and $\{\tilde{h}_{n_{x}}$}, $\tilde{V}$ are $8N_{y}\times 8N_{y}$ matrices. The reflection matrix, $\tilde{r}(\omega)$ is given by Eqs. (10) and (11), with $\{h_{n_{x}}\}$ and $V$ replaced by $\{\tilde{h}_{n_{x}}\}$ and $\tilde{V}$ , respectively. The diagonal blocks $r_{{\rm ee}}$ and $r_{{\rm hh}}$ describe normal reflection, while the off-diagonal ones, $r_{{\rm he}}$ and $r_{{\rm eh}},$ describe Andreev reflection.

I.1 Local density of states

It was pointed out in the main text that the counter-propagating modes in the QAPH phase are not related by symmetry, and therefore generally do not overlap in space. This can be seen most conveniently in the local density of states. In Fig. 6 we present the local density of states, $\rho_{\pm}(x,y,\omega=0)$, corresponding to modes even and odd under the reflection symmetry, calculated according to the method described in Ref. (Potter and Lee, 2011). The results are shown for a single disorder realization having $\bar{m}=0$ and $\delta m=2$. The rest of the system parameters are the same as in Fig. 1(b) of the main text. For these parameters, the system is in the QAPH phase. Both the $+$ and $-$ sectors host a mode confined roughly to the edge, though the positions clearly do not coincide.

II Self-consistent Born approximation

We start from Eq. (5) of the main text and focus on the $h_{+}$ block (the treatment of the $h_{-}$ block follows along the same lines). The self-energy can generally be written as

[TABLE]

Substituting the above expansion in Eq. (5) of the main text results in four coupled equations for $\Sigma_{0,x,y,z}(k,\omega)$. We are interested in the self-energy at $\omega=0$ from which we define an effective Hamiltonian,

[TABLE]

For $\omega=0$, the self-consistent equations read

[TABLE]

We have omitted above the $\omega$ argument for brevity.

First, we notice that Eq. (16a) is satisfied by $\Sigma_{0}(\boldsymbol{k})=0$. Second, in accordance with our long-wavelength analysis of the system, we expand $\Sigma_{x,y,z}(\boldsymbol{k})$ to second order in $\boldsymbol{k}$,

[TABLE]

The form of the above expansion derives from the symmetry properties of the self-consistent equations, Eq. (16), whereby $(\Sigma_{x},\Sigma_{y})$ transform as $(k_{x},k_{y})$, and $\Sigma_{z}$ transforms as a scalar under rotations and reflections of $\boldsymbol{k}$. Plugging this expansion in Eq. (16), we have

[TABLE]

where $v^{\prime}\equiv v-\Sigma^{(1)}$, and $t_{2}^{\prime}=t_{2}+\Sigma^{(2)}$.

We are looking for the critical disorder strength, $\delta m_{\rm c}$, for which the $h_{+}$ block is at a transition between a trivial and a topological phase. This happens when the $\boldsymbol{k}=0$ coefficient of the $\sigma_{z}$ term in Eq. (15) vanishes, namely when $\Sigma_{z}(\boldsymbol{k}=0)=-\bar{m}-t_{0}$. We substitute this condition in Eq. (18), and solve the integrals to second order in $\boldsymbol{k}$, which results in the following equations

[TABLE]

Finally, requiring that Eqs. (19a,19b) are obeyed order by order results in the expression for the critical disorder strength,

[TABLE]

and the following equations for $v^{\prime},t_{2}^{\prime}$,

[TABLE]

Note that Eqs. (21a,21b) are solved by $v^{\prime}\simeq v$, $t_{2}^{\prime}\simeq t_{2}$ when either $(\bar{m}+t_{0})t_{2}/v^{2}$ or $v\xi/t_{2}$ are sufficiently small.

III The long-range-disorder limit

As explained in the main text, in the large $\xi$ limit the results can be understood by considering the effect of $t_{2}$ on the position of the boundary mode between two domains. Consider, for example, the sector $h_{+}$ and let us start from the case $t_{0}=0$. We examine the boundary ($x=0$) between two domains of opposite-sign mass,

[TABLE]

where $\bar{m}(x)=\bar{m}+\delta m\cdot{\rm sgn}(x),$ and $\delta m>|\bar{m}|$. For simplicity of presentation we assume $t_{2}>0$, which means that $x<0$ is topological, while $x>0$ is trivial.

We look for a mid gap solution to the Schrödinger equation,

[TABLE]

where we have set $k_{y}=0$. The solution is given by

[TABLE]

where $\gamma_{0}=\frac{v}{2t_{2}}$ and $\gamma_{\pm}=\frac{v}{2t_{2}}\sqrt{1+4t_{2}(\bar{m}\pm\delta m)/v^{2}}$.

The wave function decays away from the boundary at $x=0$, however, the decay lengths to the left and to the right are not equal. They are given by

[TABLE]

When $\lambda_{-}>\lambda_{+}$, the wave function is shifted towards the left, effectively decreasing the size of the topological region and overall driving $h_{+}$ to the trivial phase. When $\lambda_{+}>\lambda_{-}$, on the other hand, the wave function is shifted towards the right thereby increasing the size of the topological region and driving $h_{+}$ to the topological phase. This shift can be quantified by $\Delta x=\lambda_{+}-\lambda_{-}$, which in the limit of $\delta mt_{2}/v^{2}\ll 1$ is given by

[TABLE]

The transition between the trivial and topological phases of $h_{+}$ occurs when $\Delta x=0$, which to leading order in $\bar{m}t_{2}/v^{2}$ (recall $|\bar{m}|<\delta m$) is at the critical disorder strength

[TABLE]

Finally, the introduction of a finite $t_{0}$ is accounted for by taking $\bar{m}\to\bar{m}+t_{0}$, while for the odd symmetry sector $h_{-}$ one takes $\bar{m}\to\bar{m}-t_{0}$ and $t_{2}\to-t_{2}$. This reproduces Eq. (7) of the main text in the limit of $\xi\gg t_{2}/v$ and $(t_{0}\pm\bar{m})t_{2}/v^{2}\ll 1$.