Magnetotransport in Layered Dirac Fermion System Coupled with Magnetic Moments

Dirac materials have attracted much interest because of their intriguing topological characteristics. Unconventional half-integer quantum Hall effect was observed^1,2) in graphene³⁾ due to the Landau level structures of Dirac fermions. Unlike conventional metals, the backward scattering is strongly suppressed⁴⁾ due to the Berry phase π, which makes Dirac fermions extremely high mobility carriers.

Recently first-principles calculations predicted that manganese pnictide SrMnBi₂ is a Dirac material, because of the electronic properties of the Bi square net, and subsequently, its dispersion and Fermi surface were observed by angle-resolved photoemission spectroscopy and quantum oscillations.^5–7) Interestingly, there are structural and physical similarities between SrMnBi₂ and the 112-type iron-based superconductors.⁶⁾ Antiferromagnetic (AF) ordering of the magnetic moments of Mn occurs below 290 K as suggested by the temperature dependences of magnetization, resistivity, and specific heat.⁶⁾ From the perspective of spintronics application,⁸⁾ it is important to investigate the interplay between Dirac fermions and magnetic moments. In this regard, the Dirac material EuMnBi₂, which is isostructural with SrMnBi₂, provides a plausible platform.⁹⁾ Half-integer quantum Hall effect was observed, and the Berry phase π of Dirac fermions was found from the analysis of Shubnikov–de-Haas oscillations.¹⁰⁾

In EuMnBi₂, the layer of Eu²⁺ with spin \(S=7/2\) is closer to the Bi square net than the layer of Mn–Bi edge sharing tetrahedra. AF ordering of Eu moments around \(T_{\text{N}} = 22\) K is suggested from the magnetic susceptibility measurements.⁹⁾ Compared to SrMnBi₂, the Néel temperature associated with the AF ordering of Mn moments is enhanced to 310 K,⁹⁾ which is presumably due to the interaction between Mn moments and Eu moments. Transport measurements demonstrated¹⁰⁾ that the Dirac fermion transport strongly couples to Eu moments. The ordering of Eu moments is AF in the direction of the c-axis, which is perpendicular to the layers of Dirac fermions, and ferromagnetic in the ab plane. Below 120 K, both the in-plane resistivity \(\rho_{xx}\) and the interlayer resistivity \(\rho_{zz}\) show metallic behavior down to \(T_{\text{N}}\) with the large anisotropy of \(\rho_{zz}/\rho_{xx}\sim 480\) at 50 K.¹⁰⁾ A small drop in \(\rho_{xx}\) and an enhancement in \(\rho_{zz}\) were observed at \(T=T_{\text{N}}\). The effect of coupling between Dirac fermions and Eu moments is seen much clearly under a magnetic field. When the magnetic field is applied in the c-axis, \(\rho_{zz}\) increases sharply below \(T_{\text{N}}\), and \(\rho_{zz}/\rho_{xx}\) exceeds 1000% at 9 T, while it is about 180% at 0 T. Spin-flop transition of Eu moments occurs at ∼5.3 T, and there is a steep increase in \(\rho_{zz}\) at the transition point. As the magnetic field was increased, a peak was observed in \(\rho_{zz}\) around 20 T.¹⁰⁾ Remarkably, this peak shows a hysteresis between the field-increasing and field-decreasing runs. A hysteretic anomaly was also observed in \(\rho_{xx}\).

In this letter, we theoretically study the magnetotransport of two-dimensional Dirac fermions coupled with localized moments, and discuss the features experimentally observed in EuMnBi₂. We calculate the in-plane and interlayer conductivities by using the Kubo formula, assuming a phenomenological form of the interlayer tunneling. In order to explain the hysteresis observed in EuMnBi₂, the lift of valley degeneracy is taken into account.

