Microscopic Theory of Spin Toroidization in Periodic Crystals

TL;DR

This paper develops a gauge-invariant semiclassical theory for spin toroidization in periodic crystals, revealing its components, relation to magnetoelectric effects, and demonstrating its bulk nature through a tight-binding model.

Contribution

It introduces a novel gauge-invariant expression for spin toroidization, linking it to magnetoelectric polarizability and establishing its bulk property in crystals.

Findings

Derived a gauge-invariant formula for spin toroidization.

Linked spin toroidization to antisymmetric magnetoelectric polarizability.

Confirmed spin toroidization as a genuine bulk quantity in a model.

Abstract

Using the semiclassical theory of electron dynamics, we derive a gauge-invariant expression for the spin toroidization in a periodical crystal. We show that the spin toroidization is comprised of two contributions: one is due to the configuration of a classical spin array, while the other comes from the coordinate shift of the electron as spin carrier in response to the inhomogeneous magnetic field. We then establish a direct and elengant relation between our spin toroidization and the antisymmetric magnetoelectric polarizability in insulators. Finally, we demonstrate our spin toroidization in a tight-binding model and show that it is a genuine bulk quantity.

Equations

$$\mathcal{T}=\frac{g{\mu }_B}{2\hslash V}\sum _i{\mathbf{r}}_i\times {\mathbf{s}}_i,$$

$$\begin{array}{rl}F(\mathbf{r})& =F_0(\mathbf{r})-\mathbf{M}\cdot \mathbf{B}(\mathbf{r})-{\mathcal{Q}}_{ij}{\partial }_i{B}_j(\mathbf{r})+\cdots ,\end{array}$$

$${\mathcal{T}}_k=\frac{1}{2}{ϵ}_{ijk}{\mathcal{Q}}_{ij},$$

$$\bm{\mathcal{T}}(\bm{r})=-\lim_{\bm{B}(\bm{r})\to 0}\frac{\partial F(\bm{r})}{\partial(\bm{\nabla}\times\bm{B})}\biggr{|}_{\bm{B}(\bm{r})}\;.$$

$${\hat{H}}_{\text{F}}=\hat{H}(i\hslash {\partial }_{\mathbf{r}},\mathbf{r})-\frac{g{\mu }_B}{\hslash }\mathbf{B}(\mathbf{r})\cdot \hat{\mathbf{s}}.$$

$$\mathcal{D}({\mathbf{r}}_c,{\mathbf{k}}_c)=1+\mathrm{Tr}({\Omega }_{\mathbf{k},\mathbf{r}}),$$

$$({\Omega }_{\mathbf{k},\mathbf{r}}{)}_{ij}=-2\mathrm{Im}⟨{\partial }_{k_{ci}}{\tilde{u}}_0|{\partial }_{r_{cj}}{\tilde{u}}_0⟩$$

$${\epsilon }_0^{\prime }={\tilde{\epsilon }}_0+\mathrm{Im}⟨{\partial }_{k_{ci}}{\tilde{u}}_0|({\tilde{\epsilon }}_0-{\hat{H}}_c)|{\partial }_{r_{ci}}{\tilde{u}}_0⟩.$$

$$\delta F=-{\int }^{\mu }\frac{d{\mathbf{k}}_c}{(2\pi {)}^3}\mathrm{Im}⟨{\partial }_{k_{ci}}{\tilde{u}}_0|({\tilde{\epsilon }}_0+{\hat{H}}_c-2\mu )|{\partial }_{r_{ci}}{\tilde{u}}_0⟩.$$

$$\displaystyle\bm{\mathcal{T}}=\frac{1}{2}\int^{\mu}\frac{d\bm{k}}{(2\pi)^{3}}\operatorname{Im}\langle\partial_{\bm{k}}\tilde{u}_{0}|\times(\tilde{\varepsilon}_{0}+\hat{H}_{c}-2\mu)|\partial_{\bm{B}}\tilde{u}_{0}\rangle\Bigr{|}_{\bm{B}\to 0}\;.$$

$$|{\tilde{u}}_0⟩=|u_0⟩-\frac{g{\mu }_B}{\hslash }\sum _{n\ne 0}\frac{\mathbf{B}\cdot {\mathbf{s}}_{n0}}{{\epsilon }_0-{\epsilon }_n}|u_n⟩.$$

$$\mathcal{T}=-\frac{g{\mu }_B}{2}\sum _{n\ne 0}{\int }^{\mu }\frac{d\mathbf{k}}{(2\pi {)}^3}({\epsilon }_0+{\epsilon }_n-2\mu )\frac{\mathrm{Im}({\mathbf{v}}_{0n}\times {\mathbf{s}}_{n0})}{({\epsilon }_0-{\epsilon }_n{)}^2},$$

