Time-reversal-symmetric topological magnetoelectric effect in three-dimensional topological insulators

TL;DR

This paper explores a novel time-reversal-symmetric topological magnetoelectric effect in three-dimensional topological insulators, where the electric charge depends solely on the magnetic field's magnitude, not its direction.

Contribution

It introduces a nonperturbative theoretical framework for a time-reversal-symmetric TME and proposes an experimental setup to observe this effect.

Findings

Electric charge depends only on magnetic field magnitude

The effect is independent of magnetic field sign

Proposed mesoscopic experimental demonstration

Abstract

One of the hallmarks of time-reversal-symmetric topological insulators in three dimensions is the topological magnetoelectric effect (TME). So far, a time-reversal breaking variant of this effect has attracted much attention, in the sense that the induced electric charge changes sign when the direction of an externally applied magnetic field is reversed. Theoretically, this effect is described by the so-called axion term. Here, we discuss a time-reversal-symmetric TME, where the electric charge depends only on the magnitude of the magnetic field but is independent of its sign. We obtain this nonperturbative result by a combination of analytic and numerical arguments, and suggest a mesoscopic setup to demonstrate it experimentally.

Figures

Equations

$$\rho ({\mathbf{E}}_z,\Phi )=\frac{e}{2}\frac{w(\Phi )}{\sqrt{V(\Phi {)}^2+\Delta (\Phi {)}^2}}V(\Phi ),$$

$$⟨\hat{\rho }(\mathbf{x}){⟩}_{A_0,\mathbf{A}}=(-e)\sum _{E_{\alpha }\le 0}|{\eta }_{\alpha }(\mathbf{x}){|}^2,$$

$$\delta ⟨\hat{\rho }(\mathbf{x}){⟩}_{A_0,\mathbf{A}}=(-e)\left[\frac{1}{2}\sum _{E_{\alpha }\le 0}|{\eta }_{\alpha }(\mathbf{x}){|}^2-\frac{1}{2}\sum _{E_{\beta }>0}|{\eta }_{\beta }(\mathbf{x}){|}^2\right].$$

$$\hat{H}=\left(\begin{array}{cc}0& h^{+}\\ h^{-}& 0\end{array}\right),h^{\pm }=\mp \left({\partial }_{\theta }+\frac{1}{2}\cot \theta \right)+\frac{i{\partial }_{\varphi }}{\sin \theta }+eR{\mathbf{A}}_{\varphi },$$

$$h^{+}{h}^{-}$$

$$+\frac{1}{1-x^2}{\left[-m+\frac{1}{2}x+N_{\varphi }\right]}^2+\frac{d{N}_{\varphi }}{dx}+\frac{1}{2},$$

