Electronic chiralization as an indicator of the anomalous Hall effect in unconventional magnetic systems

We first discuss the symmetry properties of the AHE of a crystal. The AHE is described by the anomalous Hall (pseudo)vector $(\bm{\sigma}_{\rm AH})_{\alpha}=\frac{1}{2}\epsilon_{\alpha\beta\gamma}\sigma_{\beta\gamma}$. $\bm{\sigma}_{\rm AH}$ changes sign under time reversal (TR, acting on the equilibrium state of the system) as a consequence of the Onsager relation; it rotates as a pseudovector under O(3) and is invariant under continuous translation operations because it is the response of a uniform current to a uniform electric field.

It is obvious why $\bm{\sigma}_{\rm AH}$ must transform as a vector under proper rotations. By transform we mean the comparison between $\bm{\sigma}_{\rm AH}$ of two systems that are related to each other by such a transformation. Suppose that system 2 is obtained from system 1 by rigidly rotating the former. Then $\bm{\sigma}_{\rm AH}$ for system 2 is measured by applying the external electric field and attaching the current probes in the same manner as that for 1. In other words, we have the response relation

and the linear response function is defined as

However, because the rotation is equivalent to a coordinate transformation, we immediately know that if the electric field is rotated together with the system, the result should not change. In other words

where $R$ is a rotation matrix. Therefore

It then follows that $\bm{\sigma}_{\rm AH,2}=R\bm{\sigma}_{\rm AH,1}$. A similar argument can be made for improper rotation and continuous translation. Alternatively one can use linear response theory and show that all unitary transformations due to O(3) operations cancel out when taking the trace.

In addition to the system-independent properties above, $\bm{\sigma}_{\rm AH}$ must also be invariant under any symmetry operations of the crystal, as dictated by Neumann’s principle. In particular, for symmetry operations that combine a point group operation $R$ with a spatial translation $T$, since $\bm{\sigma}_{\rm AH}$ is invariant under continuous translation, it must be invariant under $R$ even if $R$ is not a symmetry of the crystal.

Based on the discussion above, a suitable indicator of the AHE should be (1) a TR-odd pseudovector, and (2) invariant under all symmetry operations of the crystal. Then $\bm{\sigma}_{\rm AH}$ will be linearly dependent on this indicator to the lowest order of the latter. An important consequence of (1) is that the indicator must be translationally invariant.

The TR-odd property of $\bm{\sigma}_{\rm AH}$ is fundamentally due to the microscopic magnetization density $\mathbf{m}(\mathbf{r})$ in equilibrium. In ferromagnets the spatial average of $\mathbf{m}(\mathbf{r})$, $\bar{\mathbf{m}}$ serves as a suitable indicator of the AHE. When $\bar{\mathbf{m}}$ nearly vanishes, it is reasonable to associate the AHE with the spatial variation of $\mathbf{m}(\mathbf{r})$. We thus propose indicators of the AHE constructed from the spatial gradients of $\mathbf{m}(\mathbf{r})$. For definiteness we only consider the first-order spatial derivative of $\mathbf{m}(\mathbf{r})$ in this work, although indicators based on higher orders in $\mathbf{m}$ or its derivatives can be constructed similarly and may be useful in different cases. We start from a Cartesian tensor $T_{ijk}$ defined as

where $\phi$ is a TR-even scalar field observable of the crystal, which can be the charge density $\rho(\mathbf{r})$ or the nonmagnetic potential $V(\mathbf{r})$; uc stands for unit cell. We ignore any boundary contributions to $\mathbf{m}(\mathbf{r})$ and $\phi(\mathbf{r})$ so that they have the same discrete translation symmetry as the infinite crystal. The inclusion of $\phi$ ensures that $T_{ijk}$ does not become a boundary term and also signifies the role of orbital degrees of freedom in the AHE. $T_{ijk}$ is a TR-odd rank-3 pseudotensor and is translationally invariant. It is also invariant under any symmetry operations of the crystal since both $\phi(\mathbf{r})$ and $\mathbf{m}(\mathbf{r})$ are physical observables of the crystal. A pseudovector can be obtained from $T_{ijk}$ by contracting it with Kronecker $\delta$ or Levi-Civita symbol $\epsilon$. The only two independent pseudovectors obtained from this construction are