$${\mathcal{Q}}_{ij}=-g{\mu }_B\sum _{n\ne 0}{\int }^{\mu }\frac{d\mathbf{k}}{(2\pi {)}^3}({\epsilon }_0+{\epsilon }_n-2\mu )\frac{\mathrm{Im}[(v_i{)}_{0n}(s_j{)}_{n0}]}{({\epsilon }_0-{\epsilon }_n{)}^2}.$$

$$\begin{split}\bm{\mathcal{T}}&=\frac{1}{V_{\text{cell}}}\operatorname{Re}\langle w_{0}(\bm{B})|\bm{r}(\hat{H}_{c}-\mu)\times\bm{\partial}_{\bm{B}}|w_{0}(\bm{B})\rangle\Bigr{|}_{\bm{B}\to 0}-\frac{g\mu_{B}}{2\hbar V_{\text{cell}}}\langle w_{0}(\bm{B})|\bm{r}\times\hat{\bm{s}}|w_{0}(\bm{B})\rangle\Bigr{|}_{\bm{B}\to 0}\\ &\quad-\frac{g\mu_{B}}{2\hbar V_{\text{cell}}}\sum_{\bm{R}}\langle w_{0}(\bm{B})|\bm{r}|w_{0}(\bm{R},\bm{B})\rangle\Bigr{|}_{\bm{B}\to 0}\times\langle w_{0}(\bm{R},\bm{B})|\hat{\bm{s}}|w_{0}(\bm{B})\rangle\Bigr{|}_{\bm{B}\to 0}\;,\end{split}$$

$$\bm{\mathcal{T}}_{1}=\frac{g\mu_{B}}{2\hbar V_{\text{cell}}}\langle w_{0}(\bm{B})|\bm{r}\times\hat{\bm{s}}|w_{0}(\bm{B})\rangle\Bigr{|}_{\bm{B}\to 0}\;.$$

$$\bm{\mathcal{T}}_{2}=-\frac{1}{2V_{\text{cell}}}(\varepsilon_{0}-\mu)(\partial_{\bm{B}}\times\bar{\bm{r}})\Bigr{|}_{\bm{B}\to 0}\;,$$

$$\Delta P_i=e\int \frac{d\mathbf{k}d{B}_j}{(2\pi {)}^3}\mathrm{Im}⟨{\partial }_{k_i}{\tilde{u}}_0|{\partial }_{B_j}{\tilde{u}}_0⟩.$$

$$\alpha_{ij}=\frac{\partial P_{i}}{\partial B_{j}}\biggr{|}_{\bm{B}\to 0}=e\int\frac{d\bm{k}}{(2\pi)^{3}}\operatorname{Im}\langle\partial_{k_{i}}\tilde{u}_{0}|\partial_{B_{j}}\tilde{u}_{0}\rangle\Bigr{|}_{\bm{B}\to 0}\;.$$

$$e\frac{\partial {\mathcal{T}}_k}{\partial \mu }=-\frac{1}{2}{ϵ}_{ijk}{\alpha }_{ij},$$

$$\frac{\partial {\mathcal{T}}_i}{\partial \mu }=\frac{\partial \rho }{\partial (\nabla \times \mathbf{B}{)}_i}.$$

$$e\rho =\nabla \cdot \mathbf{P}={\partial }_i({\alpha }_{ij}{B}_j)=\frac{1}{2}{ϵ}_{ijk}{\alpha }_{ij}{h}_k.$$

$$e\frac{\partial {\mathcal{Q}}_{ij}}{\partial \mu }=-{\alpha }_{ij}.$$

$${\hat{H}}_{\text{TB}}=-\Delta \sum _i{\mathbf{n}}_i\cdot {\mathbf{\sigma }}_{\alpha \beta }{c}_{i\alpha }^{†}{c}_{i\beta }+\sum _{⟨i,j⟩}{t}_{ij}{c}_{i\alpha }^{†}{c}_{j\alpha },$$

$$\mathcal{T}=-\sum _{m,n}\frac{E_m+E_n-2\mu }{E_m-E_n}(⟨{\psi }_m|\mathbf{r}|{\psi }_n⟩\times ⟨{\psi }_n|\mathbf{\sigma }|{\psi }_m⟩+\text{c.c.}),$$

$$\mathcal{T}={\mathcal{T}}_{\text{bulk}}+\frac{a_1}{N}+\frac{a_2}{N^2}$$

$$\mathbf{B}=(h/2q)(\sin (qx)\hat{y}-\sin (qy)\hat{x}),$$

$$(h/2)(\cos (qx)+\cos (qy))\hat{z},$$

$$\delta F(\mathbf{r})=\sum _{n\mathbf{k}}(\delta f_{n\mathbf{k}}){\psi }_{n\mathbf{k}}^{\ast }{\hat{F}}_0{\psi }_{n\mathbf{k}}+f_{n\mathbf{k}}{\psi }_{n\mathbf{k}}^{\ast }\mathbf{B}\cdot \hat{\mathbf{s}}{\psi }_{n\mathbf{k}}+f_{n\mathbf{k}}(\delta {\psi }_{n\mathbf{k}}^{\ast }{\hat{F}}_0{\psi }_{n\mathbf{k}}+{\psi }_{n\mathbf{k}}^{\ast }{\hat{F}}_0\delta {\psi }_{n\mathbf{k}}).$$