Before investigating the magnetotransport of Dirac fermions, we first consider a model for the layers of magnetic moments in EuMnBi₂. In the absence of a magnetic field, we assume that the localized moments are antiferromagnetically ordered along the c-axis in the ground state. Therefore, we denote the even layers of localized moments as A sublayer and the odd layers of localized moments as B sublayer. As shown in Fig. 1, we introduce the angles \(\theta_{\text{A}}\) and \(\theta_{\text{B}}\) to describe the direction of the localized moments in A sublayer and B sublayer, respectively. The optimum values of \(\theta_{\text{A}}\) and \(\theta_{\text{B}}\) are determined numerically by minimizing the following energy: \begin{align} E &= - K(\cos^{2}\theta_{\text{A}}+\cos^{2}\theta_{\text{B}})-\mu_{B}B(\cos\theta_{\text{A}}-\cos\theta_{\text{B}})\notag\\ &\quad - J\cos(\theta_{\text{A}}-\theta_{\text{B}}). \end{align} (1) Here, K is the anisotropic energy, \(\mu_{B}\) is the Bohr magneton, B is the magnetic field, and \(J>0\) is the AF interaction between localized moments. The magnetic field dependences of \(\theta_{\text{A}}\) and \(\theta_{\text{B}}\) are shown in Fig. 2. In the spin-flop phase, we found that the energy is given by \(E = - J - (\mu_{B}B)^{2}/[2(J - K)]\) and \(\cos\theta_{\text{A}} = -{\cos\theta}_{\text{B}}=\mu_{B}B/[2(J- K)]\). In terms of the spin flop field \(B_{f}\) and the critical field \(B_{c}\), above which the spins are fully polarized, K and J are given by \(K =\mu_{B}B_{f}^{2}/(2B_{c})\) and \(J =\mu_{B}(B_{f}^{2} + B_{c}^{2})/(2B_{c})\), respectively. Substituting the values \(B_{f}\simeq 5\) T and \(B_{c}\simeq 22\) T, which were observed in the experiment,¹⁰⁾ we obtain \(K=0.43\) K and \(J=7.8\) K. In the following, we use this result for the magnetic field dependence of the directions of localized moments.

»View larger version

Figure 1. Layers of localized moments and the definition of the angles, \(\theta_{\text{A}}\) and \(\theta_{\text{B}}\). In each layer, the localized moments are ferromagnetically ordered, while they are antiferromagnetically ordered along the c-axis in the ground state. In between the layers of localized moments, there is a single layer of Dirac fermions, which is not shown in the figure.

»View larger version

Figure 2. (Color online) Magnetic field dependence of the angles of localilzed moments determined by minimizing the energy (1). Because of the anisotropic energy with the parameter K, there is a spin-flop transition at \(B\sim 5\) T. For \(B > 22\) T, the localized moments are fully polarized.

For the purpose of calculating the in-plane and interlayer conductivities, we need the Landau level wave functions of Dirac fermions. For this, we consider a single valley of Dirac fermions. The presence of another valley and the effect of the lift of its degeneracy are examined later. The Hamiltonian is given by \begin{equation} H_{D} = v \begin{pmatrix} 0 &\pi_{x} - i\pi_{y}\\ \pi_{x} + i\pi_{y} &0 \end{pmatrix}, \end{equation} (2) where v is the Fermi velocity of Dirac fermions and \(\pi_{\alpha} = p_{\alpha} + eA_{\alpha}\) (\(\alpha = x, y\)), with \(p_{\alpha}\) and \(A_{\alpha}\) being the momentum operators and the vector potential, respectively. Here, e denotes the electron charge. Based on a first-principles calculation, it was pointed out in Ref. 5 that the Dirac cone in the Bi square net is anisotropic. The largest value of the Fermi velocity is \(1.51\times 10^{6}\) m/s while its smallest value is \(1.91\times 10^{5}\) m/s. Taking the geometrical mean of these values, we assume, \(v=5\times 10^{5}\) m/s.