We name $\bm{\chi}_{1,2}$ generally as “electronic chiralization” to emphasize their electronic origin and pseudovector nature, analogous to the optical chirality (flow) in optics [33, 34, 32, 35]. Several comments are in order:

In this subsection we discuss how to efficiently calculate EC in realistic systems. Since crystallographic models are usually represented by localized atomic charge and magnetic moments, we consider such cases first by assuming that the atomic charge and spin densities are described by Gaussians, which gives

where $Q_{n}$ and $\mathbf{M}_{n}$ are the charge and magnetic moment on a lattice site located at $\mathbf{r}_{n}$ in the unit cell, $\mathbf{R}$ stands for Bravais lattice vectors, and $g(\mathbf{r})=(2\pi\sigma^{2})^{-\frac{3}{2}}e^{-\frac{r^{2}}{2\sigma^{2}}}$ is the Gaussian function whose Fourier transform is $g_{\mathbf{k}}=e^{-\frac{\sigma^{2}k^{2}}{2}}$. The Fourier transform of $\rho(\mathbf{r})$ and $\mathbf{m}(\mathbf{r})$ are therefore

from which we can obtain $\bm{\chi}_{1,2}$

where

Because of the $g_{\mathbf{K}}^{2}$ in Eq. (10) the summation will converge quickly if $\sigma$ is comparable to the lattice constant.

Although the momentum space expressions above can be directly applied to scattering data (with the $g_{\mathbf{K}}$ replaced by appropriate form factors) and plane-wave-based DFT calculations (see Sec. III.2 for an example), they do not necessarily provide a transparent picture on the real-space charge and magnetization distributions. Moreover, when the charge and magnetization densities are very localized around individual atoms, and have fine structures at small length scales (or equivalently at high energy scales) that can be irrelevant to transport phenomena, the values of $\bm{\chi}_{1,2}$ can be sensitive to the cutoff in $\mathbf{K}$. We therefore discuss real-space expressions of $\bm{\chi}_{1,2}$ next, still assuming that the local densities for each atom are described by Gaussians.

For $\bm{\chi}_{1}$ we have

in which the integral

Because of the Gaussian factor in $\overleftrightarrow{I}$ one can potentially truncate the summation in Eq. (12) at e.g. nearest neighbor. Note that the second term in $\overleftrightarrow{I}$ has the form of an electric quadrupole. The main difference between it and the quadrupole moment in [23] is that the origin of the former is different for different pairs of ions.

Using the same $\overleftrightarrow{I}$ one can express $\bm{\chi}_{2}$ as

The contributions due to the traceless part of $\overleftrightarrow{I}$ to $\bm{\chi}_{1,2}$ are therefore the same. The trace of $\overleftrightarrow{I}$ gives rise to a weighted sum of electric charge surrounding a given magnetic moment $\mathbf{M}_{m}$. A compensated ferrimagnet will have a nonzero $\bm{\chi}_{1}$ mainly due to this contribution. We will focus the traceless part of $\overleftrightarrow{I}$ in Sec. III below since it is less trivial. This also makes it sufficient to consider $\bm{\chi}_{1}$ only. Real-space formulas for the generalized EC introduced in Eq. (7) are given in Appendix A.

We first calculate $\bm{\chi}_{1}$ using crystallographic models for the noncollinear AHE antiferromagnets Mn₃X (X = Ir, Pt, Sn, Ge, etc.) and the real-space formula Eq. (12). Considering nearest neighbors and the traceless part of $\overleftrightarrow{I}$ only, we have $\bm{\chi}_{1}=C\sum_{m}\mathbf{M}_{m}\cdot\mathcal{Q}_{m}$, where $C=-e^{\frac{r_{\rm nn}^{2}}{4\sigma^{2}}}/(32\pi^{\frac{3}{2}}\sigma^{7}V_{\rm uc})$ is a system dependent constant and $r_{\rm nn}$ is the distance between a magnetic atom and its nearest neighbors. $\mathcal{Q}_{m}$ is the total electric quadrupole relative to the position of site $m$:

where $\{i\}_{m}$ stands for the set of nearest neighbors of site $m$. $\mathcal{Q}_{m}$ has the same symmetry as a second-order, i.e. easy-axis or easy-plane, magnetic anisotropy. The origin of weak ferromagnetism and hence the AHE in systems with a nonzero $\bm{\chi}_{1}$ is thus a site-dependent second order anisotropy, which applies to the known examples of Mn₃X and collinear AHE antiferromagnets [20, 21]. If such a site-dependent second order anisotropy is forbidden by symmetry, as in the case of hematite, one needs to consider the generalizations of EC as mentioned in Sec. II.3.

For the structure of cubic Mn₃X, e.g. Mn₃Ir [Fig. 1 (a)], $\mathcal{Q}_{m}$ is diagonal with the principal axis of the largest eigenvalue along the four-fold axis on each Mn atom. One can then obtain the dependence of $\bm{\chi}_{1}$ on rigid rotations of the sublattice moments, which is very different from rotating the pseudovector $\bm{\chi}_{1}$ directly since the latter is not required to transform as a pseudovector under separate rotations of the lattice and magnetic moments. For example, the length of $\bm{\chi}_{1}$ depends on the rotation about its direction as $|\cos\gamma|$ [Fig. 1 (b)], where $\gamma$ is the rotation angle, similar to $\bm{\sigma}_{\rm AH}$ (see below) and the orbital magnetization [50]. For the structure of hexagonal Mn₃X, e.g. Mn₃Sn [Fig. 1 (c)], one can follow the same procedure and obtain a compact expression:

where $\alpha,\beta,\gamma$ are Euler angles about the $z,y^{\prime},z^{\prime\prime}$ axes, respectively. Here $z$ and $y$ are the original crystalline axes $[0001]$ and $[01\bar{1}0]$, while the primed and double-primed axes correspond to those co-rotated with the local moments successively by $\alpha$ and $\beta$, respectively. The expression captures the phenomenon that a counterclockwise rotation about $[0001]$ of all sublattice moments leads to clockwise rotation of the weak magnetization and $\bm{\sigma}_{\rm AH}$ [12]. Moreover, rotation about the $[01\bar{1}0]$ axis in Fig. 1 (c) by $\pi$ makes $\bm{\chi}_{1}$ vanish [Fig. 1 (d)], since in this case the magnetic order becomes triangular rather than inverse triangular, and the AHE or weak ferromagnetism are forbidden by a $C_{3}$ symmetry.

To compare the angular dependence of $|\bm{\chi}_{1}|$ with that of $|\bm{\sigma}_{\rm AH}|$, we use the following generic tight-binding model [50] adapted to the structures of cubic and hexagonal Mn₃X to calculate the intrinsic contribution to the anomalous Hall conductivity:

where $i,j$ label lattice sites, $\langle\rangle$ means nearest neighbor, $\alpha,\beta$ label spin, $t>0$ is the spin-independent hopping amplitude and is chosen as the energy unit, $\lambda_{\rm so}$ is the spin-orbit coupling strength, $\hat{r}_{ij}$ is a unit vector along the position vector $\mathbf{r}_{j}-\mathbf{r}_{i}$, $\bm{\eta}_{ij}$ is the electric field or electric dipole moment vector at the center of the nearest-neighbor $ij$ bond [50] (normalized using the largest $|{\eta}_{ij}|$), $J$ is the strength of a local exchange field along $\hat{n}_{i}$ mimicking the noncollinear magnetic order in a given material.

The anomalous Hall conductivity is calculated as

where $\sigma^{\alpha\beta}(\omega)$ is the optical conductivity tensor, $N_{\mathbf{k}}$ is the number of points of the $k$-mesh, and $\eta$ is a band broadening parameter that depends on disorder. To facilitate rapid convergence of the Brillouin zone integration we have used $\eta=0.1$ (in units of $t$) and a thermal smearing with $k_{\rm B}T=0.3$. Such smearing parameters also help to eliminate any spurious abrupt changes of the zero-temperature intrinsic AHC in a perfect crystal versus smooth changes of tuning parameters such as the rotation angle, since the former is sensitively dependent on small band splittings at the Fermi energy. In reality the dependence of transport coefficients on orientations of the magnetic order parameter are expected to consist of low-order Fourier components due to both thermal and disorder effects. The Brillouin zone integration is performed using a $31\times 31\times 31$ mesh.