$${\mathcal{T}}_z=-\frac{2}{Vh}\int dx\delta F(\mathbf{r})(\cos (qx)+\cos (qy)).$$

$$\delta\psi_{n\bm{k}}=-{hg\mu_{B}\over 4i\hbar q}\Bigl{[}\sum_{n^{\prime}}{e^{i(\bm{k}+\bm{q})\cdot\bm{r}}|n^{\prime}\bm{k}+\bm{q}\rangle\langle n^{\prime}\bm{k}+\bm{q}|\hat{s}_{y}|n\bm{k}\rangle\over\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}+\bm{q}}}-(\bm{q}\rightarrow-\bm{q})\Bigr{]}\;,$$

Full text

There has been continuous interest in toroidal moments in crystals Dubovik1990 ; Gorbatsevich1994 ; Ederer2007 ; Spaldin2008 , mainly due to their intriguing role in various magnetoelectric effects Fiebig2005 ; Arima2005 ; Sawada2005 ; Aken2007 ; Zimmermann2014 ; Hayami2014 ; Fiebig2016 . A toroidal moment is generally associated with a vortex-like structure of magnetic moments. Its spontaneous ordering characterizes a ferrotoroidal state that may exhibit a non-vanishing magnetoelectric effect. The density of the toroidal moment, called toroidization, also constitutes an essential building block in the free energy expansion in inhomogeneous fields. However, despite its importance, a microscopic theory of the toroidization based on quantum mechanical wave functions is still missing.

In crystals the toroidization can arise from two sources, the orbital and spin moments. Here we will focus only on the contribution from spins foot1 . By treating the spins as classical vectors, it has been proposed that the spin toroidization can be written as Ederer2007 ; Spaldin2008

[TABLE]

where $g$ is the gyromagnetic factor, $\mu_{B}$ is the Bohr magneton, $\bm{r}_{i}$ and $\bm{s}_{i}$ are the position and spin of each lattice site, and $V$ is the volume of the sample. There is also a recent attempt to obtain a microscopic theory of the spin toroidization by treating $\bm{r}$ and $\bm{s}$ as operators and directly evaluating the expectation of Eq. (1) using Wannier functions Thole2016 . However, the resulting expression is not gauge-invariant.

In this work, we develop a quantum theory of spin toroidization in crystals. Using the semiclassical theory of electron dynamics Sundaram1999 ; Xiao2010 , we obtain a gauge-invariant expression for the spin toroidization in terms of bulk Bloch functions, which is amenable to implementation in first-principles codes. By considering the molecular insulator limit, we find that the contributions to the spin toroidization consists of two parts with clear physical interpretations: one is due to the configuration of a classical spin array, similar to Eq. (1), while the other comes from the coordinate shift of the electron as spin carrier in response to the inhomogeneous magnetic field.

Using our theory, we are able to establish a direct and elegant relation between the spin toroidization and the antisymmetric magnetoelectric polarizability in the case of insulators [see Eq. (19)]. Such a relation is dictated by general thermodynamic principles. Finally, using a tight-binding toy model, we show that our spin toroidization is a genuine bulk quantity.

I General Formalism

In this section, we first introduce the definition of the toroidization as a response function of the free energy to the derivative of the magnetic field. Then we use the semiclassical theory to derive the spin toroidization.

I.1 Toroidization as a response function

Our starting point is the free energy density $F(\bm{r})$ in an inhomogeneous magnetic field $\bm{B}(\bm{r})$. Suppose that $\bm{B}(\bm{r})$ is small and varies slowly in space. At a given point $\bm{r}$, we can perform a gradient expansion of $F(\bm{r})$ up to first order with respect to the derivatives of the magnetic field

[TABLE]

where $F_{0}(\bm{r})$ is the free energy density at $\bm{B}=0$ and $\bm{M}$ is the magnetization. Here and hereafter the Einstein summation convention is implied for repeated indices. The quantity $\mathcal{Q}_{ij}$ is the magnetic quadrupole moment density. The toroidization $\bm{\mathcal{T}}$ is the antisymmetric part of $\mathcal{Q}_{ij}$:

[TABLE]

where $\epsilon_{ijk}$ is the total antisymmetric tensor.

Based on Eq. (2), we can obtain a linear-response expression for $\bm{\mathcal{T}}(\bm{r})$ by treating $\bm{B}(\bm{r})$ and $\bm{\nabla}\times\bm{B}(\bm{r})$ as independent variables, arriving at

[TABLE]

Here the subscript $\bm{B}(\bm{r})$ in taking the derivative with respect to $\bm{\nabla}\times\bm{B}$ means that the magnetic field at the point $\bm{r}$ has to be kept fixed as $\bm{\nabla}\times\bm{B}(\bm{r})$ is varied.