$${\psi }_{\uparrow }(x)=\{\begin{array}{l}(1-x{)}^{\frac{1}{2}(m-\frac{1}{2})}(1+x{)}^{\frac{1}{2}(m+\frac{1}{2})-\frac{N_0}{\delta }}{g}_0(x)\\ (1-x{)}^{\frac{1}{2}(m-N_0-\frac{1}{2})}(1+x{)}^{\frac{1}{2}(m-N_0+\frac{1}{2})}g(x)\\ (1-x{)}^{\frac{1}{2}(m-\frac{1}{2})-\frac{N_0}{\delta }}(1+x{)}^{\frac{1}{2}(m+\frac{1}{2})}{g}_1(x),\end{array}$$

$$g$$

$$+\beta {\xi }^{1-c}F\left(1-E,1+E,2-c,\xi \right)$$

$$g^{\prime }(\epsilon )/g(\epsilon )-g_0^{\prime }(\epsilon )/g_0(\epsilon )=N_0/\epsilon$$

$$E=\pm \{\begin{array}{ll}c+n& \text{ if }c>0\\ 0\text{ or }1-c+n& \text{ if }c<0,\end{array}$$

$${\psi }_{\uparrow }(x)=\frac{1}{\sqrt{\mathcal{N}}}(1-x{)}^{\frac{1}{2}(-\Phi /{\Phi }_0)}(1+x{)}^{\frac{1}{2}(1-\Phi /{\Phi }_0)}$$

$$\hat{H}=-2\lambda {\tau }_z\sum _{\mu =1}^3{\sigma }_{\mu }\sin (k_{\mu })+{\tau }_x\left(\kappa -2t\sum _{\mu =1}^3\cos (k_{\mu })\right),$$

$$⟨\hat{\rho }(\mathbf{x}){⟩}_{A_0,\mathbf{A}}=(-e)\sum _{E_{\alpha }\le 0}|{\eta }_{\alpha }(\mathbf{x}){|}^2,$$

$$⟨\hat{\rho }(\mathbf{x}){⟩}_{A_0,\mathbf{A}}=⟨\hat{\rho }(\mathbf{x}){⟩}_{A_0,-\mathbf{A}}.$$

$$h^{+}{h}^{-}=$$

$$\alpha =\eta \frac{1}{2}\left(m-\frac{1}{2}-A-B\right),\beta =\eta \frac{1}{2}\left(m+\frac{1}{2}+A-B\right),\eta =\pm 1$$

$$\displaystyle\Bigg{[}\xi(1-\xi)\frac{d^{2}}{d\xi^{2}}+\eta\Bigg{(}m-B-A+\eta-\frac{1}{2}$$

$$\xi (1-\xi )\frac{d^{2}g}{d{\xi }^2}+[c-(a+b+1)\xi ]\frac{dg}{d\xi }-abg=0$$

$$d:=m-B+1/2,c=d-A,a+b=2d,ab=d^2-A^2-E^2.$$

$${\psi }_{\uparrow }(x)=\{\begin{array}{ll}(1-x{)}^{\frac{1}{2}(m-\frac{1}{2})}(1+x{)}^{\frac{1}{2}(m+\frac{1}{2})-\frac{N_0}{\delta }}{g}_0(x)& \text{ if }1-\delta \le x\le 1\\ (1-x{)}^{\frac{1}{2}(m-N_0-\frac{1}{2})}(1+x{)}^{\frac{1}{2}(m-N_0+\frac{1}{2})}g(x)& \text{ if }-1+\delta \le x\le 1-\delta \\ (1-x{)}^{\frac{1}{2}(m-\frac{1}{2})-\frac{N_0}{\delta }}(1+x{)}^{\frac{1}{2}(m+\frac{1}{2})}{g}_1(x)& \text{ if }-1\le x\le -1+\delta .\end{array}$$

$$\frac{g^{\prime }(\epsilon )}{g(\epsilon )}-\frac{g_0^{\prime }(\epsilon )}{g_0(\epsilon )}=\frac{N_0}{\epsilon },\text{ and }\frac{g_1^{\prime }(1-\epsilon )}{g_1(1-\epsilon )}-\frac{g^{\prime }(1-\epsilon )}{g(1-\epsilon )}=\frac{N_0}{\epsilon }.$$

$$g_0(\xi )$$

$$g(\xi )$$

$$=\alpha F\left(c-E,c+E,c,\xi \right)+\beta {\xi }^{1-c}F\left(1-E,1+E,2-c,\xi \right)$$

$$g_1(\xi )$$

$$g(\xi )$$

$$\tilde{\alpha }=\alpha \frac{\Gamma (c)\Gamma (-c)}{\Gamma (E)\Gamma (-E)}+\beta \frac{\Gamma (2-c)\Gamma (-c)}{\Gamma (1-c+E)\Gamma (1-c-E)},\tilde{\beta }=\alpha \frac{\Gamma (c)\Gamma (c)}{\Gamma (c-E)\Gamma (c+E)}+\beta \frac{\Gamma (2-c)\Gamma (c)}{\Gamma (1+E)\Gamma (1-E)}.$$

$$0$$

$$0=\alpha \left[-r_0{E}^2\right]+\beta (1-c){\epsilon }^{-c},0=\tilde{\alpha }+\tilde{\beta }\left[-\frac{c}{p_1-N_0}+1\right]{\epsilon }^{-c}$$

$$0=\frac{\Gamma (c)\Gamma (c)}{\Gamma (c-E)\Gamma (c+E)}\text{textperiodcentered}\left[-\frac{c}{p_1-N_0}+1\right].$$

$$E=\pm (c+n),\text{ with }n=0,1,2,\dots \text{ for }c=m+1/2-N_0>0.$$

Full text

Introduction. Time-reversal-symmetric (TRS) topological insulators (TIs) Hasan and Kane (2010); Qi and Zhang (2011) are a fascinating class of electronic materials with insulating bulk and topologically protected surface states, which are either gapless, break a symmetry, or feature topological order Vishwanath and Senthil (2013). Evidence for the existence of such surface states comes from spin textures observed in photoemission experiments Xia et al. (2009); Zhang et al. (2009) and from the observation of a half-integer quantum Hall effect Brüne et al. (2011); Chang et al. (2013); Xu et al. (2014); Yoshimi et al. (2015).

From a theoretical point of view, the hallmark response of TRS TIs in three dimensions (3D) is the topological magnetoelectric effect Qi et al. (2008); Essin et al. (2010). So far, a time-reversal (TR) breaking variant of this effect has attracted much attention. When TR is broken by, say, a magnetic coating with Zeeman coupling to the TI surface, a Hall conductivity $\sigma_{xy}=\frac{\tilde{\theta}}{2\pi}\frac{e^{2}}{2\pi}$ arises, where $\tilde{\theta}$ is quantized to $\tilde{\theta}=\pm\pi$. Then, the insertion of a magnetic flux tube gives rise to the accumulation of a charge $|Q|=e/2$ per flux quantum $\Phi_{0}=h/e$. Importantly, $\mathop{\text{sgn}}(Q)$ depends on the direction of the magnetic field inside the flux tube, i.e., the response is not TRS. A consequence of the surface Hall conductivity is quantized Kerr and Faraday rotations Tse and MacDonald (2010); Maciejko et al. (2010), which have recently been confirmed experimentally Wu et al. (2016); Okada et al. (2016); Dziom et al. (2017).

In the presence of TRS, the linear magnetoelectric response vanishes Landau and Lifschitz (1984); Note (17). Thus, strictly speaking, all variants of the topological magnetoelectric effect are nonlinear effects, as they require an additional perturbation, say, a Zeeman coupling on the surface as above. The absence of a linear magnetoelectric response may seem to be at odds with the fact that the bulk of a 3D TRS TI has been characterized by the so-called axion action Qi et al. (2008); Essin et al. (2010) $S_{{\theta}}=\frac{{\theta}}{2\pi}\frac{e^{2}}{2\pi}\int d^{3}x\,dt\,{\mathbf{E}}\cdot{\mathbf{B}}$. Under a TR transformation, ${\mathbf{E}}\to{\mathbf{E}}$, ${\mathbf{B}}\to-{\mathbf{B}}$, and $S_{{\theta}}\to-S_{{\theta}}$. Classically, this action breaks TRS, but quantum mechanically, only the Feynman amplitude $e^{iS_{{\theta}}}$ needs to be symmetric. If the electronic wave functions and electromagnetic fields satisfy periodic boundary conditions, one can show that this integral is quantized to integer multiples of $4\pi^{2}/e^{2}$ Milnor and Stasheff (1974); Vazifeh and Franz (2010). This implies that $S_{{\theta}}={\theta}$ modulo $2\pi$, hence $\theta=\pm\pi$ would respect TRS also. Now, for a TI with boundaries, by using partial integration, $S_{{\pm\pi}}$ can be converted into the surface quantum Hall term discussed above, and it seems naively that the axion action indeed describes a magnetoelectric effect. However, this is not the case, as the surface response due to $S_{\pm\pi}$ is canceled by additional contributions from the surface states, in this case the parity anomaly, to restore TRS Mulligan and Burnell (2013); Zirnstein and Rosenow (2013); Witten (2016).

In this work, we describe a TR-symmetric topological magnetoelectric effect. This means that the accumulated electric charge depends only on the magnitude of the magnetic flux but is independent of its sign (magnetic field direction). In particular, we consider a spherical TI threaded by a thin magnetic flux tube, and subject to a small uniform electric field [see Fig. 1(a)]. If $Q_{\text{top}}(\mathbf{E}_{z},\Phi)$ denotes the total charge on the top half, we show that the additional charge due to the insertion of one flux quantum is $\Delta Q(\mathbf{E}_{z},\pm\Phi_{0})\equiv Q_{\text{top}}(\mathbf{E}_{z},\pm\Phi_{0})-Q_{\text{top}}(\mathbf{E}_{z},0)=+(e/2)\mathop{\text{sgn}}\mathbf{E}_{z}$. This response is TRS, in contrast to the TR-breaking surface quantum Hall effect Brüne et al. (2011); Chang et al. (2013); Xu et al. (2014); Yoshimi et al. (2015); Nomura and Nagaosa (2011); König et al. (2014); Lee (2009); Vafek (2011); Nogueira et al. (2016) described above. In addition, we numerically analyze a lattice model, and suggest a mesoscopic setup to demonstrate this effect experimentally. These results are not related to the wormhole effect Rosenberg et al. (2010), which only occurs when the diameter of the flux tube is much smaller than the lattice spacing. Here, we instead consider a flux tube that covers many plaquettes, so that the surface states of the TI can be described by a continuum two-dimensional (2D) Dirac Hamiltonian. Then, on an infinite planar surface, inserting one flux quantum $\Phi_{0}$ will give rise to a single, spin-polarized zero-energy state localized (power law) at the tube Aharonov and Casher (1979); Jackiw and Schrieffer (1981); Polychronakos (1987). Our treatment of a finite geometry goes significantly beyond this result.

Physical picture. The fractional charge $e/2$ arises because an occupied delocalized eigenstate is transformed into an occupied localized eigenstate when adiabatically inserting one flux quantum: In a spherical geometry, the delocalized state has equal weight in both hemispheres, whereas the localized state is contained in one of them; this corresponds to a change of charge in that hemisphere by $e/2$. The presence of a small background field $\mathbf{E}_{z}$ is important for this to happen: Flux insertion generates two localized states, one at each pole, which would otherwise hybridize with each other, resulting in delocalized eigenstates again, and no charge difference would be observable. Only when the energy difference between the localized states, due to $\mathbf{E}_{z}$, is larger than the hybridization, one localized state will be occupied, and the other one empty.

In order to evaluate the induced charge, we compute the two eigenstates with energies closest to the Fermi level for each value of the flux $\Phi$ between zero and $\Phi_{0}$, and then evaluate the effect of an external electric field $\mathbf{E}_{z}$ in the subspace spanned by these states. This projection to a low-energy subspace becomes exact if one first takes the limit of an infinitesimally thin flux tube [see Fig. 2(a)], and then considers an infinitesimal electric field $\mathbf{E}_{z}$. In this way, the infinitesimal $\mathbf{E}_{z}$ does not polarize the system for zero flux, but due to the order of limits described above, at $\Phi=\Phi_{0}$ the level splitting caused by $\mathbf{E}_{z}$ is much larger that the hybridization energy. In a numerical calculation, both the flux tube diameter and the electric field are small but finite, and the exact analytic result is recovered by finite-size scaling [see Fig. 3(b)].

The two lowest-energy eigenstates are superpositions of spin-polarized wave functions $\eta_{\text{n}}(\Phi,{\mathbf{x}})$ and $\eta_{\text{s}}(\Phi,{\mathbf{x}})$ located at the north and south pole; explicit expressions are derived later. Projecting to these states, the Hamiltonian becomes ${\hat{H}}=\Delta s_{x}+Vs_{z}$ where $s_{x}$, $s_{z}$ are Pauli matrices, $\Delta$ is the hybridization energy, and $V$ is the projected potential energy $V(\Phi)=ed(\Phi)\mathbf{E}_{z}$ with $d(\Phi)=\langle\eta_{\text{n}}|{\mathbf{x}}_{3}|\eta_{\text{n}}\rangle$ the dipole moment. The charge response is $\Delta Q(\mathbf{E}_{z},\Phi)=\rho(\mathbf{E}_{z},\Phi)-\rho(\mathbf{E}_{z},0)$ with

[TABLE]

and $w(\Phi)=\int_{\text{top}}d{\mathbf{x}}\,[\absolutevalue{\eta_{\text{n}}(\Phi,{\mathbf{x}})}^{2}-\absolutevalue{\eta_{\text{s}}(\Phi,{\mathbf{x}})}^{2}]$. In the limit of a thin flux tube, the hybridization equals the energy $\Delta(\Phi)=(\Phi/\Phi_{0}-1)v_{F}/R$ for $\Phi<\Phi_{0}$ where $v_{F}$ is the Fermi velocity and $R$ the radius of the sphere. The corresponding charge response is depicted in Fig. 1(b). Since we have $\Delta(\Phi)\to 0$ and $w(\Phi)\to 1$ for $\Phi\to\Phi_{0}$, an infinitesimally small electric field is sufficient to lift the degeneracy, and we find that $\Delta Q(\mathbf{E}_{z},\Phi_{0})=+(e/2)\mathop{\text{sgn}}\mathbf{E}_{z}$ for $\mathbf{E}_{z}\to 0$.

Fractional charge. The appearance of a half-integer charge due to localized zero-energy states can be interpreted as an instance of charge fractionalization Jackiw and Schrieffer (1981); Krive and Rozhavskiĭ (1987); Leinaas (2009); Leinaas et al. (2009). To make the connection to condensed matter, we consider a general lattice Hamiltonian ${\hat{H}}(A_{0},{\mathbf{A}})$ which depends on an electric potential $A_{0}$ and a vector potential ${\mathbf{A}}$, and has $n$ internal degrees of freedom per site ${\mathbf{x}}$. At zero temperature and chemical potential, the expectation value of the charge density is

[TABLE]

where $E_{\alpha}$ are the energy eigenvalues of ${\hat{H}}(A_{0},{\mathbf{A}})$, and $\eta_{\alpha}({\mathbf{x}})$ the corresponding eigenfunctions. We now focus on Hamiltonians with a chiral pseudosymmetry, described by a local operator ${\hat{C}}$ anticommuting with the Hamiltonian at zero electric field, $\{{\hat{H}}(A_{0}=0,{\mathbf{A}}),{\hat{C}}\}=0$. Microscopically, such a pseudosymmetry can be realized for Hamiltonians defined on a bipartite lattice, with hoppings only between sublattices. Anticommutation implies that eigenstates come in pairs with positive and negative energy $E=E_{\alpha}=-E_{\beta}$, and locality means that the corresponding eigenfunctions have equal probabilities, $|\eta_{\beta}({\mathbf{x}})|^{2}=|{\hat{C}}\eta_{\alpha}({\mathbf{x}})|^{2}=|\eta_{\alpha}({\mathbf{x}})|^{2}$. In particular, each probability can be expressed as half the sum of a positive and a negative energy eigenstate, $|\eta_{\alpha}({\mathbf{x}})|^{2}=\frac{1}{2}(|\eta_{\alpha}({\mathbf{x}})|^{2}+|\eta_{\beta}({\mathbf{x}})|^{2})$. In the absence of zero-energy states, whose consequences will be discussed shortly, $E_{\alpha}\neq 0$, we can apply the above decomposition to Eq. (2), and find that the charge density can be expressed as $1/2$ times the sum over all eigenstates. But this sum corresponds to all bands completely filled, hence the electron density in the absence of an electric field is spatially uniform and has the constant value $\langle{\hat{\rho}}({\mathbf{x}})\rangle_{A_{0}=0,{\mathbf{A}}}=(-e)n/2$.

The charge deviation from a reference configuration without fields is $\delta\langle{\hat{\rho}}({\mathbf{x}})\rangle_{A_{0},{\mathbf{A}}}=\langle{\hat{\rho}}({\mathbf{x}})\rangle_{A_{0},{\mathbf{A}}}-\langle{\hat{\rho}}({\mathbf{x}})\rangle_{A_{0}=0,{\mathbf{A}}=0}$. Being constant, the reference charge density is $1/2$ times the sum over any complete set of states, giving

[TABLE]

In high-energy physics, this expression defines the vacuum charge Jackiw and Schrieffer (1981); Krive and Rozhavskiĭ (1987). Here, an occupied state contributes a total charge of $\delta Q=-e/2$, an empty state contributes $\delta Q=+e/2$, while zero-energy states are ambiguous and can be attributed to either sum. In our setup, the magnetic flux gives localized zero-energy states, and this formula now describes a physical, localized charge deviation $\delta Q=\mp e/2$ whose sign is determined by the occupation of the states selected by the infinitesimal electric field $\mathbf{E}_{z}$.

Analytical results. Since the low-energy states of a TI are localized on the surface 171717See Supplemental Material at … for details regarding time-reversal symmetry and the calculation of the eigenstates in the presence of a thin magnetic flux tube., we consider the Dirac Hamiltonian on a sphere with radius $R$, Imura et al. (2012); Lee (2009); Parente et al. (2011)

[TABLE]

where $\phi\in[0,2\pi)$, $\theta\in[0,\pi]$ are spherical coordinates, and energy is measured in units of $v_{F}/R$ with $v_{F}$ denoting the velocity of Dirac electrons. We have specialized to a vector potential ${\mathbf{A}}$ which has only an azimuthal component ${\mathbf{A}}_{\phi}(\theta)$, reflecting rotational symmetry around the flux tube. We can decompose $\psi(\phi,\theta)=\tilde{\psi}(\theta)e^{im\phi}/\sqrt{R}$ with half-integer angular momentum $m=\pm\frac{1}{2},\pm\frac{3}{2},\dots$ for spin-$1/2$ electrons. In the absence of an external field, $\mathbf{A_{\phi}}=0$, the energy eigenvalues of the spherical Dirac operator are known to be nonzero integers $E=\pm 1,\pm 2,\dots$ whose multiplicities increase with angular momentum Newman and Penrose (1966); Abrikosov (2002). It is tempting to incorporate the flux tube by using the Aharonov-Bohm effect and shifting the angular momentum, $m\to m-N_{0}$ Lee (2009), where $N_{0}=\Phi/\Phi_{0}$ is the number of flux quanta, but this approach does not allow the implementation of the correct boundary condition that the wave function $\tilde{\psi}(\theta)$ stays finite as $\theta\to 0,\pi$ for spatial coordinates within the region of the flux tube. 181818Ref. Lee (2009) does not mention boundary conditions and concludes that even in the absence of an external magnetic field, the spherical Dirac operator has zero-energy states, $E=0$, which does not agree with our results.

To model the flux tube, we express the vector potential as $N_{\phi}(\theta)=eR{\mathbf{A}}_{\phi}(\theta)\sin\theta$ and substitute $x=\cos\theta$. Then, we choose $N_{\phi}(x)=N_{0}\min\{1,(1-|x|)/\delta\}$, which is equal to the total flux $N_{0}$ in most of the sphere, but vanishes at the poles. For a thin flux tube, we have $0<\delta\ll 1$. Using that ${\hat{H}}^{2}=\text{diag}(h^{+}h^{-},h^{-}h^{+})$, we only need to solve the eigenvalue equation $h^{+}h^{-}{\psi_{\uparrow}}=E^{2}{\psi_{\uparrow}}$. Then, the eigenvectors of the original Hamiltonian are obtained as $\psi_{\pm E}=(\pm E{\psi_{\uparrow}},h^{-}{\psi_{\uparrow}})^{T}$ if $E\neq 0$. Zero modes, $E=0$, are obtained from $h^{-}\psi_{\uparrow}=0$, or $h^{+}\psi_{\downarrow}=0$. One finds that

[TABLE]

resembling the Legendre differential operator. This is a special case of the Schrödinger–Lichnerowicz formula Lichnerowicz (1963). We now use a piecewise ansatz Abrikosov (2002); Imura et al. (2012), here shown for angular momentum $m\geq 1/2$:

[TABLE]

where the pieces are defined for $1-\delta\leq x\leq 1$, $|x|<1-\delta$, and $-1\leq x\leq-1+\delta$, respectively. We emphasize that the wave function stays finite near the poles. In coordinates $\xi=(1-x)/2$, the eigenvalue equation is equivalent to a set of hypergeometric equations for the functions $g_{0}(\xi),g(\xi),g_{1}(\xi)$, solved by hypergeometric functions $F(a,b,c;\xi)$ Olver et al. (2016). Abbreviating $c=(m+1/2-N_{0})$ and assuming $c\not\in\mathbb{Z}$, the general solution for the middle part is

[TABLE]

with parameters $\alpha$, $\beta$ that have to be determined by two jump conditions for derivatives, the first one being

[TABLE]

with $\varepsilon=\delta/2$, and the other similar Note (17). To make progress, we now take the limit of a thin flux tube, $\varepsilon\to 0$. Expanding the solutions $g_{0},g_{1}$ to leading order in $\varepsilon$ while ignoring powers of order $\varepsilon^{|c|}$ or higher, we find Note (17) that the jump conditions can only be satisfied for energies $E$ that fulfill

[TABLE]

where $n=0,1,2,\dots$. This spectrum is visualized in Fig. 2(a).

We can now give explicit wave functions for the two eigenstates with energy closest to zero, used in Eq. (1). For flux $0<\Phi<\Phi_{0}$, they have angular momentum $m=1/2$ and are superpositions $\eta_{\mp}=(\eta_{\text{n}}\pm\eta_{\text{s}})/\sqrt{2}$ of states $\eta_{\text{n}}({\mathbf{x}})=[{\psi_{\uparrow}}(x),0]^{T}$ and $\eta_{\text{s}}({\mathbf{x}})=[0,{\psi_{\uparrow}}(-x)]^{T}$ with

[TABLE]

in the region $|x|<1-\delta$ and with normalization $\mathcal{N}=R[\sqrt{\pi}\Gamma(c)/\Gamma(c+1/2)+O(\delta^{c})]$.

Numerical results. To confirm our analytic result, we have performed a numerical calculation using a minimal lattice model for a topological insulator Qi et al. (2008). This model concerns a four-component fermionic wave function on a cubic lattice. In momentum (Bloch) space, the TR-symmetric Hamiltonian is

[TABLE]

where $\lambda,\kappa,t$ are real parameters, and $\sigma_{\mu}$ and $\tau_{\mu}$ are Pauli matrices acting on the spin (respectively, orbital) degrees of freedom. This model is a (strong) topological insulator in the parameter range $2t\leq\kappa\leq 6t$ Qi et al. (2008). Coupling to the electromagnetic field is achieved by Peierls substitution Note (17).

The induced charge in the presence of a flux tube and a small electric field is shown in Fig. 3. The charge distribution $\delta\langle{\hat{\rho}}({\mathbf{x}})\rangle_{A_{0},{\mathbf{A}}}$ is TR symmetric, localized on the surface, and concentrated at the flux tube. To compare with the analytical results, we need to take the limit of a thin flux tube by a scaling analysis where we fix a system size and shrink the size of the flux tube. The electric potential has to be smaller than the level spacing, but larger than the hybridization, as the order of limits is important. We find that the extrapolated value of the charge difference in the top half is $\Delta Q(\Phi_{0})=+(0.49\pm 0.02)e$, in excellent agreement with our analytical results. For a topologically trivial insulator (e.g., $6t<\epsilon$, not shown) no significant charge is accumulated at the flux tube.

Experimental realization. A mesoscopic setup for measuring the TR-symmetric charge response is illustrated in Fig. 4. If we use Bi2Se3 as an example TI Xia et al. (2009) with a Fermi velocity $v_{F}\sim 5\times 10^{5}\,\text{m}\,\text{s}^{-1}$ and assume that the system has a diameter of $1\textmu\text{m}$, then the level spacing (9) of the surface states, $\Delta E=\hbar v_{F}/R$, should be on the order of half a meV. This is well below the bulk band gap $E_{\text{gap}}\sim 0.35\,\text{eV}$, but large enough to be comparable to an externally applied voltage. A thin magnetic flux tube could be generated by using a superconductor to focus the magnetic field, or by pinning two magnetic vortices, or a giant vortex Kanda et al. (2004); Cren et al. (2011). The numerical results, Fig. 3, indicate that larger tube diameters still yield an appreciable charge response, and the value for a thin tube can be obtained by finite-size scaling. The key signature of a TI is that the charge is a half-integer multiple of the elementary charge, as opposed to an integer multiple for purely 2D materials like graphene. Such an experiment is expected to be challenging but within reach of present day technology.

Conclusion. While the electromagnetism of TRS TIs in 3D is commonly associated with the axion action, we have argued that this action is inadequate for describing a physical response, in particular one that is TRS. In search of the latter, we have adapted the idea of charge fractionalization from high-energy physics to a condensed-matter setting. Our main result is that the insertion of a thin flux tube leads to a pair of localized states whose hybridization at one flux quantum is much smaller than the level spacing. Combined with small background electric field, this allows us to adiabatically transform a delocalized eigenstate into a localized eigenstate, giving a TRS charge response of $e/2$ that is, in principle, amenable to experimental detection.

.1 Time-reversal symmetric response in a finite lattice model

In this section, we show that any time-reversal symmetric (TRS) finite lattice model has a TRS electromagnetic response. This demonstrates that the linear magnetoelectric effect vanishes for a TRS topological insulator, and that the axion action cannot not describe a TRS response of a finite system.

We consider any finite lattice, with or without boundary, where sites are labelled by ${\mathbf{x}}$. We consider a non-interacting Hamiltonian ${\hat{H}}(A_{0},{\mathbf{A}})$ in the presence of an electric potential $A_{0}$ and a vector potential ${\mathbf{A}}$. Single-particle wave functions $\psi({\mathbf{x}})$ may have $n$ different components (spin, sublattice, $\dots$) per site. Time-reversal can be represented by an anti-unitary operator ${\hat{T}}=\hat{U}K$, where the unitary operator $\hat{U}$ acts locally on the $n$ internal degrees of freedom, and $K$ denotes complex conjugation. TRS of a Hamiltonian is defined by ${\hat{T}}{\hat{H}}(A_{0},{\mathbf{A}}){\hat{T}}^{-1}={\hat{H}}(A_{0},-{\mathbf{A}})$, i.e.  reversing the direction of electronic motion has the same effect as reversing the magnetic field.

We choose a chemical potential $\mu=0$, such that at zero temperature, the expectation value of the charge density is given by

[TABLE]

where $E_{\alpha}$ are the energy eigenvalues of ${\hat{H}}(A_{0},{\mathbf{A}})$, and $\eta_{\alpha}({\mathbf{x}})$ the corresponding eigenfunctions. Due to the locality property of ${\hat{T}}$ as described above, the probability at every site is unchanged under a time-reversal transformation, that is $|{\hat{T}}\eta_{\alpha}({\mathbf{x}})|^{2}=|\eta_{\alpha}({\mathbf{x}})|^{2}$. For a TRS Hamiltonian, ${\hat{T}}\eta_{\alpha}({\mathbf{x}})$ is an eigenstate of ${\hat{H}}(A_{0},-{\mathbf{A}})$ with energy $E_{\alpha}$, and applying time-reversal to Eq. (12), we find that the charge density is independent of the sign of the magnetic field,

[TABLE]

In other words, the charge response is TRS.

.2 Low-energy states of a bulk TI in the presence of a magnetic field

Usually, the low-energy states of a TI are assumed to be localized on the surface. However, in the literature, a wormhole effect has been discussed Rosenberg et al. (2010), whereby a flux tube threaded through a TI can lead to zero-energy states in the bulk. However, this effect requires that the diameter of the flux tube is much smaller than the lattice spacing, and we now argue that it does not occur when the diameter is larger and a continuum approximation applies.

For simplicity, we consider a geometry of a partially infinite TI that is bounded in $z$-direction by two coordinate planes. Then, the representative continuum bulk Hamiltonian can be written as a sum of two anticommuting self-adjoint operators, ${\hat{H}}=\hat{X}+\hat{Y}$, $\{\hat{X},\hat{Y}\}=0$, where $\hat{X}=\tau_{z}(\sigma_{x}(p_{x}-e{\mathbf{A}}_{x})+\sigma_{y}(p_{y}-e{\mathbf{A}}_{y}))$ and $\hat{Y}=\tau_{z}\sigma_{z}p_{z}+\tau_{x}m(z)$. Here, ${\mathbf{A}}=({\mathbf{A}}_{x},{\mathbf{A}}_{y},0)$ is the vector potential modeling a flux tube along the $z$-axis, $p_{\mu}$ is the momentum operator, $\sigma_{\mu}$, $\tau_{\mu}$ are Pauli matrices acting on spin and valley degrees of freedom and $m(z)$ is the mass profile that interpolates between $-M$ inside and $+M$ outside the material. It follows that $\hat{H}^{2}=\hat{X}^{2}+\hat{Y}^{2}$ and that $\hat{X}^{2}$ and $\hat{Y}^{2}$ commute. Thus, any low-energy state of ${\hat{H}}^{2}$ must be a low-energy state of $\hat{Y}^{2}$, but the latter describes a TI in 1D whose low-energy states are localized on the boundary in $z$-direction. Hence, there are no low-energy bulk states for the 3D Hamiltonian. We expect that this result is generic and also applies to a spherical geometry. It certainly agrees with a numerical calculation of the eigenstates in a cubical geometry, Fig. 2(b) of the main text.

.3 Continuum model: Sphere with flux tube

In this section, we detail the calculation of the eigenstates of the Hamiltonian (4) for a thin flux tube.

First, note that the Dirac Hamiltonian (4) is written in a spin basis that rotates with the angular coordinates. Here, the wave function acquires a phase of $\pi$ when the azimuthal coordinate makes a full rotation. Hence, the wave function has to satisfy antiperiodic boundary conditions, $\psi(\theta,\phi+2\pi)=-\psi(\theta,\phi)$, which means that the angular momentum is half-integral, $m=\pm\frac{1}{2},\pm\frac{3}{2},\dots$. In principle, it would also be possible to use a fixed spinor basis and periodic boundary conditions, but then the Hamiltonian would be different from Eq. (4). However, the latter form of the Hamiltonian informs the boundary conditions, which are that $|\psi(\theta,\phi)|^{2}$ says finite as $\theta\to 0,\pi$.

The Hamiltonian ${\hat{H}}$ is purely off-diagonal, hence its square is a diagonal matrix, ${\hat{H}}^{2}=\text{diag}(h^{+}h^{-},h^{-}h^{+})$. We find

[TABLE]

if the vector potential ${\mathbf{A}}_{\phi}={\mathbf{A}}_{\phi}(\theta)$ does not depend on the azimuthal angle $\phi$. This formula is a special case of the Schrödinger–Lichnerowicz formula Lichnerowicz (1963). The first two terms correspond to the spinor Laplacian, where the term $\frac{1}{2}\cot\theta$ captures the spinorial nature of the wave function. The third term is a Zeeman term and the constant $\frac{1}{2}$ corresponds to the scalar curvature of the sphere.

We now proceed with the eigenvalue problem for the operator (5), $h^{+}h^{-}{\psi_{\uparrow}}(x)=E^{2}{\psi_{\uparrow}}(x)$, where we have substituted $x=\cos\theta$. In any region where the flux is a linear function, $N_{\phi}(x)=Ax+B$, the eigenvalue equation can be brought into the form of a hypergeometric differential equation by making the ansatz ${\psi_{\uparrow}}(x)=(1-x)^{\alpha}(1+x)^{\beta}g(x)$. Choosing the exponents

[TABLE]

and using new coordinates $\xi=(1-x)/2$, we find that the eigenvalue equation is equivalent to

[TABLE]

This has the form of a hypergeometric differential equation Olver et al. (2016)

[TABLE]

with parameters $a,b,c\in\mathbb{R}$. For the positive sign, $\eta=+1$, we can identify the parameter combinations

[TABLE]

Since we have modeled the flux tube by a piecewise linear function $N_{\phi}(x)$, we can make a piecewise ansatz for the wave function. For positive angular momentum $m\geq 1/2$, we use $\eta=+1$ and consider

[TABLE]

The eigenvalue equation for the wave function turns into one hypergeometric equation for each of the functions $g,g_{0},g_{1}$. The choice $\eta=+1$ ensures that the wave function ${\psi_{\uparrow}}(x)$ is finite at the poles if and only if the functions $g_{0}(x)$ and $g_{1}(x)$ are finite at $x=+1$ and $x=-1$, respectively.

The continuity of the wave function ${\psi_{\uparrow}}$ and its derivative ${\partial}_{x}{\psi_{\uparrow}}$ lead to jump conditions for the functions $g_{0},g,g_{1}$ and their derivatives at the points $x=\pm(1-\delta)$. We change coordinates to $\xi=(1-x)/2$, so that we consider functions $g(\xi),g_{0}(\xi),g_{1}(\xi)$ with derivative $g^{\prime}(\xi)={\partial}_{\xi}g(\xi)=-2{\partial}_{x}g(x)|_{\xi=(1-x)/2}$. Setting $\varepsilon=\delta/2$, we find that

[TABLE]

We can now solve Eq. (16) using hypergeometric functions. Let $F(a,b,c;\xi)\equiv{}_{2}F_{1}(a,b,c;\xi)$ denote the Gauss hypergeometric function. We introduce the abbreviation $k=m+\frac{1}{2}=1,2,3,\dots$ and let $c=k-N_{0}$. For simplicity, we consider only the case where the latter quantity is not an integer, $c\not\in\mathbb{Z}$. Then, the general solution is

[TABLE]

where $\alpha_{0},\alpha,\beta,\beta_{0}\in\mathbb{C}$ are complex parameters to be determined, and where we have already used that the solutions $g_{0}(\xi)$ and $g_{1}(\xi)$ have to be finite at the points $\xi=0$ and $\xi=1$, respectively. Here, $w_{1}(\xi)$ and $w_{2}(\xi)$ are two linearly independent solutions of the hypergeometric differential equation. We can also consider two other, linearly independent solutions $w_{3}(\xi)$ and $w_{4}(\xi)$, and write the middle part as

[TABLE]

where the parameters $\tilde{\alpha},\tilde{\beta}$ are related to $\alpha,\beta$ by Kummer’s relations Olver et al. (2016)

[TABLE]

Here, $\Gamma(z)$ is the Gamma function.

The parameters are determined by the jump conditions. Expanding to leading order in $\varepsilon$, we obtain

[TABLE]

To make further progress, we need to know the asymptotic expansions of the quotients $g_{0}^{\prime}(\varepsilon)/g_{0}(\varepsilon)$ and $g_{1}^{\prime}(1-\varepsilon)/g_{1}(1-\varepsilon)$ in the limit $\varepsilon\to 0$. These will be derived in the next section. Using Eqs. (41) and (40), we find that to leading order in $\varepsilon$, the energy $E$ is an eigenvalue if and only if the system of linear equations

[TABLE]

has a non-trivial solution, where $r_{0}$ and $p_{1}$ are coefficients featuring in the asymptotic expansion. We can express $\tilde{\alpha},\tilde{\beta}$ in terms of $\alpha,\beta$. Then, the existence of a solution corresponds to the vanishing of a determinant to leading order in $\varepsilon$. Depending on the sign of $c$, there are two possibilities for what the leading order in $\varepsilon$ is.

Case $c>0$. We have $\varepsilon^{-c}\gg\varepsilon^{0}$. To leading order, the determinant vanishes if

[TABLE]

According to Eq. (40), the second factor is non-zero for $c>0$, so this can only happen if the product of Gamma functions in the denominator has a pole, that is only if

[TABLE]

In this case, the first equation implies that the second parameter is of order $\beta\sim\alpha\varepsilon^{c}$. Hence, the actual solution (22) is dominated by the first term,

[TABLE]

This is actually a polynomial in the variable $\xi$.

Case $c<0$. We have $\varepsilon^{0}\gg\varepsilon^{-c}$. To leading order, the determinant vanishes if

[TABLE]

Since Eq. (41) establishes that $r_{0}\neq 0$ whenever $c<0$, the solutions to this condition are

[TABLE]

In other words, there is a zero-energy mode and a set of equally spaced energy modes. For the latter, we see that the second coefficient $\beta\sim\alpha\varepsilon^{c}$ dominates the first, so the solution has the form

[TABLE]

On the other hand, for the zero-energy mode $E=0$, we see that to leading order, the first equation implies that the parameter $\beta$ has to vanish, while the parameter $\alpha$ can be any constant. If this is the case, then the second equation is satisfied to any order lower than $\varepsilon^{-c}$. Hence, the zero-energy mode of the operator $h^{+}h^{-}$ has the shape

[TABLE]

.4 Asymptotics of hypergeometric functions

In this section, we want to establish asymptotic expansions for the quotients $g^{\prime}_{0}(\varepsilon)/g_{0}(\varepsilon)$ and $g^{\prime}_{1}(1-\varepsilon)/g_{1}(1-\varepsilon)$, Eqs. (21) and (23), in the limit $\varepsilon\to 0$. These quotients feature in the jump conditions for the eigenfunctions.

The derivative of a hypergeometric function is given by $\frac{d}{d\xi}F(a,b,c;\xi)=\frac{ab}{c}F(a+1,b+1,c+1;\xi)$ Olver et al. (2016). Setting

[TABLE]

we find that

[TABLE]

In the limit $\varepsilon\to 0$, there are two cases: If $N_{0}>0$, then we have $a=k+O(\varepsilon)$ and $b=-\frac{N_{0}}{\varepsilon}+k+O(\varepsilon)$, whereas if $N_{0}<0$, these two equations have to be interchanged. In either case, one of the parameters $a,b$ diverges, while the argument $\xi$ of the hypergeometric functions tends to zero. This limit can be expressed in terms of Kummer’s confluent hypergeometric function $M(b,c;\xi)$, which satisfies $\lim_{b\to\infty}F(a,b,c;\xi/b)=M(a,c;\xi)$. Olver et al. (2016) Note that this limit is finite, so the leading term in the asymptotic expansions of the quotients $g_{0}^{\prime}/g_{0}$ and $g_{1}^{\prime}/g_{1}$ must be of order $\varepsilon^{-1}$. In general, we expect that the asymptotic expansion has the form $\frac{p}{\varepsilon}+q+o(1)$ where the expansion coefficients $p,q$ may still depend on the values of $k,N_{0}$ and $E$.

First, let us show that the term of order $\varepsilon^{-1}$ in the expansion of $g^{\prime}_{0}/g_{0}$ at the north pole is precisely $N_{0}/\varepsilon$. If $N_{0}>0$, the parameter $b$ diverges as $-N_{0}/\varepsilon$. We write $F(a,b,k;\varepsilon)\sim F(a,b,k;-N_{0}/b)$ and apply the confluent hypergeometric function, finding

[TABLE]

The last equality follows from the identity $M(n,n;z)=e^{z}$ for any real $n$. Since the hypergeometric function is symmetric in the parameters $a$,$b$, the same argument applies to the case $N_{0}<0$. This establishes the leading order.

For the expansion at the south pole, a similar argument shows that

[TABLE]

In this case, the confluent hypergeometric function can no longer be expressed in terms of elementary functions. We are not so much interested in the precise form of the expansion, but only that the coefficient is distinct from $N_{0}$. We want to establish that it has the form

[TABLE]

For this, we resort to a numerical calculation: The following two plots of $p_{1}(N_{0},k)$ for different values of $k=1,2,\dots$ confirm that the coefficient satisfies $p_{1}\neq N_{0}$ as long as $N_{0}\neq 0$ and that it also satisfies $p_{1}-N_{0}<c$ if $c>0$.

At the north pole, we also need the subleading term. In this case, the limit of the hypergeometric functions is more difficult to compute, and we resort to a numerical computation again. It turns out that the asymptotic expansion depends on the sign of $c=k-N_{0}$, i.e. whether the flux exceeds the angular momentum or not. We have

[TABLE]

This is confirmed by the following plot, which shows the remainder $r_{0}(k,N_{0},E)$ for several angular momenta $k=1,2,\dots$:

Whenever $c<0$, this value is clearly negative, while for $c>0$, it can approach zero. The plot was calculated for a fixed value of $\varepsilon=10^{-4.5}$, which is a good representation of the limit $\varepsilon\to 0$ if the other quantities are significantly larger. Due to numerical difficulties, the value of $r_{0}$ at $E=0$ has to be extrapolated from its values at finite energies, but this can be done visually, as $r_{0}(k,N_{0},E)$ is a constant function of $E$ for $c<0$.

The existence of the zero-energy state, Eq. (33), crucially depends on the fact that the asymptotic expansion contains the non-vanishing term $r_{0}E^{2}$ whenever the flux exceeds the angular momentum, $c<0$.

.5 Lattice model: Coupling to electromagnetic gauge field

We now couple the lattice Hamiltonian (11) from the main text to the electromagnetic gauge field $(A_{0},{\mathbf{A}})$ in order to numerically calculate its charge response.

Peierls substitution provides a minimal gauge-invariant coupling. With link variables ${\mathbf{A}}_{{\mathbf{x}},{\mathbf{x}}+ae_{\mu}}=\int_{{\mathbf{x}}}^{{\mathbf{x}}+ae_{\mu}}{\mathbf{A}}({\mathbf{x}}^{\prime})\cdot d{\mathbf{x}}^{\prime}$, we can introduce covariant lattice derivatives

[TABLE]

where $e_{\mu}$ denotes the unit vector in direction $\mu=1,2,3$ and $a$ is the lattice spacing, which we set to $a=1$. Our single-particle Hamiltonian ${\hat{H}}(A_{0},{\mathbf{A}})$ is then given by

[TABLE]

Here, $\lambda,\kappa,t$ are real parameters and $a^{2}\mathbf{\tilde{B}}_{\mu}({\mathbf{x}})$ is the magnetic flux through a square-shaped surface of side length $a$ centered at the lattice position ${\mathbf{x}}$ and perpendicular to the direction $e_{\mu}$. We imagine that the four components of the fermion wave function are a product of two “orbital” (sublattice) and two “spin” degrees of freedom. Then, the Pauli matrices $\sigma_{\mu}$ and $\tau_{\mu}$ act only on the spin, resp. orbital degrees of freedom. In the parameter range $2t\leq\kappa\leq 6t$, this model describes a (strong) topological insulator Qi et al. (2008).

Note, however, that the coupling to the electromagnetic gauge field is not unique: Different terms can be added as long as they vanish when ${\mathbf{A}}=0$; they correspond to different $g$-factors for the electron. Here, we have chosen a contribution involving the field strength $\mathbf{\tilde{B}}_{\mu}({\mathbf{x}})$ such that the mass term (the term multiplying $\tau_{x}$) commutes with the square of the kinetic term, i.e. with $(\sum_{\mu=1}^{3}\sigma_{\mu}D_{\mu})^{2}$, in the continuum limit ($a\to 0$).