Figures. 2 (b) and (d) agree qualitatively with the angular dependence of the EC in Figs. 1 (b) and (d). Such an agreement to some extent depends on the parameter values used, but is generally expected based on symmetry arguments: For an arbitrary rotation of the local magnetic moments along a closed path one can generally expand the $\bm{\sigma}_{\rm AH}$ and EC as Fourier series. The high-symmetry points on the rotation path, either corresponding to exact vanishing of the $\bm{\sigma}_{\rm AH}$ and EC (e.g., $\pi$ rotation about $[111]$ in Mn₃Ir), or equivalent to certain magnetic space group operations (e.g., $2\pi/3$ rotation about $[0001]$ in Mn₃Sn), place identical constraints on the Fourier coefficients of $\bm{\sigma}_{\rm AH}$ and EC. As a result, when the angular dependence is smooth so that only a few low-order Fourier coefficients are relevant, it is expected that $\bm{\sigma}_{\rm AH}$ and EC should behave similarly. Such a symmetry analysis is analogous to that commonly used in ferromagnetic crystals [51]. Nonetheless, $\bm{\sigma}_{\rm AH}$ and EC are not required to have the same angular dependence, which is also similar to the relation between $\bm{\sigma}_{\rm AH}$ and the net magnetization in a ferromagnet.

In this subsection we give an example of a microscopic calculation of $\bm{\chi}_{1}$ by using a 2D Rashba-type continuum model. The procedure can be straightforwardly applied to plane-wave-based density functional theory calculations. The model also serves as an example of the applicability of EC for conventional ferromagnets. The Hamiltonian describes itinerant electrons subject to periodic scalar, Zeeman, and Rashba spin-orbit coupling potentials on a square lattice:

where $\lambda_{R}$ is the strength of the Rashba spin-orbit coupling, $J$ is that of a Zeeman field along $\hat{z}$, and $V$ is that of a scalar confinement potential. $\{,\}$ stands for anti-commutation and is needed to ensure that the position-dependent Rashba term is Hermitian. $G(\mathbf{r})$ is a periodic Gaussian-like function:

where uc stands for unit cell defined by the lattice vectors $\mathbf{a}_{1}=a\hat{x},\mathbf{a}_{2}=a\hat{y}$, $\mathbf{K}$ are reciprocal lattice vectors, and $g(\mathbf{r})$ is the 2D Gaussian $g(\mathbf{r})=\frac{1}{2\pi\sigma^{2}}e^{-\frac{r^{2}}{2\sigma^{2}}}$. In practice we will set $G_{\mathbf{K}}=0$ when $\max(|K_{x}|,|K_{y}|)>K_{G}$, $K_{G}$ being a parameter. We choose $a$ as the length unit, and $E_{0}=\hbar^{2}/(2ma^{2})$ as the energy unit. An illustration of $G(\mathbf{r})$ with $K_{G}=10\left(\frac{2\pi}{a}\right)$ is shown in Fig. 3 (a). In the plane wave basis an arbitrary eigenfunction can be written as

where $c$ are $2\times 1$ column vectors. Substituting this wavefunction into the eigen-equation $H\psi=\epsilon\psi$ for the dimensionless Hamiltonian and using the orthogonality between plane waves we obtain

where

Equation (22) represents infinite coupled linear equations for a given $\mathbf{k}$, or a matrix equation for the column vector $c_{\mathbf{K}}(\mathbf{k})$, with $\mathbf{K}$ understood as a row or column index. We truncate the Hamiltonian matrix by requiring $\max(|K_{x}|,|K_{y}|)\leq K_{H}$. Note that $K_{H}$ is generally different from $K_{G}$ defined above. The dimension of the Hamiltonian matrix for our model is therefore $N\times N$, $N=2(2K_{H}+1)^{2}$ ($K_{H}$ is in units of $2\pi/a$). Figure 3 (b) shows a typical band structure of the model [$K_{H}=10\left(\frac{2\pi}{a}\right)$, only the lowest six bands are plotted].

The spin density (at position $\bm{\tau}$) operator has the following matrix elements (taking $\frac{\hbar}{2}$ as the units of spin)

whose Fourier transform at reciprocal lattice vector $\mathbf{K}_{0}$ is

The expectation value of the spin density for a given chemical potential is therefore

where $f_{n\mathbf{k}}$ is the Fermi-Dirac distribution function at eigen-energy $\epsilon_{n\mathbf{k}}$. Its Fourier transform at reciprocal lattice vector $\mathbf{K}_{0}$ is