I.2 Semiclassical theory of the spin toroidization

With the above definition of the toroidization, we now formulate its microscopic theory. We will focus only on the spin toroidization. Therefore, we take $\bm{B}(\bm{r})$ as the Zeeman field, which couples to the spin operator $\hat{\bm{s}}$. Then the full Hamiltonian can be written as

[TABLE]

For definiteness we can consider the context to be that of a spinor implementation of density functional theory with spin-orbit interactions included. The first term, $\hat{H}(i\hbar\partial_{\bm{r}},\bm{r})$, describes a perfect crystal in the absence of a Zeeman field, while the second term is inhomogeneous and breaks the translational symmetry, making it difficult to diagonalize the Hamiltonian (5) analytically. Here we take a different route by using the semiclassical theory of electron dynamics Sundaram1999 ; Xiao2010 , which is designed to study Bloch electrons subject to perturbations varying slowly in space.

In the spirit of the semiclassical theory, each Bloch electron responds to the external Zeeman field in the form of a wave packet, which has a specified center of mass position $\bm{r}_{c}$ and momentum $\bm{k}_{c}$. To construct the wave packet, we make a local approximation and assume that the system can be described by a set of local Hamiltonians $\hat{H}_{c}[\bm{B}(\bm{r}_{c})]=\hat{H}(i\hbar\partial_{\bm{r}},\bm{r})-(g\mu_{B}/\hbar)\bm{B}(\bm{r}_{c})\cdot\bm{s}$. Since $\hat{H}_{c}$ respects the lattice translational symmetry, its eigenstate has the form of a Bloch function $e^{i\bm{k}\cdot\bm{r}}|\tilde{u}_{n}(\bm{k},\bm{B}(\bm{r}_{c}))\rangle$ with the eigenenergy $\tilde{\varepsilon}(\bm{k}_{c},\bm{B}(\bm{r}_{c}))$, where $n$ is the band index. In the limit $\bm{B}(\bm{r}_{c})\to 0$, $|\tilde{u}_{n}(\bm{k},\bm{B}(\bm{r}_{c}))\rangle$ reduces to $|u_{n}(\bm{k})\rangle$, the periodic part of the Bloch function of $\hat{H}$, and $\tilde{\varepsilon}_{n}(\bm{k}_{c},\bm{B}(\bm{r}_{c}))$ reduces to $\varepsilon_{n}(\bm{k}_{c})$, the eigenenergy of the unperturbed Hamiltonian $\hat{H}$. For simplicity, we hereafter drop the argument of $|\tilde{u}_{n}\rangle$, $|u_{n}\rangle$, $\tilde{\varepsilon}_{n}$ and $\varepsilon_{n}$. For illustrative purposes, we consider a single band with index [math], and the wave packet is thus the superposition of $e^{i\bm{k}\cdot\bm{r}}|\tilde{u}_{0}\rangle$.

The wave-packet dynamics can be properly formulated as a set of semiclassical equations of motion in the phase space spanned by $\bm{r}_{c}$ and $\bm{k}_{c}$ Sundaram1999 ; Xiao2010 . The spatial inhomogeneity of $\bm{B}(\bm{r})$ introduces two essential ingredients for the purpose of evaluating the spin toroidization in Eq. (4). First, the phase space density of states $\mathcal{D}$ is modified. It has the form Xiao2005 (see also Sec. VI.B of Ref. Xiao2010, )

[TABLE]

where

[TABLE]

is the mixed Berry curvature between the real and momentum space. This modified density of states has been applied to derive the polarization in inhomogeneous crystals Xiao2009 . Secondly, the band energy $\tilde{\varepsilon}$ is also affected by the spatial inhomogeneity Sundaram1999

[TABLE]

With the above two ingredients we are ready to evaluate the free energy density $F$. For simplicity we set $T=0$. The free energy density is given by $F=\int\frac{d\bm{k}_{c}}{(2\pi)^{3}}\mathcal{D}(\bm{r}_{c},\bm{k}_{c})(\varepsilon_{0}^{\prime}-\mu)\Theta(\mu-\varepsilon_{0}^{\prime})$, where $\Theta$ is the Heaviside function. At first order with respect to the derivative of $\bm{B}$, the correction to the free energy density is

[TABLE]

Here the upper limit $\mu$ means that the integration is taken up to $\varepsilon_{0}=\mu$.

The toroidization defined in Eq. (4) can be obtained from the above free-energy correction. Since $|\tilde{u}_{0}\rangle$ depends on $\bm{r}_{c}$ through $\bm{B}$, we make the substitution $\partial_{r_{ci}}|\tilde{u}_{0}\rangle=\partial_{r_{ci}}B_{\ell}\partial_{B_{\ell}}|\tilde{u}_{0}\rangle$. We then collect terms involving the antisymmetric part of $\partial_{r_{ci}}B_{\ell}$ and take the derivative as in Eq. (4). The final expression is

[TABLE]