Taking the Landau gauge, \(\boldsymbol{{A}}=(0,Bx,0)\), we assume a plane wave function in the y-direction. The energy of the Landau levels is given by¹¹⁾ \begin{equation} E_{n} = \text{sgn}(n)v\sqrt{2e\hbar B | n |}, \end{equation} (3) with integer n. The wave function of the Landau level with n is denoted by \(\boldsymbol{{F}}_{n}(\boldsymbol{{r}})\), where \begin{equation} \boldsymbol{{F}}_{0}(\boldsymbol{{r}}) = \begin{pmatrix} 0\\ h_{0}(\boldsymbol{{r}}) \end{pmatrix} \end{equation} (4) and \begin{equation} \boldsymbol{{F}}_{n}(\boldsymbol{{r}}) = \frac{1}{\sqrt{2}} \begin{pmatrix} \text{sgn}(n) h_{| n| -1}(\boldsymbol{{r}})\\ h_{| n |}(\boldsymbol{{r}}) \end{pmatrix} \end{equation} (5) for \(n\neq 0\). The function \(h_{|n|}(\boldsymbol{{r}})\) is given by \begin{align*} h_{| n |}(\boldsymbol{{r}}) &= i^{n}\sqrt{\frac{1}{\sqrt{\pi}2^{| n|}| n| !\ell L}}H_{| n|}\left(\frac{x-k\ell^{2}}{\ell}\right)\\ &\quad \times \exp\left[-\frac{1}{2}\left(\frac{x-k\ell^{2}}{\ell}\right)^{2}\right] \exp(iky). \end{align*} Here, \(H_{| n|}(x)\) is the Hermite polynomial, \(\ell=\sqrt{\hbar/e|B|}\) is the magnetic length, and L is the system size in the y-direction.

The interlayer conductivity is computed by the Kubo formula.¹²⁾ The result is, \begin{equation} \sigma_{zz} = \frac{2t_{c}^{2} a_{c} e^{3}\tau}{\pi\hbar^{3}}B\rho(\mu). \end{equation} (6) Here, \(t_{c}\) is the interlayer hopping, \(a_{c} = 22.5\) Å is the lattice constant in the c-axis,¹⁰⁾ τ is a constant relaxation time, μ is the chemical potential, and \(\rho(\epsilon)\) is the spectral function, which has the form of a Lorentz function. From the analysis of the temperature dependence of the Shubnikov–de Haas oscillation amplitude in \(\rho_{xx}\), τ was estimated as \(\tau = 3.5\times 10^{-14}\) s in SrMnBi₂.⁵⁾ We found that this value is too small for EuMnBi₂, and therefore we assumed that \(\tau = 2.5\times 10^{-13}\) s to reproduce the Shubnikov–de Haas oscillation observed in Ref. 10. For the value of μ, we assume \(\mu=1000\) K to reproduce the Shubnikov–de Haas oscillation period that was observed experimentally in \(\rho_{xx}\) and \(\rho_{zz}\).¹⁰⁾

In order to apply the formula (6) to EuMnBi₂, we need to determine the interlayer hopping parameter, \(t_{c}\). In between the Dirac fermion layers of EuMnBi₂, there are three layers: two layers of Eu moments and one layer composed of Mn²⁺ ions and Bi³⁻ ions. The parameter \(t_{c}\) depends on not only the directions of Eu moments and Mn moments but also the interaction between Dirac fermions and Bi³⁻ ions. We note that those Bi³⁻ ions are different from Bi¹⁻ ions forming the layers of Dirac fermions.^9,10) In order to avoid this complication, here we assume that the magnetic layers between Dirac fermion layers consist of a single component of localized moments, whose AF order is described by \(\theta_{\text{A}}\) and \(\theta_{\text{B}}\), above. Then, in the spin-flop phase, \(t_{c}\) depends on \(\theta_{\text{B}} =\pi -\theta_{\text{A}}\). Apparently, the interlayer hopping \(t_{c}\) is maximum when the localized moments are fully polarized, while it is minimum in the AF phase. As a simple phenomenological model, we assume the following form \begin{equation} t_{c} \propto 1-\cos\frac{\theta_{\text{B}}}{2}. \end{equation} (7) The proportionality constant is chosen such that \(t_{c}\) takes the value of 200 K in the fully polarized state of localized moments.