The ferromagnetic Rashba model [49] is known to have the AHE with an out-of-plane $\bm{\sigma}_{\rm AH}$ [3]. Therefore we only consider the out-of-plane components of the EC. Taking the spatial dependence of the Rashba coefficient as $\partial_{z}\phi(\mathbf{r})$ in the definition of $\bm{\chi}_{1}$, we can finally obtain

Moreover, the spatial dependence of the Gaussian potentials allows us to calculate $\bm{\chi}_{2}$ and $\bm{\chi}_{3}$ as well, for which we will directly use the electron density $\rho(\bm{\tau})$ as the scalar field, whose expressions are similar to Eqs. (24)–(27) but with $\bm{\sigma}$ replaced by $\sigma_{0}$. Then

Note that in the present case $\bm{\chi}_{3}$ is not a linear combination of $\bm{\chi}_{1,2}$ any more.

Figures 3 (c-f) show representative results of the spin densities and the EC from this model. The in-plane spin components have a nonzero divergence near the center of the unit cell, and become large when the gradient of the confinement potential is significant, which leads to a finite $\chi_{1}^{z}$. The out-of-plane spin component is largest near the center of the unit cell, where the Laplacian of the confinement potential is also largest, leading to a finite $\chi_{2}^{z}$. The summands of Eqs. (28) and (29) are plotted in Figs. 3 (d) and (f), respectively. The bright dots represent the (Bragg) peaks in diffraction experiments and the color coding in log scale suggests fast decay versus increasing crystal momentum. The final results are $\chi_{1}^{z}=0.02685$, $\chi_{2}^{z}=-0.03058$ [in the units implied in Eqs. (28) and (29)] for the parameter values listed in the caption of Fig. 3.

We first consider a minimal tight-binding model having the essential ingredients for magnetic-charge-induced AHE. The model describes $s$ electrons hopping between nearest neighbors on a honeycomb lattice, with magnetic charge of opposite signs residing on the two sublattices [Fig. 4 (a)]:

where the four terms respectively correspond to spin-independent hopping, spin-dependent hopping due to the magnetic charge, Rashba spin-orbit coupling, and an on-site potential breaking the sublattice symmetry; $\eta_{ij}=+1(-1)$ if $\hat{r}_{ij}$ points from sublattice A (B) to B (A); $\eta_{ij}$ together with $\bm{\sigma}\cdot\hat{r}_{ij}$ capture the spin-dependent hopping due to the magnetic field ($\mathbf{H}$ field) lines between neighboring magnetic charges; $\gamma_{i}=+1(-1)$ on A (B) sublattice. The Rashba term is needed to provide the direction of $\bm{\chi}_{1}$ along $z$ as suggested by the expression of $\bm{\chi}_{1}$, and the sublattice potential is needed to break the degeneracy between the opposite magnetic charges on the two sublattices.

Figure 4 (b) shows the typical band structure and Berry curvature of model (IV.1). The Berry curvature is generally nonzero and is largest near the Brillouin zone corners $K,K^{\prime}$. However, the Berry curvatures at opposite momenta do not cancel each other due to the broken time-reversal symmetry by the $t_{M}$ term, which is also evident from the 2D plot in Fig. 4 (c). Therefore the AHE is generally nonzero, except when the Fermi energy is in the gap opened by $\Delta$ which we explain next.

After Fourier transform the toy model becomes [for convenience we rotate the honeycomb lattice in Fig. 4 (a) clockwisely by $\pi/6$]

where $C_{\mathbf{k}}^{\dagger}=(c_{\mathbf{k}A\uparrow}^{\dagger},c_{\mathbf{k}A\downarrow}^{\dagger},c_{\mathbf{k}B\uparrow}^{\dagger},c_{\mathbf{k}B\downarrow}^{\dagger})$,

and $\mathbf{r}_{1}=\frac{\sqrt{3}}{2}\hat{x}-\frac{1}{2}\hat{y},\mathbf{r}_{2}=\hat{y},\mathbf{r}_{3}=-\frac{\sqrt{3}}{2}\hat{x}-\frac{1}{2}\hat{y}$. $a$ is the lattice constant. Near $\pm\mathbf{K}=\pm\frac{4\pi}{3a}\hat{x}$ we have