Here and hereafter we drop the subscript $c$ of $\bm{k}_{c}$. Note that Eq. (10) can be straightforwardly generalized to the multiband case by summing over all occupied states [see Eq. (35)], but we continue to focus on the single-band case here.

The structural similarity between Eq. (10) and the orbital magnetization formula Xiao2005 ; Thonhauser2005 ; Ceresoli2006 ; Shi2007 is striking. In fact, by making the substitution $\bm{\partial}_{\bm{B}}\to\bm{\partial}_{\bm{k}}$, Eq. (10) exactly coincides with the expression of the orbital magnetization. This similarity has its root in the nature of spin toroidization and orbital magnetization: they both measure the moment of some observable, which is spin for spin toroidization and velocity for orbital magnetization.

Equation (10) can be cast in a form involving only unperturbed Bloch states $|u_{0}\rangle$ instead of $|\tilde{u}_{0}\rangle$. Using the perturbation theory, up to the first order in the Zeeman field we have

[TABLE]

Then Eq. (10) can be rewritten as

[TABLE]

where $\bm{v}_{0n}=\langle u_{0}|\hat{\bm{v}}|u_{n}\rangle$ and $\bm{s}_{n0}=\langle u_{n}|\hat{\bm{s}}|u_{0}\rangle$ are the interband elements of the velocity and spin operators, respectively. Both Eq. (10) and (12) are amenable to implementation in a first-principles calculation. To further check the validity of our result, we have also carried out a linear response calculation (see Appendix A for details), similar to the derivation of the orbital magnetization in Ref. Shi2007, , and obtained the same result.

We comment that in general the spin magnetic quadrupole moment density can obtained in a similar way. The result reads

[TABLE]

One can easily check that $\bm{\mathcal{T}}$ and $\mathcal{Q}_{ij}$ satisfy Eq. (3).

It is clear that our expression (12) for the spin toroidization is gauge-invariant since it does not change if an arbitrary phase factor is applied to $|u_{n}\rangle$. As a consequence, the spin toroidization does not have any quantum of uncertainty, and it always vanishes for a system with either time-reversal or inversion symmetry. This is in sharp contrast to both the electric polarization Kingsmith1993 ; Resta1994 and the previous theory of the spin toroidization Ederer2007 ; Thole2016 .

It is also worth mentioning that our toroidization cannot be used to predict a surface magnetization density, unlike the electric polarization, which has a definitive relation to the surface charge density Vanderbilt1993 . This difference can be traced to the fact that charge is conserved but spin is not.

II Interpretation of spin toroidization

In this section, we explore the physical meaning of the spin toroidization in the Wannier representation and discuss the difference between our result and the classical definition of spin toroidization in Eq. (1). Finally, we show that the spin toroidization can be directly related to the spin magnetoelectric polarizability.

II.1 Molecular Insulator Limit

To shed light on the physical meaning of the spin toroidization in Eq. (10), we rewrite it for an insulator using the Wannier function representation. We label the Wannier function defined from the local Hamiltonian $\hat{H}_{c}$ by $|w_{0}(\bm{R},\bm{B})\rangle$, with [math] being the band index and $\bm{R}$ being the lattice site. In this representation Eq. (10) becomes (see Appendix B for details)

[TABLE]

where $|w_{0}(\bm{B})\rangle=|w_{0}(\bm{R},\bm{B})\rangle$ with $\bm{R}=0$, and $V_{\text{cell}}$ is the unit cell volume.

The meaning of Eq. (14) can be clarified further by taking the molecular insulator limit. Since the spin toroidal moment arises from a vortex-like arrangement of spins, there must be multiple atoms in a unit cell, which we call a molecule. The molecular insulator limit is then taken by letting the distance between neighboring molecules go to infinity while the relative structure of each molecule is unchanged. In this limit, $|w_{0}(\bm{R},\bm{B})\rangle$ is just the energy eigenstate of the molecule, translated to sit in cell $\bm{R}$. We will further assume that the system respects the combined time reversal and inversion symmetry such that $\langle w_{0}|\hat{\bm{s}}|w_{0}\rangle$ vanishes.

In the molecular insulator limit Eq. (14) consists of two parts (see Appendix C for details). The first part is

[TABLE]

It is clear that this term is due to the configuration of an array of classical spins, similar to the equation appearing as Eq. (1) in the classical picture. The second part, coming from the modified density of states, is

[TABLE]

where $\varepsilon_{0}$ refers to the molecular electronic energy levels and $\bar{\bm{r}}=\langle w_{0}(\bm{B})|\bm{r}|w_{0}(\bm{B})\rangle$ is the electron position under the external Zeeman field. Here $\varepsilon_{0}-\mu$ is the free energy for state [math].