The Kubo formulae for the in-plane conductivity and Hall conductivity are given by \begin{align} \sigma_{xx} &= \frac{e^{2}\hbar}{2 \ell^{2}a}\sum_{n,n'} \rho(\mu-E_{n})\rho(\mu-E_{n'})|\langle \boldsymbol{{F}}_{n}| \sigma_{x}| \boldsymbol{{F}}_{n'}\rangle|^{2}, \end{align} (8) \begin{align} \sigma_{xy} &= \frac{\hbar e^{2}v^{2}}{\pi a \ell^{2}}\sum_{n,n'} \frac{f(\mu-E_{n})}{(E_{n}-E_{n'})^{2}}\text{Im}\langle\boldsymbol{{F}}_{n}| \sigma_{x}| \boldsymbol{{F}}_{n'}\rangle\langle\boldsymbol{{F}}_{n'}| \sigma_{y}| \boldsymbol{{F}}_{n}\rangle, \end{align} (9) respectively. Here, \(a=4.5\) Å is the lattice constant in the plane,⁹⁾ and \(\sigma_{x}\) and \(\sigma_{y}\) are Pauli matrices. The function \(f(x)\) is defined by \(f(x)=1/2+\arctan(\frac{x}{\Gamma})/\pi\), where \(\Gamma=\hbar/\tau\) is a constant parameter describing the broadening of Landau levels due to impurity scatterings. The resistivity is computed numerically using the formula given by \(\rho_{xx}=\sigma_{xx}/(\sigma_{xx}^{2}+\sigma_{xy}^{2})\), \(\rho_{yx} =\sigma_{xy}/(\sigma_{xx}^{2} +\sigma_{xy}^{2})\), and \(\rho_{zz}=1/(\sigma_{zz}+\sigma_{0})\). Here, parameter \(\sigma_{0}=0.01\), which is associated with impurity scattering, is introduced to reproduce the experimentally observed \(\rho_{zz}\).¹⁰⁾

The result obtained by including the Zeeman energy splitting is shown in Fig. 3. Here, the magnetic field is normalized by introducing parameter \(B_{F}\equiv \mu^{2}/(2e\hbar v^{2})\simeq 23\) T, which is the frequency of Shubnikov–de-Haas oscillations. The inverse of the Hall resistivity \(\rho_{yx}\) shows a half-integer quantum Hall effect, as shown in Fig. 3(a). Some deviations from the ideal quantum Hall plateaus occur due to the broadening factor, Γ, and similar features were seen in the experiment.¹⁰⁾ \(\rho_{xx}\) and \(\rho_{zz}\) are shown in Figs. 3(b) and 3(c), respectively. The positions of the peaks of \(\rho_{xx}\) and the minimum of \(\rho_{zz}\) reflect the Landau level structure.¹⁰⁾ The density of states takes a large value when the chemical potential is equal to a Landau level energy upon varying the magnetic field. Reflecting this, \(\rho_{zz}\) shows a dip while \(\rho_{xx}\) shows a peak around integer values of \(B/B_{F}\). The splitting of peaks in \(\rho_{xx}\) and that of dips in \(\rho_{zz}\) are due to the lift of the level degeneracy. Here, the degeneracy is associated with the spin degrees of freedom. The result shown in Fig. 3 qualitatively agrees with the experiment.¹⁰⁾ In the experiment, the splitting of the peak in \(\rho_{xx}\) is observed at \(B_{F}/B=1\). Furthermore, the peak at \(B_{F}/B<1\) is larger than that at \(B_{F}/B>1\). This is consistent with the experiment as well. In the experiment, the splitting of the peak in \(\rho_{xx}\) is not observed at \(B_{F}/B=2\). However, the second derivative of \(\rho_{xx}\) with the minus sign, obtained from the experimental data, shows the splitting of the peak. The positions of the minimum in \(\rho_{zz}\), shown in Fig. 3(c), is also consistent with the experiment.¹⁰⁾ However, there is an important difference: the theory suggests the presence of two dips and one small peak around \(B_{F}/B=1\), whereas this peak is much larger in the experiment. It is unlikely that the peak appears as a result of splitting of the dip.

»View larger version

Figure 3. (Color online) Normalized inverse Hall resistivity \(\rho_{yx}^{0}\rho_{yx}/\) (a), in-plane resistivity \(\rho_{xx}\) (b), and interlayer resistivity \(\rho_{zz}\) (c) versus the inverse of the normalized magnetic field. Here, we take \(\rho_{yx}^{0} = 1.1\times 10^{-3}\) Ω·cm. \(B_{F}\) is the frequency of Shubnikov–de-Haas oscillations defined in the main text. Here, we plot the values in \(B<B_{c}\) where \(t_{c}\), which is given by Eq. (7), is finite. In the panel (a) we plot \(\rho_{yx}\) for different values of Γ and the idealized case is shown in the dashed line.

Now we consider the effect of valley splitting. The Hamiltonian is given by \begin{equation} H = v \begin{pmatrix} \tilde{\Delta} &\pi_{x} - i\pi_{y} &0 &0\\ \pi_{x} + i\pi_{y} &\tilde{\Delta} &0 &0\\ 0 &0 &- \tilde{\Delta} &\pi_{x} + i\pi_{y}\\ 0 &0 &\pi_{x} - i\pi_{y} &- \tilde{\Delta} \end{pmatrix}. \end{equation} (10) Here, \(\tilde{\Delta}=\Delta/v\). The parameter Δ is the energy gap created by the valley splitting. The Landau level wave functions for another valley are obtained by simply multiplying \(\tau_{x}\), which is the Pauli matrix in the sublattice space in the layer of Dirac fermions, to the two-component spinors, Eqs. (4^,) and (5^,), from the left hand side. Experimentally, it is suggested that Δ depends linearly on the magnetic field.¹³⁾ We assume that the valley splitting occurs through the interaction with the localized moments, and the origin of the hysteresis is in the system of localized moments. Therefore, we use different proportionality constants for the field-decreasing run and the field-increasing run. We use 1.6\(\mu_{B}\) for the former and 1.5\(\mu_{B}\) for the latter. We used these values to reproduce the two-peak structure of \(\rho_{xx}\) that was observed around \(B_{F}/B = 1\) in the experiment.¹⁰⁾ The difference of 0.1\(\mu_{B}\) in these values is taken such that the numerical calculation reproduces the experimentally observed difference in peak values, which are associated with the hysteresis, of \(\rho_{xx}\) at \(B_{F}/B\sim 1.2\). The numerical calculation result is shown in Fig. 4. The difference between the two cases is seen in \(\rho_{xx}\) around \(B_{F}/B\sim 1\), which is similar to that observed in the experiment.¹⁰⁾ However, the theory failed to reproduce the experimentally observed feature in \(\rho_{zz}\) around \(B_{F}/B\sim 1\).¹⁰⁾ This discrepancy is probably due to the form of the interlayer hopping Eq. (7^,). We need a more realistic \(t_{c}\) than Eq. (7) to explain the experiment, which is left for future research.

»View larger version

Figure 4. (Color online) In-plane resistivity \(\rho_{xx}\) and interlayer resistivity \(\rho_{zz}\) (inset) versus the inverse of the normalized magnetic field. The green lines indicate the result with \(\Delta = 1.6\)\(\mu_{B}\), and the blue lines indicate the result with \(\Delta = 1.5\)\(\mu_{B}\).

An important question, which needs to be considered, is the mechanism of valley splitting. If we recall that bismuth plays a major role in topological insulators^14,15) and the key is strong spin–orbit coupling, one possible scenario is that valley splitting occurs through spin–orbit coupling, \(i{\boldsymbol{{\lambda}}}\cdot {\boldsymbol{{\sigma}}}\) with \({\boldsymbol{{\lambda}}}=(\lambda_{x},\lambda_{y},\lambda_{z})\). Generally, \({\boldsymbol{{\lambda}}}\) depends on the wave vector of Dirac fermions; however, if we focus on a single Dirac point, then the Hamiltonian is given by \begin{equation} H_{\lambda} = \begin{pmatrix} - \mu_{B}B_{\text{eff}} &0 &i\lambda_{z} &i\lambda_{x} + \lambda_{y}\\ 0 &\mu_{B}B_{\text{eff}} &i\lambda_{x} - \lambda_{y} &- i\lambda_{z}\\ - i\lambda_{z} &- i\lambda_{x} - \lambda_{y} &- \mu_{B}B_{\text{eff}} &0\\ - i\lambda_{x} + \lambda_{y} &i\lambda_{z} &0 &\mu_{B}B_{\text{eff}} \end{pmatrix}. \end{equation} (11) Here, we have included in B, the effect of the molecular fields created by localized spins, and we denote the effective magnetic field by \(B_{\text{eff}}\). The eigenvalues of this Hamiltonian is obtained exactly as, \(E^{(\pm,\pm)} =\pm\sqrt{\smash{(\lambda_{z}\pm \mu_{B}B_{\text{eff}})^{2} +\lambda_{x}^{2} +\lambda_{y}^{2}}\mathstrut}\). For \(|\lambda_{z}|\gg |\mu_{B} B_{\text{eff}}|\), we obtain, \(E^{(\pm,\pm)}\simeq \pm\lambda \pm\frac{\lambda_{z}}{\lambda}\mu_{B}B_{\text{eff}}\). Therefore, in this scenario, the parameter Δ is given by \begin{equation} \Delta = \frac{2\lambda_{z}}{\lambda}\mu_{B}B_{\text{eff}}. \end{equation} (12) If \(2\lambda_{z}/\lambda\sim 1\), then the above calculation is justified. Furthermore, the slight change in Δ in the field-decreasing run and the field-increasing run is associated with the hysteresis in the system of localized moments.

To conclude, we have investigated the magnetotransport of two-dimensional Dirac fermions coupled with localized moments, and compared the theoretical result with the experimental result for EuMnBi₂. Most of the features observed in EuMnBi₂ is understood by our model with interlayer hopping Eq. (7). The hysteresis observed in \(\rho_{xx}\) and \(\rho_{zz}\) can be associated with the valley splitting resulting from the spin–orbit coupling and the coupling between Dirac fermions and localized moments. However, the theory failed to explain the feature around the Landau level with the index \(n=1\). Presumably, this discrepancy arises from the complicated interaction between Dirac fermions and the magnetic layer in EuMnBi₂.

References

• 1 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005). 10.1038/nature04233 Crossref, Google Scholar
• 2 Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005). 10.1038/nature04235 Crossref, Google Scholar
• 3 A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007). 10.1038/nmat1849 Crossref, Google Scholar
• 4 T. Ando, T. Nakanishi, and R. Saito, J. Phys. Soc. Jpn. 67, 2857 (1998). 10.1143/JPSJ.67.2857 Link, Google Scholar
• 5 J. Park, G. Lee, F. Wolff-Fabris, Y. Y. Koh, M. J. Eom, Y. K. Kim, M. A. Farhan, Y. J. Jo, C. Kim, J. H. Shim, and J. S. Kim, Phys. Rev. Lett. 107, 126402 (2011). 10.1103/PhysRevLett.107.126402 Crossref, Google Scholar
• 6 J. K. Wang, L. L. Zhao, Q. Yin, G. Kotliar, M. S. Kim, M. C. Aronson, and E. Morosan, Phys. Rev. B 84, 064428 (2011). 10.1103/PhysRevB.84.064428 Crossref, Google Scholar
• 7 G. Lee, M. A. Farhan, J. S. Kim, and J. H. Shim, Phys. Rev. B 87, 245104 (2013). 10.1103/PhysRevB.87.245104 Crossref, Google Scholar
• 8 I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004). 10.1103/RevModPhys.76.323 Crossref, Google Scholar
• 9 A. F. May, M. A. McGuire, and B. C. Sales, Phys. Rev. B 90, 075109 (2014). 10.1103/PhysRevB.90.075109 Crossref, Google Scholar
• 10 H. Masuda, H. Sakai, M. Tokunaga, Y. Yamasaki, A. Miyake, J. Shiogai, S. Nakamura, S. Awaji, A. Tsukazaki, H. Nakao, Y. Murakami, T. H. Arima, Y. Tokura, and S. Ishiwata, Sci. Adv. 2, e1501117 (2016). 10.1126/sciadv.1501117 Crossref, Google Scholar
• 11 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). 10.1103/RevModPhys.81.109 Crossref, Google Scholar
• 12 T. Osada, J. Phys. Soc. Jpn. 77, 084711 (2008). 10.1143/JPSJ.77.084711 Link, Google Scholar
• 13 H. Masuda, private communication. Google Scholar
• 14 M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010). 10.1103/RevModPhys.82.3045 Crossref, Google Scholar
• 15 X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011). 10.1103/RevModPhys.83.1057 Crossref, Google Scholar