The Dirac Hamiltonian at each valley $\pm\mathbf{K}\equiv\eta_{v}\mathbf{K}$ is therefore

where $\bm{\tau}$ is the Pauli matrix vector in the sublattice space. One can see that the magnetic-charge contribution only changes the magnitude of the Rashba terms in each valley. Since the Rashba terms do not gap the Dirac Hamiltonian by themselves [52], the magnetic charge contribution cannot lead to topological phase transitions. Separately, since $t_{M}$ effectively changes the relative strength of the Rashba spin-orbit coupling at the two valleys, when the Fermi energy is not within the gap, the Berry curvature at the two valleys will not cancel out, leading to the finite anomalous Hall conductivity.

The minimal model can be connected with the magnetic-charge-ordered state of the kagome spin ice [37, 38, 45] by the duality between honeycomb and kagome lattices. The background magnetic field connecting neighboring magnetic charges can be regarded, as a first approximation, as the homogenized effect of the fluctuating magnetic dipole moments in the charge-ordered state of the kagome spin ice on itinerant electrons. In comparison with the models studied in [53, 54] where the local spins on the kagome lattice is noncoplanar, the present model has a vanishing net magnetization. More importantly, the essential symmetry breaking in the non-magnetic part of model (IV.1) is already present in the pyrochlore iridate Pr₂Ir₂O₇ [9], in which an AHE in the absence of long-range dipolar order and of net magnetization has been observed [55]. Although the ground state of Pr₂Ir₂O₇ may be elusive and the direct measurement [56] of scalar spin chirality by scattering techniques is challenging, it is possible to alternatively measure the EC which, if nonzero, can help to solve the puzzle of the zero-field AHE in [55].

In this subsection we predict an extrinsic contribution to the AHE by magnetic charge through skew scattering. Again motivated by the expression of $\bm{\chi}_{1}$, we consider the following model of 2D Dirac electrons with Rashba-type spin-momentum locking scattered by a magnetic charge whose magnetic field is truncated at finite radius $R$:

where $\Theta(x)$ is the step function. $H_{\Delta}$ represents the Zeeman coupling between the electron spin and the magnetic field $\mathbf{h}(\mathbf{r})=-\mathbf{m}(\mathbf{r})=\Delta\hat{r}/(2\pi r)$ generated by a magnetic charge located at the origin within a radius $R$. Here we consider the analytically simpler case of $\alpha\equiv\Delta/(2\pi\hbar\lambda)=1$ and relegate the more general solution to Appendix B. Assuming a positive chemical potential, the solution for $r<R$ with energy $E=\hbar\lambda k_{0}>0$ is

while that for $r>R$ with an incident plane wave traveling along $\hat{x}$ is

where $H_{n}=J_{n}+\imath Y_{n}$ is the Hankel function of the first kind, and $Y_{n}$ is the Bessel function of the second kind. The appearance of the Hankel function is because $\Psi_{>}$ does not include the origin where $Y_{n}$ diverges. Also the Hankel function of the first kind represents outgoing waves [5]. Solving $b_{n}$ from the boundary condition $\Psi_{<}|_{r=R}=\Psi_{>}|_{r=R}$ and taking the large-distance asymptotic form of $\Psi_{\rm scatt}$, we obtain the scattering cross section

When $k_{0}R\ll 1$ one can consider up to $p$-wave contributions to $\sigma(\theta)$. The Hall angle due to scattering by the magnetic charge only can be calculated as [57]

The result in the last line above turns out to be a good approximation even at $k_{0}R\sim 1$, for which $\tan\theta_{\rm H}\approx 0.39$. We also found that when $\alpha=-1$, $\tan\theta_{\rm H}\approx-\frac{\pi}{8}(k_{0}R)^{2}$, confirming the time-reversal-odd property of the AHE.

Experimental detection of such an effect may be performed using topological insulator surface states [58] or 2D electron gas with large Rashba spin-orbit coupling that are proximate-coupled to magnetic textures [59, 60] having a nonzero 2D magnetic charge density. With $\hbar\lambda\sim 10^{1}\,{\rm eV\cdot\AA}$ and $R\sim 1/k_{0}\sim 100$ Å, the Zeeman coupling $\Delta/(2\pi R)\sim 0.1$ eV is reasonable to achieve experimentally.