The $\bm{\mathcal{T}}_{2}$ term can be intuitively understood as follows. In the spirit of the molecular insulator limit, if each molecule is simply a cluster of classical spins, under an inhomogeneous magnetic field the spins on each site can rotate but cannot move. However, the spins are carried by electrons, and the inhomogeneous Zeeman field will exert a spin force on the electron. Therefore, the electron will shift to a new equilibrium position due to the balance between the spin force and the restoring force that binds electrons to ions. The corresponding energy change gives rise to $\bm{\mathcal{T}}_{2}$. In a semiclassical picture, $\partial_{\bm{B}}\times\bar{\bm{r}}$ counts the change of the number of electronic states within a volume element located at $\bm{r}$.

Equation (16) also provides a strong hint connecting the toroidization with the magnetoelectric polarizability. Taking the derivative with respect to $\mu$ in Eq. (16) yields $(1/2V_{\text{cell}})\bm{\partial}_{\bm{B}}\times\bm{\bar{}}{\bm{r}}$. Since $\bar{\bm{r}}$ is proportional to the electric polarization, its derivative with respect to the Zeeman field $\bm{B}$ is exactly the magnetoelectric polarizability. We show below that this is a general relation born out from our theory.

II.2 Connection to magnetoelectric polarizability

It is well known that the toroidization and the antisymmetric part of the magnetoelectric polarizability transform in the same way under symmetry operations Spaldin2008 . However, an explicit relation between these two quantities has not previously been identified. Here we show that for an insulator, the spin toroidization admits a direct and elegant connection to the spin magnetoelectric polarizability.

According to the modern theory of polarization Kingsmith1993 ; Resta1994 , as we vary the $j$-th component of the Zeeman field, the change of the polarization is given by

[TABLE]

Therefore, the magnetoelectric polarizability has the form Essin2010

[TABLE]

On the other hand, note that for an insulator the Fermi-surface contribution to $\bm{\mathcal{T}}$ vanishes. If we take the derivative of Eq. (10) with respect to $\mu$, we find the desired connection

[TABLE]

where $\epsilon_{ijk}$ is the Levi-Civita symbol.

There is a heuristic derivation of the relation (19). Equation (2) suggests that the differential form of the free energy is $dF=-\bm{\mathcal{T}}\cdot d(\bm{\nabla}\times\bm{B})-\rho d\mu$, where $\rho$ is the particle density. We can then obtain via the Maxwell relation

[TABLE]

For a given point $\bm{r}$, we write $\bm{B}=(1/2)\bm{h}\times\bm{r}$ such that $\bm{\nabla}\times\bm{B}=\bm{h}$ and $\bm{B}$ vanishes exactly at $\bm{r}$. This choice ensures that $\partial_{i}B_{j}$ only has the antisymmetric component. On the other hand, the application of an inhomogeneous Zeeman field $\bm{B}(\bm{r})$ will induce an inhomogeneous polarization $\bm{P}(\bm{r})$, which in turn leads to a charge density change, i.e.,

[TABLE]

Combining Eq. (20) and (21) then yields Eq. (19).

In fact, similar relation exists between the spin magnetic quadrupole moment density $\mathcal{Q}_{ij}$ and each component of the spin magnetoelectric polarizability. Using similar derivations, one can find that

[TABLE]

This relation is also implied from Eq. (19) based on the definition (3).

We note that the above argument is thermodynamic in nature and does not depend on microscopic details. Therefore, it is valid for any physical systems, and also for orbital toroidization, where the Zeeman field $\bm{B}$ is replaced by a magnetic field.

III Spin toroidization as a bulk quantity

For a response function to reflect the bulk properties of the sample, it is essential that the response function has a valid thermodynamics limit, i.e., it has a well-defined limit as the sample size grows to infinity. This is well illustrated in the modern theory of electric polarization and orbital magnetization Bianco2013 ; Marrazzo2016 ; Ceresoli2006 ; Thonhauser2005 ; Kingsmith1993 ; Vanderbilt1993 ; Shi2007 ; Resta1994 .

However, the spin toroidization in Eq. (1) in the classical picture is not a bulk quantity. To show this, we consider the bulk sample in Fig. 1, which has uniform but opposite surface magnetization on the left and right surfaces. We choose the origin to be the center of the sample and label the left and right surface magnetization, the surface area, and the distance between two surfaces by $M_{\text{s}}$, $-M_{\text{s}}$, $S$, and $L$, respectively. Then from Eq. (1), the contribution from the surface reads $g\mu_{B}/2[(L/2)SM_{\text{s}}/V+(L/2)SM_{\text{s}}/V]=g\mu_{B}M_{\text{s}}/2$, which obviously does not vanish as $V\rightarrow\infty$. Moreover, this contribution can point to any direction as the direction of the surface magnetization varies.

In comparison, our spin toroidization in Eq. (10) is indeed a genuine bulk quantity. In the following, we demonstrate this point by studying a toy model.

III.1 Toy model

We consider the model Hamiltonian

[TABLE]

where $\Delta$ is the local exchange field, $\bm{n}_{i}$ is the exchange field direction, $\alpha$ and $\beta$ label the spin components, and $t_{ij}$ is the spin-independent nearest neighbor hopping strength alternating between $t_{1}$ and $t_{2}$ as shown in Fig. 2a.

Since $\bm{\mathcal{T}}$ transforms as a vector, it is useful to first analyze the symmetry of this system. The system has a 4-fold rotational symmetry about the vertical axis. Therefore, the toroidization cannot have any in-plane component. Moreover, if $t_{1}=t_{2}$, the system also respects the combined symmetry of the mirror operation $\sigma_{h}$ (i.e., $z\to-z$) followed by a translation across the diagonal direction. This requires that the out-of-plane component of the toroidization vanishes. Both results have been verified in our numerical calculations. We thus focus on the case $t_{1}\neq t_{2}$ for which an out-of-plane toroidization is expected.

Figure 2b shows the toroidization calculated from Eq. (10) as a function of the chemical potential. The system has four bands separated by three global band gaps, and each band is doubly degenerate. The curve is symmetric with respect to $\mu=0$ because of the particle-hole symmetry of our model. When the chemical potential falls inside the lowest and highest gap region, the toroidization varies linearly, and we confirm that its slope coincides with the magnetoelectric polarizability, consistent with our Eq. (19). When the chemical potential falls inside the middle gap region, the toroidzation is a nonzero constant while the magnetoelectric polarizability vanishes, in accordance with the particle-hole symmetry. Our toy model thus represents an interesting scenario with a vanishing magnetoelectric polarizability but a finite spin toroidization when the Fermi energy is in the middle gap.

III.2 Irrelevance of the boundary

Now we consider a finite sample with the model Hamiltonian in Eq. (23) using open boundary conditions. We label the $n$-th eigenstate by $E_{n}$ and the corresponding wave function by $|\psi_{n}\rangle$. Eq. (12) reduces to

[TABLE]

where $m$ labels occupied states and $n$ labels unoccupied states, and $\bm{\mathcal{T}}$ is in units of $g\mu_{B}/4S$ with $S$ being the area of the sample.

The finite-sample value of the spin toroidization is calculated based on Eq. (24) and plotted in Fig. 3. The bulk value is $-0.126$ in units of $g\mu_{B}/4a$ and is displayed as the straight black line in Fig. 3. The red dots and blue diamonds are the finite-sample results with an even and odd number of lattice sites along the sample edge respectively. They fall on different curves due to the different surface terminations as shown in the inset. For large $N$, the finite-sample results should asymptotically satisfy Ceresoli2006

[TABLE]

where $\mathcal{T}_{\text{bulk}}$ is the bulk value. The second and third terms are due to the edge and corner contributions respectively. We find that the red and blue dots indeed fit Eq. (25) very well. This clearly demonstrates that our spin toroidization in Eq. (10) and (12) is a genuine bulk property.

To further demonstrate the irrelevance of the boundary condition, we rotate the spin direction on the surface of our finite sample as shown in the inset of Fig. 4. This type of perturbation obviously contributes to Eq. (1) in the thermodynamic limit. In Fig. 4, we plot the calculated spin toroidization based on Eq. (12) for this type of finite samples. We find that although the toroidizations of finite samples are changed by the surface perturbation, they still follow the asymptotic rule in Eq. (25) very well and converge to the same bulk value.

IV Summary

In this work, we develop a quantum theory of spin toroidization in crystals. Using the semiclassical theory of electron dynamics, we obtain a gauge-invariant expression for the spin toroidization in terms of bulk Bloch functions that is amenable to implementation in first-principles codes. Using our theory, we are able to establish a direct and elegant relation between the spin toroidization and the antisymmetric magnetoelectric polarizability in the case of insulators, as dictated by general thermodynamic principles. Finally, using a tight-binding toy model, we show that our spin toroidization is a genuine bulk quantity.

Acknowledgements.

This work is supported by the Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division, Grant No. DE-SC0012509 (Y.G. and D.X.). D.V. is supported by the National Science Foundation, Grant No. DMR-1408838.

Appendix A Spin toroidization in linear response theory

In this section we derive the spin toroidization using linear response theory, similar to the calculation of the orbital magnetization in Ref. Shi2007, . Without loss of generality, we only calculate the $z$-th component of the spin toroidization.

As given by Eq. (4), the spin toroidization is the response of the free energy density to an external Zeeman field. Let us consider a Zeeman field of the following form

[TABLE]

where $h$ is small. This Zeeman field has a curl

[TABLE]

which reduces to $h$ in the limit $q\rightarrow 0$. Note that the symmetric part of the derivative of the Zeeman field (26) is $\partial_{x}B_{y}+\partial_{y}B_{x}=(h/2)(\cos(qx)-\cos(qy))$, which vanishes in the limit $q\rightarrow 0$. As a result, the response of the free energy density to this Zeeman field in the limit $q\rightarrow 0$ is purely due to its curl.

At zero temperature, the free energy density reads $\hat{F}=\hat{H}-\mu\hat{N}$. With the above Zeeman field, the change in $F$ can be divided into four parts:

[TABLE]

Here $\psi_{n\bm{k}}=e^{i\bm{k}\cdot\bm{r}}|u_{n}(\bm{k})\rangle$ is the Bloch function of the unperturbed Hamitonian $\hat{H}$ with $\varepsilon_{n\bm{k}}$ being the corresponding eigenenergy, $f_{n\bm{k}}$ is the Fermi function, and $\hat{F}_{0}$ is the unperturbed part of the free energy density. The spin toroidization can be obtained from the appropriate Fourier component of $\delta F(\bm{r})$:

[TABLE]

We first calculate the contribution from the $y$-component of the Zeeman field. The perturbation to the wave function is

[TABLE]

where $\bm{q}=q\hat{x}$, and $|n\bm{k}\rangle$ is short for $|u_{n}(\bm{k})\rangle$. Inserting this expression into Eq. (28). We can see that the first two terms cancel each other. The last two terms read

[TABLE]

Now we take the limit $q\rightarrow 0$ in the above expression. Terms in Eq. (31) with $n\neq n^{\prime}$ reads

[TABLE]

where $\bm{A}_{nn^{\prime}}=\langle n\bm{k}|i\bm{\partial}_{\bm{k}}|n^{\prime}\bm{k}\rangle$ is the interband Berry connection and $\bm{s}_{n^{\prime}n}=\langle n^{\prime}\bm{k}|\hat{\bm{s}}|n\bm{k}\rangle$ is the interband element of the spin operator. Terms in Eq. (31) with $n=n^{\prime}$ reads

[TABLE]

Note that to obtain the last equality in the above equation, we have used the fact that at $T=0$, $f_{n\bm{k}}^{\prime}=\delta(\varepsilon_{n\bm{k}}-\mu)$. The total contribution from the $y$-component of the Zeeman field thus is

[TABLE]

We can also calculate the contribution from the $x$-component of the Zeeman field. The final result reads

[TABLE]

where $v_{nn^{\prime}}=\langle n\bm{k}|\hat{\bm{v}}|n^{\prime}\bm{k}\rangle$ is the interband element of the velocity operator. This is the multi-band formula for the spin toroidization. It reduces to Eq. (12) in the main text for a single band.

Appendix B Wannier representation

In this section we express the spin toroidization in terms of the Wannier functions. Denote by $|w_{0}(\bm{R},\bm{B})\rangle$ the Wannier function located at the lattice site $\bm{R}$ from band [math], derived from the local Hamiltonian $\hat{H}_{c}$. The periodic part of the Bloch function $|\tilde{u}_{0}\rangle$ is given by

[TABLE]

where $N$ is the number of unit cells.

We begin with the spin toroidization formula in Eq. (35), which can be recast as

[TABLE]

The first term in Eq. (37) can be expressed in terms of $|\tilde{u}_{0}\rangle$,

[TABLE]

where we have used the identity: $|\tilde{u}_{0}\rangle\langle\tilde{u}_{0}|+\sum_{n\neq 0}|\tilde{u}_{n}\rangle\langle\tilde{u}_{n}|=I$. Inserting Eq. (36) into the above expression yields

[TABLE]

Here $|w_{0}(\bm{B})\rangle=|w_{0}(\bm{R},\bm{B})\rangle$ with $\bm{R}=0$.

Next we turn to the second term of Eq. (37). Using $-g\mu_{B}\hat{\bm{s}}/\hbar=\bm{\partial}_{\bm{B}}\hat{H}_{c}$, we have,

[TABLE]

and

[TABLE]

The second term in Eq. (37) then becomes

[TABLE]

Combining Eq. (39) and Eq. (B), we obtain the final expression in Eq. (14) fo the spin toroidization in the Wannier representation.

Appendix C Molecular insulator limit

Under this limit the first two terms in Eq. (14) becomes:

[TABLE]

where $\bar{\bm{r}}=\langle w_{0}(\bm{B})|\bm{r}|w_{0}(\bm{B})\rangle$ is the expectation value of the position, and $\bar{\bm{s}}=\langle w_{0}(\bm{B})|\hat{\bm{s}}|w_{0}(\bm{B})\rangle$ is the expectation value of the spin. Since the combined time reversal and space inversion symmetry is respected, we must have $\bar{\bm{s}}=0$. Therefore,

[TABLE]

Now we consider the remaining term in Eq. (14). In the molecular insulator limit, its form does not change. Note that $|w_{0}(\bm{B})\rangle$ and $|w_{n}(\bm{B})\rangle$ becomes the molecular eigenfunctions and $\varepsilon_{0}$ and $\varepsilon_{n}$ become the molecular eigenenergy. We further manipulate this term as follows:

[TABLE]

where $\bar{\bm{r}}$ has been defined before, and stands for the position of electron under external magnetic field. Here $\varepsilon_{0}-\mu$ is the free energy for state [math].

Therefore, the total toroidization in the molecular insulator limit reads:

[TABLE]