Frustration and chiral orderings in correlated electron systems

As we have seen in the previous section, certain spin orderings lead to non-trivial charge effects in Mott Insulators. In particular, orbital currents arise from nonzero scalar spin chirality: $\langle {{\chi}_{jkl}}\rangle =\langle {{\mathbf{S}}_{j}}\cdot {{\mathbf{S}}_{k}}\times {{\mathbf{S}}_{l}}\rangle \ne 0$. While chiral ordering can exist even in absence of magnetic ordering ($\langle {{\mathbf{S}}_{j}}\rangle =0$), it is clear that non-coplanar magnetic ($\langle {{\mathbf{S}}_{j}}\rangle \ne 0$) orderings have a non-zero local chirality $\langle {{\chi}_{jkl}}\rangle \ne 0$. The purpose of this section is to show how non-coplanar magnetic orderings can emerge from frustration in local moment systems. We will then consider the simple example of a triangular spin S Heisenberg model with a FM nearest-neighbor exchange, J₁  <  0, and an AFM third-nearest-neighbor exchange J₃:

$$\begin{eqnarray}&&{{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}={{J}_{1}}\underset{\langle \,j,l\rangle}{\sum}{{\mathbf{S}}_{j}}\cdot {{\mathbf{S}}_{l}}+{{J}_{3}}\underset{\langle \langle \,j,l\rangle \rangle}{\sum}{{\mathbf{S}}_{j}}\cdot {{\mathbf{S}}_{l}}-H\underset{j}{\sum}S_{j}^{z}.\end{eqnarray} \tag{ 10 }$$

The first sum runs over nearest-neighbor bonds $\langle \,j,l\rangle$, while the second sum runs over third-nearest-neighbor bonds $\langle \langle \,j,l\rangle \rangle$. From now on, we will use the lattice parameter, a, as our length unit.

The relative position vectors of the nearest-neighbors of a given atom in the triangular lattice are $\pm {{\mathbf{e}}_{\nu}}$ ($\nu =1,2,3$), with ${{\mathbf{e}}_{1}}=\hat{\mathbf{x}}$, ${{\mathbf{e}}_{2}}=-\hat{\mathbf{x}}/2+\sqrt{3}\hat{\mathbf{y}}/2$ and ${{\mathbf{e}}_{2}}=-\hat{\mathbf{x}}/2-\sqrt{3}\hat{\mathbf{y}}/2$. For what follows, the main role of the competing real space interactions, J₁ and J₃, is to induce an exchange interaction in momentum space,

$$\begin{eqnarray}&&J\left(\mathbf{q}\right)=\underset{j=1,3}{\sum}\left({{J}_{1}}\cos \mathbf{q}\cdot {{\mathbf{e}}_{j}}+{{J}_{3}}\cos 2\mathbf{q}\cdot {{\mathbf{e}}_{j}}\right),\end{eqnarray} \tag{ 11 }$$

with a local maximum at $\mathbf{q}=\mathbf{0}$ and global minima at six wave-vectors $\pm {{\mathbf{Q}}_{\nu}}$ ($\nu =1,2,3$) connected by the C₆ symmetry group of the triangular lattice. The magnitude of these vectors is:

$$\begin{eqnarray}&&Q=\,|{{\mathbf{Q}}_{\nu}}|\,=2\arccos \left[\frac{1}{4}\left(1+\sqrt{1-\frac{2{{J}_{1}}}{{{J}_{3}}}}\right)\right].\end{eqnarray} \tag{ 12 }$$

According to equation (12), the system is a ferromagnet (Q  =  0), for ${{J}_{3}}<-{{J}_{1}}/4$. The ordering wave-vector becomes nonzero for ${{J}_{3}}>-{{J}_{1}}/4$. Therefore, ${{J}_{3}}=-{{J}_{1}}/4$ is the Lifshitz transition point from FM (Q  =  0) to incommensurate ($Q\ne 0$) AFM ordering. Correspondingly, the saturation field (magnetic field required to fully saturate the magnetic moments) also becomes finite for ${{J}_{3}}>-{{J}_{1}}/4$:

$$\begin{eqnarray}&&{{H}_{\text{sat}}}=-SJ\left({{\mathbf{Q}}_{\mu}}\right)+3S\left({{J}_{3}}+{{J}_{1}}\right).\end{eqnarray} \tag{ 13 }$$

Clearly, ${{H}_{\text{sat}}}=0$ for ${{J}_{3}}<-{{J}_{1}}/4$ because $J\left(\mathbf{0}\right)=3\left({{J}_{1}}+{{J}_{3}}\right)$. Right above ${{J}_{3}}=-{{J}_{1}}/4$, we can expand in the small parameter $\delta =3\left({{J}_{1}}+4{{J}_{3}}\right)/\left(2{{J}_{3}}\right)$:

$$\begin{eqnarray}&&{{H}_{\text{sat}}}\simeq 6S{{\delta}^{2}}J_{3}^{2}/\left({{J}_{1}}+32{{J}_{3}}\right).\end{eqnarray} \tag{ 14 }$$

In addition, we can expand equation (12) to obtain:

$$\begin{eqnarray}&&Q\simeq \frac{2}{3}\sqrt{\delta},\end{eqnarray} \tag{ 15 }$$

which implies that H_sat is proportional to Q⁴ near the Lifshitz point ${{J}_{3}}=-{{J}_{1}}/4$.

We will assume now that ${{J}_{3}}>-{{J}_{1}}/4$ and that the magnetic field H is larger than H_sat. We want to know the effect of replacing one spin by a non-magnetic impurity. It has been shown recently that non-magnetic impurities lead to non-coplanar spin structures in the triangular Heisenberg model with a nearest-neighbor AFM interaction [41–43]. In contrast, here we will consider the effect of non-magnetic impurities on frustrated magnets with competing FM and AFM interactions. Given the FM nature of J₁, the saturation field of the spins near the impurity must be higher than the saturation field H_sat far away from the impurity. This simple observation implies that the spins should cant away from the z-axis around the non-magnetic impurity.

The spin texture far away from the non-magnetic impurity can be obtained by taking the continuum limit $Q=\,|{{\mathbf{Q}}_{\nu}}|\ll 1$. In this limit, ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$ becomes

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}} & = & \mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{iso}}+\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{ani}}, \end{array}\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{iso}} & = & {\int}^{}\text{d}{{r}^{2}}\left[-\frac{\delta}{2}{{J}_{3}}{{\left(\nabla \mathbf{S}\right)}^{2}}+\frac{3}{64}\left({{J}_{1}}+16{{J}_{3}}\right){{\left({{\nabla}^{2}}\mathbf{S}\right)}^{2}}-\mathbf{H}\cdot \mathbf{S}\right], \end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{ani}}= & \frac{1}{7680}{\int}^{}\text{d}{{r}^{2}}({{J}_{1}}+64{{J}_{3}}) \\ {} & \times [11{{(\partial _{x}^{3}\mathbf{S})}^{2}}+15(\partial _{x}^{2}{{\partial}_{y}}\mathbf{S}){{}^{2}}+45{{({{\partial}_{x}}\partial _{y}^{2}\mathbf{S})}^{2}}+9{{(\partial _{y}^{3}\mathbf{S})}^{2}}], \end{array}\end{eqnarray} \tag{ 18 }$$

up to sixth order in a gradient expansion and $A=({{J}_{1}}+64{{J}_{3}})/7680$. Given that the leading order amplitude of the C₆ anisotropy is Q⁶, we can neglect ${{\mathcal{H}}^{\text{ani}}}$ for $Q\ll 1$. In addition, we will rescale our energy, length and magnetic field units in order to normalize the coefficients [44]

$$\begin{eqnarray}&&\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{iso}}={\int}^{}\left[-\frac{\delta}{2}{{\left(\nabla \mathbf{S}\right)}^{2}}+\frac{1}{2}{{\left({{\nabla}^{2}}\mathbf{S}\right)}^{2}}-\mathbf{H}\cdot \mathbf{S}\right]\text{d}{{r}^{2}}.\end{eqnarray} \tag{ 19 }$$

In the new units the saturation field is ${{H}_{\text{sat}}}={{\delta}^{2}}/4$ and $Q=\sqrt{\delta}/\sqrt{2}$. Equation (19) is the universal theory for centrosymmetric frustrated magnets near a Lifshitz point.

To model the effect of the non-magnetic impurity located at the origin, we simply need to modify the stiffness $-\delta$ near the origin. We note that $\delta =3\left({{J}_{1}}+4{{J}_{3}}\right)/\left(2{{J}_{3}}\right)$ is enhanced near the origin because of the suppression of J₁ on the bonds connecting the non-magnetic impurity and its nearest neighbors. Therefore, we will replace the spatially uniform δ with

$$\begin{eqnarray}&&\tilde{\delta}(r)=\delta \left[1+{{\Delta}_{0}}\Theta\left({{r}_{0}}-r\right)\right],\end{eqnarray} \tag{ 20 }$$

where $\Theta(x)$ is the Heaviside step function, r₀ is the range over which the local stiffness is modified by the presence of the impurity, and ${{\Delta}_{0}}$ is a positive dimensionless parameter. We note that a similar continuum model is valid for other types of impurities. The choice of the step function is arbitrary because details of this non-universal function do not affect the asymptotic behavior of the solution for $r\gg {{r}_{0}}$.

As we explained above, the saturation field near the impurity, $H_{\text{sat}}^{I}$, is bigger than the bulk saturation field ${{H}_{\text{sat}}}={{\delta}^{2}}/4$. The nature of the magnetic distortion around the impurity is obtained from the structure of the lowest energy internal mode (localized mode around the impurity), which becomes gapless at $H=H_{\text{sat}}^{I}$. We will then compute the normal modes of the fully polarized phase in order to determine this instability. To describe the small oscillations around the ground state configuration, we introduce the bosonic field $\boldsymbol{\varphi}\left(\mathbf{r}\right)$ via the Holstein-Primakoff transformation

$$\begin{eqnarray}\begin{array}{*{35}{l}} S_{\mathbf{r}}^{\dagger} & = & S_{\mathbf{r}}^{x}+\text{i}S_{\mathbf{r}}^{y}=\sqrt{2S-\boldsymbol{\varphi}^{{\dagger}}\left(\mathbf{r}\right)\boldsymbol{\varphi}\left(\mathbf{r}\right)}\boldsymbol{\varphi}\left(\mathbf{r}\right), \\ S_{\mathbf{r}}^{z} & = & S-\boldsymbol{\varphi}^{{\dagger}}\left(\mathbf{r}\right)\boldsymbol{\varphi}\left(\mathbf{r}\right). \end{array}\end{eqnarray} \tag{ 21 }$$

By expanding $\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{iso}}$ up to quadratic order in $\boldsymbol{\varphi}\left(\mathbf{r}\right)$, we obtain the spin-wave Hamiltonian density:

$$\begin{eqnarray}&&\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{sw}}=H\boldsymbol{\varphi}^{{\dagger}}\boldsymbol{\varphi}-\tilde{\delta}(r)S\nabla \boldsymbol{\varphi}^{{\dagger}}\nabla \boldsymbol{\varphi}+S{{\nabla}^{2}}\boldsymbol{\varphi}^{{\dagger}}{{\nabla}^{2}}\boldsymbol{\varphi}.\end{eqnarray} \tag{ 22 }$$

The resulting equation of motion for $\boldsymbol{\varphi}\left(\mathbf{r}\right)$ is:

$$\begin{eqnarray}&&-\text{i}{{\partial}_{t}}\boldsymbol{\varphi}=H\boldsymbol{\varphi}-\tilde{\delta}(r)S{{\nabla}^{2}}\boldsymbol{\varphi}+S{{\nabla}^{4}}\boldsymbol{\varphi}.\end{eqnarray} \tag{ 23 }$$

In absence of the impurity (translationally invariant case), the eigenmodes are magnons with a dispersion relation,

$$\begin{eqnarray}&&{{\omega}_{\mathbf{k}}}=-\delta S{{k}^{2}}+S{{k}^{4}}+H=S{{\left({{k}^{2}}-\frac{\delta}{2}\right)}^{2}}+H-\frac{S{{\delta}^{2}}}{4},\end{eqnarray} \tag{ 24 }$$

obtained by Fourier transforming the equation of motion equation (23). This expression shows explicitly that ${{H}_{\text{sat}}}=S{{\delta}^{2}}/4$, i.e. the magnon gap is ${{\Delta}_{s}}=H-{{\delta}^{2}}/4$. We note that the magnetic field enters as an additive constant that shifts the whole spectrum. This is so because it couples to the conserved quantity $S_{T}^{z}={\int}^{}S_{\mathbf{r}}^{z}\text{d}{{r}^{2}}$ (total magnetization along the z-axis).

In the presence of the impurity, we need to replace the spatially uniform stiffness $-\delta$ by the non-uniform stiffness, $-\tilde{\delta}$, given by equation (20). The eigenmodes, $\boldsymbol{\psi}^{{\dagger}}|0\rangle$, are now created by operators of the form $\boldsymbol{\psi}^{{\dagger}}={\int}^{}{{\psi}^{\ast}}\left(\mathbf{r}\right)\boldsymbol{\varphi}^{{\dagger}}\left(\mathbf{r}\right)\text{d}{{r}^{2}}$, where $\psi \left(\mathbf{r}\right)$ is an eigenfunction of the operator $H-S\tilde{\delta}(r){{\nabla}^{2}}+S{{\nabla}^{4}}$. While translational symmetry is broken by the presence of the impurity, rotational symmetry is still preserved. Consequently, it is useful to introduce polar coordinates $\mathbf{r}=r\left(\cos \phi,\sin \phi \right)$ and work in the basis of eigenstates of $\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{sw}}$ with well defined angular momentum l. These circular waves are described by the Bessel functions of the first kind:

$$\begin{eqnarray}&&{{\nabla}^{2}}\left[{{J}_{l}}(kr){{\text{e}}^{\text{i}l\phi}}\right]=-{{k}^{2}}{{J}_{l}}(kr){{\text{e}}^{\text{i}l\phi}}.\end{eqnarray} \tag{ 25 }$$

For $r\leqslant {{r}_{0}}$, the stiffness is constant and equal to $-\delta \left(1+{{\Delta}_{0}}\right)$. Therefore, the eigenmodes in this region are linear combinations of two circular waves

$$\begin{eqnarray}&&\psi \left(r\leqslant {{r}_{0}}\right)={{A}_{1}}{{J}_{l}}\left({{q}_{+}}r\right){{\text{e}}^{\text{i}l\phi}}+{{B}_{1}}{{J}_{l}}\left({{q}_{-}}r\right){{\text{e}}^{\text{i}l\phi}},\end{eqnarray} \tag{ 26 }$$

with the two momenta, ${{q}_{\pm}}$, that produce the same energy eigenvalue ω:

$$\begin{eqnarray}&&{{q}_{\pm}}=\frac{\delta \left(1+{{\Delta}_{0}}\right)\pm \sqrt{{{\delta}^{2}}-4H/S-4\omega /S+2{{\delta}^{2}}{{\Delta}_{0}}+{{\delta}^{2}}\Delta _{0}^{2}}}{\sqrt{2}}.\end{eqnarray} \tag{ 27 }$$

Because we are looking for bound states, the corresponding eigenmodes must decay exponentially for $r\to \infty$. Therefore, we need to use the modified Bessel functions of the second kind for r  >  r₀, which are also eigenvectors of the Laplacian operator:

$$\begin{eqnarray}&&{{\nabla}^{2}}\left[{{K}_{l}}(kr){{\text{e}}^{\text{i}l\phi}}\right]={{k}^{2}}{{K}_{l}}(kr){{\text{e}}^{\text{i}l\phi}}.\end{eqnarray} \tag{ 28 }$$

Once again, the eigenvector for r  >  r₀ is a linear combination of two of these functions

$$\begin{eqnarray}&&\psi \left(r>{{r}_{0}}\right)={{A}_{2}}{{K}_{l}}\left({{k}_{+}}r\right){{\text{e}}^{\text{i}l\phi}}+{{B}_{2}}{{K}_{l}}\left({{k}_{-}}r\right){{\text{e}}^{\text{i}l\phi}},\end{eqnarray} \tag{ 29 }$$

with momenta ${{k}_{\pm}}$ that produce the energy eigenvalue ω and lead to an exponential decay for $r\to \infty$:

$$\begin{eqnarray}&&{{k}_{\pm}}=\frac{-\delta \pm \sqrt{{{\delta}^{2}}-4H/S-4\omega /S}}{\sqrt{2}}.\end{eqnarray} \tag{ 30 }$$

We note that ${{\delta}^{2}}-4H/S-4\omega /S>0$ because we are looking for bound states, which must lie below the magnon gap $H-S{{\delta}^{2}}/4$. The last part of the calculation is to ask that the eigenmodes and their three first derivatives must be continuous at r  =  r₀. This condition arises from the fourth order nature of the differential equation (23).

At this point, it is clear that the problem under consideration is analogous to the quantum mechanical problem of a single particle moving in 2D with an effective mass that depends on its distance to the origin (it takes one value for r  >  r₀ and a different value for $r\leqslant {{r}_{0}}$). However, an important difference relative to the usual single-particle problem is that the energy minimum occurs at a finite wave-vector k  =  Q (see equation (24)): ${{\omega}_{\mathbf{k}}}$ is flat on the circle for radius k  =  Q and it has a finite dispersion along the radial direction. This implies that the density of single magnon states, $\rho \left(\omega \right)$, has the same singularity as a one-dimensional (1D) system when ω approaches the bottom of the single magnon dispersion $\omega ={{\Delta}_{s}}$: $\rho \left(\omega \right)\propto 1/\sqrt{\omega -{{\Delta}_{s}}}$ for $\omega \to {{\Delta}_{s}}$. As it is well known, this singular behavior of $\rho \left(\omega \right)$ leads to the formation of a bound state for an infinitesimal attractive potential well. In particular, the 1D-like divergence in $\rho \left(\omega \right)$ produces a binding energy, ${{\Delta}_{b}}={{\omega}_{b}}-\Delta$, which is proportional to the square of the amplitude of the attractive interaction, ${{\Delta}_{b}}\propto -\Delta _{0}^{2}$, in contrast to the essential singularity that appears for the usual 2D problem with a single minimum at k  =  0.

In the case under consideration, the effective attraction is produced by an increase in the absolute value of the stiffness for r  <  r₀, which is felt by all the circular waves regardless of their angular momenta l. Given that the density of states exhibits the same singularity for any value of l, we expect the formation of a bound state for each value of l. The l-value of the lowest energy bound state determines the nature of the magnetic distortion around the impurity right below $H=H_{\text{sat}}^{I}$.

For weak attraction, ${{\Delta}_{0}}\ll 1$, the amplitude of the attractive interaction in the l channel is

$$\begin{eqnarray}&&{{g}_{l}}\left({{r}_{0}}\right)={{Q}^{3}}S{{\Delta}_{0}}{\int}_{0}^{{{r}_{0}}}J_{l}^{2}(Qr)\text{d}r.\end{eqnarray} \tag{ 31 }$$

Given the asymptotic form of the Bessel functions for small argument, ${{J}_{l}}(z)\simeq {{(z/2)}^{l}}/\Gamma(l+1)$ for $0\ll z\ll \sqrt{l+1}$, g_l is maximized for l  =  0 if ${{r}_{0}}\ll 1$. However, it is easy to verify that this is not necessarily true when r₀ Q becomes bigger than one (see inset of figure 3). In particular, g_l acquires its maximum value for $l=\pm 1$ when r₀Q becomes of the order or bigger than the first root of J₁(z), as it is shown in figure 3. Moreover, a generalized impurity that increases $|\delta (r)|$ in the ring $R-{{r}_{0}}/2<r<R+{{r}_{0}}/2$ will maximize g_l for $|l|$-values, which increase monotonically with R. This type of impurity corresponds to removing spins from a ring of sites, instead of the single-site (R  =  r₀/2) that we are currently considering [45].

Zoom In Zoom Out Reset image size

Figure 4 shows the bound state energies of the l  =  0, 1, 2 channels as a function of ${{\Delta}_{0}}$ for r₀  =  2 and for r₀  =  4. Interestingly enough, the l  =  1 bound state becomes the ground state above a critical value of the impurity potential ${{\Delta}_{0}}$ for r₀  =  2, while it is already the ground state for arbitrarily small values of ${{\Delta}_{0}}$ if r₀  =  4. This result implies that the leading instability around an impurity in a frustrated magnet can be a magnetic vortex with winding number $l=\pm 1$. Moreover, strong impurity potentials (larger values of ${{\Delta}_{0}}$) can produce lowest energy bound states with even higher values of l, i.e. vortices with winding numbers larger than on will emerge below $H=H_{\text{sat}}^{I}$.

Zoom In Zoom Out Reset image size

The main conclusion of our semi-classical stability analysis is that impurities can nucleate magnetic vortices in frustrated magnets well above the bulk saturation field. This conclusion can be explicitly verified in the classical limit. To find the actual distortion induced by the non-magnetic impurity below $H=H_{\text{sat}}^{I}$, we solve numerically the Landau–Lifshitz–Gilbert equation of motion for the classical limit of ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$:

$$\begin{eqnarray}&&{{\partial}_{t}}\mathbf{S}=-\mathbf{S}\times {{\mathbf{H}}_{\text{eff}}}+\alpha \mathbf{S}\times {{\partial}_{t}}\mathbf{S}.\end{eqnarray} \tag{ 32 }$$

Here α is the Gilbert damping and ${{\mathbf{H}}_{\text{eff}}}\equiv -\delta \mathcal{H}/\delta \mathbf{S}$ is the effective magnetic field acting on each spin. The non-magnetic impurity is introduced by setting the spin at the origin to 0: $\mathbf{S}\left(\mathbf{r}=0\right)=0$. Consistently with our previous analysis, a vortex is nucleated around the impurity once the system is allowed to relax from an initial fully polarized state.

As it is shown in figure 5, the linear vortex size increases upon approaching the bulk saturation field $H={{H}_{\text{sat}}}$. Indeed, by solving the Euler-Lagrange equations of the continuum model for $H_{\text{sat}}^{I}>H>{{H}_{\text{sat}}}$, one can verify that the vortex amplitude (tilting of the spins away from the z-axis) decays exponentially over the magnetic correlation length ξ, which diverges as $\xi \propto 1/\sqrt{H-{{H}_{\text{sat}}}}$ upon approaching the bulk saturation field H_sat. In other words, the vortex radius diverges at the critical point $H={{H}_{\text{sat}}}$, where the exponential decay is replaced by an algebraic decay $1/\sqrt{r}$ [45]. This algebraic decay signals a second order transition into a conical single-$\mathbf{Q}$ magnetic ordering that will be considered in the next sections.

Zoom In Zoom Out Reset image size

Finally, we note that a similar calculation can be easily extended to quantum spin systems of arbitrary spin S. In this case, the modes for $H>H_{\text{stat}}^{I}$ are obtained by diagonalizing the quantum version of ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$ in the $S_{T}^{z}=NS-1$ sector, i.e. in the subspace of states with a single-spin flip ($S\to S-1$) relative to the fully polarized ground state. As it is illustrated in figure 6, the flipped spin can be regarded as a single particle moving in the central potential generated by the impurity at the origin. If we assume that the impurity consists of a smaller magnetic moment ${{S}^{\prime}}$ (the non-magnetic impurity corresponds to ${{S}^{\prime}}=0$), the flipped spin has a lower energy, ${{J}_{1}}\left(S-{{S}^{\prime}}\right)$ (J₁  <  0) when sitting on the first hexagon of nearest-neighbor sites around the impurity (potential well) and a higher energy ${{J}_{3}}\left(S-{{S}^{\prime}}\right)$ (potential barrier) when sitting on the six third-nearest neighbors (see figure 6). The hopping amplitude between a pair of sites j and l away from the origin is J_jl S. The hopping amplitude between the impurity site and a different site l is ${{J}_{0l}}\sqrt{S{{S}^{\prime}}}$. It is then clear that the total spin S is an overall scaling factor for the effective single-particle Hamiltonian. In other words, the normal modes depend only on the ratio $\alpha ={{S}^{\prime}}/S$, implying that for the case of a non-magnetic impurity (${{S}^{\prime}}=0$) the normal modes are exactly the same all the way from the S  =  1/2 quantum limit to the $S\to \infty$ classical limit. This fact remains true for any finite ${{S}^{\prime}}$ as long as we keep the ${{S}^{\prime}}/S$ ratio fixed while varying S, and it explains why the equation of the motion for the eigenmodes in the classical spin limit (23) is formally the same as the single particle Schrödinger equation.

Zoom In Zoom Out Reset image size

As expected from our previous analysis, the ground state for the single spin-flip subspace of the quantum ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$ Hamiltonian with a non-magnetic impurity (${{S}^{\prime}}=0$) corresponds to two degenerate bound states with quasi-angular momentum $l=\pm 1$ (the phase of the bound state wave function has a winding number $l=\pm 1$) (see figure 7). As we explained before, this bound is the precursor of the vortex state that appears around the impurity below the magnetic field, $H=H_{\text{sat}}^{I}$, at which the energy of the bound state becomes equal to the energy of the fully polarized state: $H_{\text{sat}}^{I}={{H}_{\text{sat}}}-{{\Delta}_{b}}$. The absolute value of the binding energy, $|{{\Delta}_{b}}|$, reaches its maximum value near ${{J}_{3}}/{{J}_{1}}\simeq -0.6$ and it decreases for larger values of ${{J}_{3}}/{{J}_{1}}$ when the system approaches the other incommensurate to commensurate transition for ${{J}_{3}}/{{J}_{1}}\to \infty$. This is indeed the expected behavior given that the potential well disappears for ${{J}_{1}}\to 0$ (${{\epsilon}_{1}}\to 0$) and only the potential barrier (${{\epsilon}_{3}}={{J}_{3}}S$) remains.

Zoom In Zoom Out Reset image size

The inset of figure 7 shows the ratio $\left(H_{\text{sat}}^{I}-{{H}_{\text{sat}}}\right)/{{H}_{\text{sat}}}$ as a function of ${{J}_{3}}/{{J}_{1}}$. The window of magnetic field values where the non-magnetic impurity is expected to bind a vortex is a sizable fraction of the bulk saturation field, H_sat, for ${{J}_{3}}\lesssim {{J}_{1}}$. While we have used a particular Hamiltonian, ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$, for describing the creation of magnetic vortices by non-magnetic impurities above the bulk saturation field, our conclusion is valid for a larger set of frustrated Hamiltonians, which exhibit a Lifshitz transition from a commensurate FM state to incommensurate helical ordering. The FM interaction, represented by J₁ in our model, is necessary to have $H_{\text{sat}}^{I}>{{H}_{\text{sat}}}$, i.e. to have a bound state around the impurity. The underlying lattice does not need to have C₆ symmetry. Vortex sates around non-magnetic impurities should also exist on tetragonal systems (C₄ symmetry) [45].

NiGa₂S₄ provides an experimental realization of $\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{sw}}$ with a ratio $-{{J}_{3}}/{{J}_{1}}=0.2(1)$ [46]. Unfortunately this small ratio leads to an extremely narrow magnetic field window above H_sat for observing the vortex-impurity bound state. NiBr₂ is an alternative realization of $\mathcal{H}_{{{J}_{1}}-{{J}_{3}}}^{\text{sw}}$ with $-{{J}_{3}}/{{J}_{1}}=0.26$ [47]. The saturation field is very high (${{H}_{\text{sat}}}=380$ KOe) in this compound because of the presence of a rather strong inter-layer AFM exchange. Non-magnetic impurities can by introduced in this compound by replacing Ni with Zn [48]. Based on the results presented in this section, we predict that Zn_xNi_1−xBr₂ should exhibit magnetic vortices bounded to the Zn-impurities right above the saturation field. Moreover, these vortices should persist below the saturation field [45]. This observation could explain the results of neutron diffraction experiments on a single-crystal of Zn_xNi_1−xBr₂ containing 8-mole% of Zn²⁺ , which show that the helical propagation wave-vector direction becomes 'disordered', i.e. the incommensurate peaks of the magnetic structure factor form a regular hexagonal structure around the commensurate wave vector [48].

CeRhAl₄Si₂ and CeIrAl₄Si are examples of tetragonal frustrated magnets with competing FM and AFM interactions, which exhibit incommensurate intra-layer magnetic ordering at zero field and low enough temperature [49, 50]. Given that these compounds are easy-axis antiferromagnets, magnetic vortices around non-magnetic impurities should be observed above the meta-magnetic transition induced by a magnetic field parallel to the c-axis. We also note that the metallic character of these materials, whose local moments are provided by localized f-electrons that interact with each other via the RKKY interaction, does not modify our previous analysis. The theory presented in this section applies both to MI's and to itinerant magnets.

It is important to mention that the nearest-neighbor FM exchange also leads to the formation of two magnon bound states above H_sat [51–53]. This implies that for finite spin S, the bulk saturation field is in general higher than the value associated with a single-magnon condensation (see equation (13)). The attractive magnon–magnon interaction can lead to a continuous transition into some multipolar ordering (e.g. nematic ordering for the condensation of magnon pairs) [54–59] or to a discontinuous transition. In any case, the important observation is that the magnitude of the magnon–magnon interaction is of order one, while the attractive interaction, ${{\epsilon}_{1}}$, between the magnon and a non-magnetic impurity is proportional to S, implying that $H_{\text{sat}}^{I}$ remains higher than H_sat for large enough S. Another factor that makes the magnon-impurity binding energy larger than the magnon–magnon binding energy is the static nature of the impurity: the reduced mass for the two-magnon problem is half of the single magnon mass relevant for the magnon-impurity bound state formation.

Finally, it is worth mentioning that impurities can also nucleate vortices above the saturation field of frustrated magnets with two or more competing AFM interactions, as long as the saturation field near the impurity remains higher than the bulk saturation field. This condition requires that the impurity must distort the lattice locally in order to enhance the the AFM exchange interactions on the neighboring bonds.

3.2.1. Periodic array of impurities.

In the previous section we found that a magnetic vortex forms around a single impurity in a finite interval of magnetic field above the saturation field. We will consider now the effect of a periodic array of impurities. This situation can be realized in a doped system with a finite concentration of holes if the inter-hole Coulomb interaction induces strong charge ordering, i.e. the holes become strongly localized in a periodic structure. This situation has been observed in some high-T_c and related materials for a hole concentration x  =  1/8, [60, 61] although the Heisenberg term of the t  −  J model that describes these materials is not frustrated [62]. The question that we will address now is: what is the effect of a charge-density-wave (CDW) ordering on the magnetic structure in frustrated systems? While the answer to this question will of course depend on the nature and period of the CDW structure, our only purpose here is to show that spontaneous chiral magnetic orderings emerge quite naturally out of the interplay with the underlying charge ordering.

To simplify the discussion we will consider the simple case of a periodic array of non-magnetic sites (holes) forming a triangular superlattice structure with primitive reciprocal lattice vectors ${{\mathbf{K}}_{+}}$ and ${{\mathbf{K}}_{-}}$. To avoid frustration between the magnetic and the charge ordering, we will choose ${{\mathbf{K}}_{+}}$ and ${{\mathbf{K}}_{-}}$ to be commensurate with the magnetic ordering wave vectors ${{\mathbf{Q}}_{\nu}}$. In particular, we will set ${{J}_{1}}/{{J}_{3}}\sim -1.171$ (for which ${{\mathbf{Q}}_{1}}=\left(\pi /2,0\right)$) and we assume a CDW with ${{\mathbf{K}}_{\pm}}=\pi /4\left(1,\pm 1/\sqrt{3}\right)$, which, as depicted in the inset of figure 8(b), corresponds to a triangular superlattice with the lattice spacing eight times larger than the original lattice parameter. We note also that ${{\mathbf{Q}}_{1}}={{\mathbf{K}}_{+}}+{{\mathbf{K}}_{-}}$ holds in this case.

Zoom In Zoom Out Reset image size

We extend the ${{J}_{1}}-{{J}_{3}}$ Hamiltonian (equation (10)) to include the nonmagnetic impurities with the prescription,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{S}}_{j}} & \to {{p}_{j}}{{\mathbf{S}}_{j}}, \end{array}\end{eqnarray} \tag{ 33 }$$

where p_j  =  0 (1) corresponds to the presence (absence) of a quenched hole or non-magnetic impurity at the site i. Figure 8 shows the thermodynamic phase diagrams obtained from Monte Carlo (MC) simulations in the absence (figure 8(a)) and the presence of the CDW (figure 8(b)). The different phases of these phase diagrams are characterized by the spin structure factor,

$$\begin{eqnarray}&&S_{s}^{\nu \nu}\left(\mathbf{q}\right)=\frac{1}{N}\underset{j,l}{\sum}\langle S_{j}^{\nu}S_{l}^{\nu}\rangle {{\text{e}}^{\text{i}\mathbf{q}\cdot \left({{\mathbf{r}}_{j}}-{{\mathbf{r}}_{l}}\right)}},\end{eqnarray} \tag{ 34 }$$

and the chirality structure factor, which we define for both the up and the down triangles ($\mu =u,d$, respectively) as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} S_{\chi}^{\mu}\left(\mathbf{q}\right)=\frac{1}{N}\underset{\mathbf{R},{{\mathbf{R}}^{\prime}}\in \mu}{\sum}\langle {\chi}_{\mathbf{R}}^{}{{\chi}_{{{\mathbf{R}}^{\prime}}}}\rangle {{\text{e}}^{\text{i}\mathbf{q}\cdot \left(\mathbf{R}-{{\mathbf{R}}^{\prime}}\right)}}, \end{array}\end{eqnarray} \tag{ 35 }$$

where $\mathbf{R},{{\mathbf{R}}^{\prime}}$ run over sites of the specified ($\mu =u,d$) sublattice of the dual honeycomb lattice. ${\chi}_{\mathbf{R}}^{}$ is the spin scalar chirality associated with a triangle centered at $\mathbf{R}$ given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\chi}_{\mathbf{R}}}={{\mathbf{S}}_{j}}\cdot \left({{\mathbf{S}}_{k}}\times {{\mathbf{S}}_{l}}\right), \end{array}\end{eqnarray} \tag{ 36 }$$

where j, k, and l are the sites aligned counterclockwise on this triangle. We also introduce the following notations for the total scalar chirality from the upward (${{\chi}^{u}}$) and the downward (${{\chi}^{d}}$) triangles and their sum:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\chi}^{\mu}} & =\frac{1}{N}\underset{\mathbf{R}\in \mu}{\sum}{\chi}_{\mathbf{R}}^{}~\text{for}~\mu =u,d, \\ {{\chi}^{\text{tot}}} & ={{\chi}^{u}}+{{\chi}^{d}}. \end{array}\end{eqnarray} \tag{ 37 }$$

The phase diagram in absence of the CDW is rather simple (see figure 8(a)), including only the single-$\mathbf{Q}$ conical spiral phase in addition to the paramagnetic state. In this phase, the xy spin components are quasi-long-range ordered in a magnetic field and they develop long range order (LRO) only at T  =  0. However, the C₆ symmetry is broken at finite T and H because of the antisymmetric bond ordering characterized by the bond parameter $\langle {{\mathbf{S}}_{i}}\times {{\mathbf{S}}_{j}}\cdot \hat{\mathbf{z}}\rangle$ (note that this order parameter breaks only discrete symmetries).

The phase diagram becomes much more complex in the presence of the CDW (figure 8(b)). The main qualitative difference relative to the single-$\mathbf{Q}$ conical spiral phase obtained in absence of the CDW is the emergence of a net chiral component, ${{\chi}^{\text{tot}}}\ne 0$, in most of the new phases. Exceptions to this rule are the vertical 1Q^M spiral ordering and the antiferrochiral 1Q^M spiral phases that appear next to the paramagnetic phase in the low field region. These are also the only two phases that exhibit single-$\mathbf{Q}$ magnetic ordering. Notably, many phases support long-wavelength modulation of local scalar chirality (chirality wave), which is not necessarily single-${{Q}^{\chi}}$ (hereafter, our convention is to use Q^M and ${{Q}^{\chi}}$, respectively, when it is necessary to make an unambiguous distinction between spin and chirality textures). The stability of each phase in thermodynamic limit is confirmed by performing a finite-size scaling analysis presented in [63].

Figure 9 shows a sequence of snapshots of the spin configurations of the low-temperature phases from zero (figure 9(a)) to the saturation field (figure 9(j)). We also show the square root of the corresponding xy ($\sqrt{S_{s}^{\bot}\left(\mathbf{q}\right)}$) and z ($\sqrt{S_{s}^{zz}\left(\mathbf{q}\right)}$) components of the spin structure factor in equation (34), where $S_{s}^{\bot}\left(\mathbf{q}\right)=S_{s}^{xx}\left(\mathbf{q}\right)+S_{s}^{yy}\left(\mathbf{q}\right)$. Similarly, figure 10 includes snapshots of the real-space chirality distribution of each phase, along with the square root of chirality structure factors, $\sqrt{S_{\chi}^{u}\left(\mathbf{q}\right)}$ and $\sqrt{S_{\chi}^{d}\left(\mathbf{q}\right)}$, for the up and down triangles of the triangular lattice which are the two sublattices of the dual honeycomb lattice), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At very low fields, the magnetic ordering consists of a double-${{Q}^{\chi}}$ chiral phase with a net chiral component (ferrochiral 3Q^M-2${{Q}^{\chi}}$ spiral). The net chirality has the same sign on the two sublattices of up and down triangles, ${{\chi}^{u}}={{\chi}^{d}}$, which is therefore denoted as ferrochiral. The non-zero chirality does not seem to be consistent with the spin configuration shown in figure 9(a), which looks similar to a vertical 1Q^M spiral ordering with the spiraling moments in the vertical plane. However, upon looking at the spin structure factors shown in figures 9(b) and (c), it becomes clear that this spin ordering contains additional $\mathbf{Q}$ components, which render the spin configuration non-coplanar. This is due to the commensurate relation between ${{\mathbf{Q}}_{\nu}}$ and ${{\mathbf{K}}_{\pm}}$. For example, ${{\mathbf{Q}}_{1}}$ and ${{\mathbf{Q}}_{2}}$ are connected by ${{\mathbf{K}}_{\pm}}$, which means that periodic impurities favor multiple- Q orderings and can extend their stability to the H  =  0 axis.

Upon increasing the magnetic field, the ferrochiral 3Q^M-2${{Q}^{\chi}}$ spiral phase is replaced by a ferrichiral 3Q^M-2${{Q}^{\chi}}$ spiral ordering. By ferrichiral, we mean that the scalar spin chirality has different magnitude for up and down triangles, i.e. ${{\chi}^{u}}\left(\mathbf{q}=0\right)\ne {{\chi}^{d}}\left(\mathbf{q}=0\right)$, as it is shown in figures 10(e) and (f). The other difference relative to the zero field phase is that $S_{s}^{zz}\left(\mathbf{q}\right)$ has two dominant $\mathbf{Q}$ components with the same intensity and a weaker third component, as it is shown in figure 9(f). This explains the anisotropic three-$\mathbf{Q}$ modulation that is observed in the real space spin configuration shown in figure 9(d). Interestingly enough, upon increasing temperature this phase transitions into a different ferrichiral 3Q^M-2${{Q}^{\chi}}$ spiral phase (denoted as ferrichiral 3Q^M-2${{Q}^{\chi}}$ spiral II in figure 8(b)) that exhibits different intensities for the two dominant $\mathbf{Q}$-components of the chirality distribution. A third phase, denoted as ' ferrochiral 3Q^M vortex crystal' ordering in figure 8(b), appears upon further increasing the temperature. This phase is characterized by a net, $\mathbf{q}=0$, component of the scalar spin chirality, ${{\chi}^{u}}\left(\mathbf{q}=0\right)={{\chi}^{d}}\left(\mathbf{q}=0\right)$.

The new phase that appears upon further increasing the magnetic field at low temperatures exhibits a single-$\mathbf{Q}$ modulation of the scalar chirality, as it is clear from the stripy chirality pattern shown in figure 10(g) and denoted as ferrochiral 3Q^M-$1{{Q}^{\chi}}$ spiral. The red/blue stripes correspond to opposite signs of the scalar chirality. A net scalar chirality arises from the fact that the holes occupy sites inside the stripes of a given sign. The suppression of the scalar chirality in their proximity leads to a net chirality of the opposite sign.

Finally, the last phase before the system becomes fully saturated is a ferrochiral 3Q^M-6${{Q}^{\chi}}$ vortex crystal similar to the one reported recently for frustrated quantum magnets [64, 65]. As is it shown in figure 10(j), vortices with the same winding number (±1) or vorticity form a triangular lattice. The lattice of holes is one of the three sublattices of the vortex triangular crystal. Like in the case of [64], vortices with opposite winding number (antivortices) appear at the boundary between two vortices and produce the opposite sign of the scalar spin chirality. This is the blue region in figure 10(j) that appears around the red chiral spots arising from the triangular vortex crystal. As it is clear from the spin and chirality structure factors, the six-fold rotational symmetry of the underlying CDW is restored in this phase. The main difference relative to the high-temperature ferrochiral 3Q^M vortex crystal phase is that six peaks of the chirality structure factor at $\mathbf{q}=\pm {{\mathbf{Q}}_{\nu}}$ in the ferrochiral 3Q^M-6${{Q}^{\chi}}$ vortex crystal exhibit long range ordering. Once again, the ferrochiral 3Q^M-6${{Q}^{\chi}}$ vortex crystal has a net scalar spin chirality because one third of the vortices with a given vorticity (sign of the scalar chirality) have a hole at the center of their cores. These holes suppress the scalar chirality of these vortices (see figure 10(j)) producing an unbalance relative to the contribution of the antivortices.

The corollary of the last two sections is that non-coplanar chiral structures can easily arise from the combination of geometric frustration and inhomogeneous structures. In particular, we have seen that a non-magnetic impurity nucleates a magnetic vortex around it and that a periodic array of non-magnetic 'impurities' (CDW) that is commensurate with the magnetic ordering wave-vectors leads to a rich field- temperature phase diagram with different chiral phases. In particular, as expected from the solution of the single non-magnetic impurity above the bulk saturation field, the highest field broken symmetry state corresponds to a ferrochiral 3Q^M-6${{Q}^{\chi}}$ vortex crystal.

3.2.2. Easy-axis anisotropy and skyrmion crystals.

We have shown in the previous sections that a single impurity can induce a magnetic vortex and that periodic arrays of impurities can stabilize a magnetic vortex crystal. For the triangular lattice under consideration, the vortex crystal corresponds to one type of non-coplanar 3-$\mathbf{Q}$ magnetic ordering. It is then natural to ask if non-coplanar 3-$\mathbf{Q}$ orderings can also be stabilized in absence of impurities or density-waves of a different order parameter that render the medium inhomogeneous.

In general, non-coplanar magnetic orderings are quite scarce among isotropic Mott insulators. In principle, they can appear as multiple-$\mathbf{Q}$ orderings whenever the Fourier transform of the exchange interaction, $J\left(\mathbf{q}\right)$, has global minima at multiple symmetry related wave-vectors ${{\mathbf{Q}}_{\nu}}$. The basic difficulty is that, unlike the case of single-$\mathbf{Q}$ magnetic ordering, the magnitude of the moment is in general modulated for the multiple-$\mathbf{Q}$ orderings. (This is not necessarily true when ${{\mathbf{Q}}_{\nu}}$ are high-symmetry wave-vectors). This modulation is penalized by the longitudinal stiffness of the magnetic moments, which becomes infinitely large in the $T\to 0$ and $S\to \infty$ limit. In other words, the length of the magnetic moments is fixed in this limit, implying that the multiple-$\mathbf{Q}$ orderings must contain higher harmonics, which necessarily increase the exchange interaction energy.

The longitudinal stiffness of the magnetic moments can be reduced by thermal or quantum fluctuations. The most favorable situation is then the proximity to a classical or a quantum critical point (larger thermal or quantum fluctuations suppress the longitudinal stiffness). Indeed, vortex crystal phases have been obtained for S  =  1/2 (quantum limit) spins in the proximity of a magnetic field induced quantum critical point of frustrated antiferromagnets [64, 65]. In addition, Okubo, Chung and Kawamura found a skyrmion crystal phase at finite T and H in the classical frustrated ${{J}_{1}}-{{J}_{3}}$ triangular lattice model of equation (10) discussed in the previous sections [66].

Skyrmion crystals [67, 68] have also been recently discovered in chiral magnets [69–73]. These crystals have triggered enormous interest because of their potential for spintronics applications [74]. Each skyrmion consists of a spin configuration that wraps the sphere once. This property is quantified by the topological skyrmion charge $Q=\frac{1}{4\pi}{{{\int}^{}}_{\text{Sk}}}\text{d}{{r}^{2}}\mathbf{S}\cdot \left({{\partial}_{x}}\mathbf{S}\times {{\partial}_{y}}\mathbf{S}\right)=\pm 1$. In the case of chiral magnets, such as MnSi, the skyrmion crystal phase is induced by competition between FM exchange and a Dzyaloshinskii–Moriya (DM) interaction. This DM interaction breaks the chiral symmetry because it arises from the non-centrosymmetric nature of the B20 crystal structure of these compounds. In contrast, the skyrmion crystal found in [66] arises from frustrated exchange interactions, which preserve the chiral symmetry of the spin Hamiltonian. Consequently, the skyrmion phase of ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$ (10) exhibits spontaneous chiral symmetry breaking.

The skyrmion crystal phase included in the phase diagram of ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$ is also a triple-$\mathbf{Q}$ ordering. This phase appears in the proximity of the thermodynamic phase transition into the paramagnetic state and it is confined to a rather narrow range of magnetic field values [66]. However, unlike the case of chiral magnets, the skyrmion crystal becomes unstable towards single-$\mathbf{Q}$ conical ordering for small enough Q, i.e. in the long wave-length regime [75]. In other words, the skyrmion crystal phase is stable only for $|{{J}_{1}}|\ll {{J}_{3}}$ (see equation (12)).

After these considerations, we can now return to our original question and ask what type of interaction is needed to further stabilize skyrmion crystals in frustrated Mott insulators. We will see here that easy-axis anisotropy is one answer for this question. In particular, we will discuss the interplay between frustration and easy-axis anisotropy and present T  =  0 variational calculations and unbiased finite-T MC simulations showing that skyrmion crystals are induced for large enough easy-axis anisotropy irrespective of the magnitude of the ordering wave vector. Recently, Leonov and Mostovoy [76] presented an alternative variational T  =  0 study of the effect of easy-axis anisotropy on a different frustrated model for classical spins. They also found that easy-axis anisotropy helps to stabilize non-coplanar spin orderings. Indeed, [76] provides a simple argument for understanding the instability of the single-$\mathbf{Q}$ conical phase towards multiple-$\mathbf{Q}$ ordering in the presence of magnetic field and a single-ion easy-axis anisotropy term,

$$\begin{eqnarray}&&{{\mathcal{H}}_{A}}=-A\underset{j}{\sum}{{\left(S_{j}^{z}\right)}^{2}},\end{eqnarray} \tag{ 38 }$$

with A  >  0. We start by considering the conical state:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathbf{S}}_{\mathbf{r}}} & = & \sqrt{1-m_{z}^{2}}\,\left[\cos \left({{\mathbf{Q}}_{1}}\cdot \mathbf{r}\right)\hat{\mathbf{x}}+\sin \left({{\mathbf{Q}}_{1}}\cdot \mathbf{r}\right)\hat{\mathbf{y}}\right]+{{m}_{z}}\hat{\mathbf{z}}, \end{array}\end{eqnarray} \tag{ 39 }$$

which is the ground state in absence of anisotropy (see figure 8(a)). m_z is the uniform spin magnetization. We will now consider the following deformation of the conical phase¹,

$$\begin{eqnarray}\begin{array}{*{35}{l}} S_{j}^{x} & = & \sqrt{{{\sin}^{2}}\tilde{\theta}-\Delta _{2}^{2}}\cos \left({{\mathbf{Q}}_{1}}\cdot {{\mathbf{r}}_{j}}\right)+{{\Delta}_{2}}\cos \left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{j}}\right), \\ S_{j}^{y} & = & \sqrt{{{\sin}^{2}}\tilde{\theta}-\Delta _{2}^{2}}\sin \left({{\mathbf{Q}}_{1}}\cdot {{\mathbf{r}}_{j}}\right)-{{\Delta}_{2}}\sin \left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{j}}\right), \\ S_{j}^{z} & = & \sqrt{{{\cos}^{2}}\tilde{\theta}-2{{\Delta}_{2}}\sqrt{{{\sin}^{2}}\tilde{\theta}-\Delta _{2}^{2}}\cos {{\mathbf{Q}}_{3}}\cdot {{\mathbf{r}}_{j}}}, \end{array}\end{eqnarray} \tag{ 40 }$$

where the amplitude of the ${{\mathbf{Q}}_{2}}$ component, ${{\Delta}_{2}}$, is a variational parameter and $\cos \theta =\cos \tilde{{\theta}}[1-\Delta {{_{2}^{2}}_{{}}}({{\sin}^{2}}\tilde{{\theta}}-\Delta _{2}^{2})/(4{{\cos}^{4}}\tilde{{\theta}})]$. A key observation is that the z-component is modulated by the third wave-vector ${{\mathbf{Q}}_{3}}$ to linear order in ${{\Delta}_{2}}$ to preserve the constraint $\mathbf{S}_{i}^{2}=1$. This is so because of the C₆ symmetry of the TL, which implies that ${{\mathbf{Q}}_{1}}+{{\mathbf{Q}}_{2}}+{{\mathbf{Q}}_{3}}=0$. Therefore, the exchange energy cost of the higher harmonics produced by the normalization condition is proportional to $\Delta _{2}^{4}[J(2{{\mathbf{Q}}_{3}})-J(\mathbf{0})]$, while the anisotropy energy gain produced by the same modulation is proportional to $-A\Delta _{2}^{2}$. This implies that ${{\Delta}_{2}}\propto \sqrt{A/[J(2{{\mathbf{Q}}_{3}})-J(\mathbf{0})]}$ for ${{\Delta}_{2}}\ll 1$ or $A\ll |J(2{{\mathbf{Q}}_{\nu}})|$ [75]. This simple analysis shows how the presence of easy-axis anisotropy can induce multiple-$\mathbf{Q}$ magnetic orderings. It is interesting to note that a double-$\mathbf{Q}$ conical state, like the one described by equation (40), has been found in a spatially anisotropic triangular lattice model near the saturation field [77]. In the rest of this section we will study the field- temperature phase diagram of ${{\mathcal{H}}_{{{J}_{1}}-{{J}_{3}}}}$ with a relatively large value (A  =  0.5 J₁) of the single-ion anisotropy term given in equation (38). A more extensive study is presented in [75].

Figure 11(a) shows a typical phase diagram in absence of single-ion anisotropy for a small Q value ($Q=2\pi /5$, J₁  =  −1.0, J₃  =  0.5). The phase diagram only includes two phases: the paramagnetic state that is obtained at high enough temperature or magnetic field and the single-$\mathbf{Q}$ conical phase that we found in the previous section. Figure 11(b) shows the phase diagram in presence of easy-axis anisotropy (A  =  0.5) for the same parameters of figure 11(a) and for a magnetic field parallel to the easy-axis. A new skyrmion crystal phase divides the helical spin ordering from the paramagnetic state in the high-field region.

Zoom In Zoom Out Reset image size

Figures 12(c) and 13(a) show typical real-space configurations for spin and scalar chirality in the skyrmion phase. The blue regions correspond to the skyrmion cores where the spins are anti-aligned with the magnetic field. As expected for a 3-$\mathbf{Q}$ ordering, the skyrmions form a triangular lattice with lattice parameter $4\pi /\left(\sqrt{3}Q\right)$. The skyrmion crystal naturally produces a net scalar chirality (each skyrmion subtends a solid angle of 4π). The scalar chirality exhibits the same modulation as the skyrmion lattice because the chiral density is maximized at the skyrmion cores. It is important to note that the net scalar chirality can be positive or negative because the Hamiltonian is invariant under the Z₂ mirror symmetry that changes the sign of χ. As we already explained, this situation is qualitatively different from the case of chiral magnets in which the skyrmion crystal is induced by the competition between FM exchange and a chiral DM interaction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 shows snapshots of the spin distribution in the helical and the skyrmion phases. The spin component parallel to the external field is indicated with a contour plot that reveals the characteristic stripes that arise from helical ordering and the triangular lattice of skyrmion cores that characterize the skyrmion crystal. The multiple-$\mathbf{Q}$ nature of the skyrmion crystal becomes evident when looking at the spin structure factors of the helical ordering (figure 12(b)) and the skyrmion crystal (figure 12(d)). While the helical ordering exhibits two Bragg peaks at  ±${{\mathbf{Q}}_{1}}$, the skyrmion crystal has six equal intensity peaks at  ±${{\mathbf{Q}}_{\nu}}$, with $\nu =1,2,3$. The $\mathbf{q}=0$ peak, which reflects the uniform magnetization component induced by the external magnetic field H, has been subtracted to emphasize the peaks that arise from spontaneous symmetry breaking. Given that the xy spin components can only be ordered at T  =  0 [78], only the z-component of the structure factor, ${{S}^{zz}}\left(\mathbf{q}\right)$, can exhibit Bragg peaks at finite-T.

Another qualitative difference between the helical ordering and the skyrmion crystal is that the former is coplanar, while the latter is non-coplanar. As we already discussed, the local order parameter that distinguishes these two situations is the local scalar spin chirality ${{\chi}_{jkl}}$. Figure 13(a) shows a snapshot of the distribution of ${{\chi}_{jkl}}$ in the skyrmion phase. Besides the expected triple-$\mathbf{Q}$ modulation, there is a uniform component, which leads to a large $\mathbf{q}\mathbf{=}\mathbf{0}$ peak in the chiral structure factor shown in figure 13(b). This uniform chiral component arises from the fact that all the skyrmions have the same skyrmion number (1 or  −1). (Note that the skyrmion number is equal to the integral of ${{\chi}_{jkl}}$ in the continuum limit.)

Our MC results suggest that the skyrmion crystal survives all the way to T  =  0. Based on our previous arguments, this would be not possible for isotropic classical spin models because the higher harmonics required to keep the moments normalized in a multiple-$\mathbf{Q}$ state are penalized by the isotropic exchange interaction. However, as we have seen above, an easy-axis spin anisotropy favors the generation of higher harmonics in an attempt of aligning the magnetic moments with the easy-axis and it opens the possibility of stabilizing the skyrmion phase down to T  =  0.

The complement the MC results, we introduce a T  =  0 variational approach. In particular, we choose the following variational state to describe the skyrmion crystal phase [66]:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathbf{S}_{\mathbf{r}}^{xy} & = & \frac{{{I}_{xy}}}{\mathcal{N}_{\mathbf{r}}^{s}}\underset{\nu =1,3}{\sum}\sin \left({{\mathbf{Q}}_{\nu}}\cdot \mathbf{r}+{{\theta}_{\nu}}\right){{\mathbf{e}}_{\nu}} \\ S_{\mathbf{r}}^{z} & = & {{m}_{z}}-\frac{{{I}_{z}}}{\mathcal{N}_{\mathbf{r}}^{s}}\underset{\nu =1,3}{\sum}\cos \left({{\mathbf{Q}}_{\nu}}\cdot \mathbf{r}+{{\theta}_{\nu}}\right){{\mathbf{e}}_{\nu}} \end{array}\end{eqnarray} \tag{ 41 }$$

where m_z is the uniform spin magnetization and $\mathcal{N}_{\mathbf{r}}^{s}$ is an $\mathbf{r}$-dependent normalization factor with enforces the $\mathbf{S}_{\mathbf{r}}^{2}=1$ constraint. The amplitudes m_z, I_z and I_xy are treated as variational parameters. Given that the single-$\mathbf{Q}$ conical phase (39) is unstable, we will only consider the helical phase (proper-screw spiral),

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathbf{S}_{\mathbf{r}}^{xy} & = & \frac{{{I}_{xy}}}{\mathcal{N}_{\mathbf{r}}^{h}}\cos \left({{\mathbf{Q}}_{\nu}}\cdot \mathbf{r}\right)\hat{\mathbf{x}}, \\ S_{\mathbf{r}}^{z} & = & {{m}_{z}}+\frac{{{I}_{z}}}{\mathcal{N}_{\mathbf{r}}^{h}}\sin \left({{\mathbf{Q}}_{\nu}}\cdot \mathbf{r}\right), \end{array}\end{eqnarray} \tag{ 42 }$$

as an alternative variational state. $\mathcal{N}_{\mathbf{r}}^{h}$ is the normalization factor of the helical ordering.

The results of the variational approach are also shown with red triangular dots in the phase diagram of figure 11(b). These dots are quite close to the $T\to 0$ extrapolation of the phase boundaries obtained from the unbiased MC simulations, confirming the fact that the skyrmion phase extends all the way to T  =  0.

The qualitative differences between the phase diagrams with and without easy-axis anisotropy (see figure 11) can be understood in simple terms. The $\mathbf{r}$-dependence of $\mathcal{N}_{\mathbf{r}}^{s}$ and $\mathcal{N}_{\mathbf{r}}^{h}$ generates higher harmonic components, which are penalized by the isotropic exchange. This is the reason why the conical phase prevails for A  =  0. However, this energy cost from the exchange interaction can be compensated by an energy gain from the single-ion anisotropy. The tendency to align the spins along the easy-axis direction replaces the sinusoidal modulation by a soliton-like modulation, which includes higher harmonic components. From the variational ansatz (41) and (42), it is clear that the helical state and the skyrmion crystal can be continuously distorted into a soliton like modulation, while the same is not true for conical phase. It is not surprising, then, that the single conical phase for A  =  0 is replaced by the helical state and the skyrmion crystal for strong enough easy-axis anisotropy.

When the anisotropy A becomes a significant fraction of $|{{J}_{1}}|$, as in the case under consideration (A  =  0.5 J₁), the transition between the spin down to the spin up regions occurs over length scale of order $\sqrt{J/A}$. This length scale can be made much shorter than $2\pi /|Q|$ in the long wave length limit $|Q|\ll 1$. We can then assume that the boundary between regions with opposite spin alignment is a line with positive tension, implying that the energy of a given state can be reduced by minimizing the perimeter of the boundary per unit of area. In the case of the helical state, the perimeter per unit of area, ${{P}_{h}}/{{A}_{h}}=Q/\pi$, does not depend on the value of the uniform magnetization m_z. This is so because the effect of applying a magnetic field to the helical state is to translate up-down boundaries to the right and the down-up boundaries to the left in order to shrink (expand) the spin down (up) stripes. In contrast, the perimeter per unit area for the skyrmion crystal, ${{P}_{s}}/{{A}_{s}}={{3}^{1/4}}Q\sqrt{1-{{m}_{z}}}/2\sqrt{\pi}$, does depend on m_z because the skyrmion cores shrink when increasing the magnetic field. We then expect a transition from the helical to the skyrmion crystal state when ${{P}_{h}}/{{A}_{h}}$ becomes approximately equal to ${{P}_{s}}/{{A}_{s}}$. This condition leads to a critical value of ${{m}_{z}}=m_{z}^{c}\simeq 0.265$. At this point it is important to note that the transition between both phases is of first order, i.e. the magnetization changes discontinuously, as it is shown in figure 14. The average between the magnetization values right below and above the transition is approximately 0.24, which is in good agreement with our simple estimate.

Zoom In Zoom Out Reset image size

The conclusion of this section is that easy-axis anisotropy can easily induce a skyrmion crystal right below the saturation field of a frustrated magnet. Unlike the case of the isotropic model [66], the skyrmion phase stabilized by easy-axis anisotropy survives in the long wave length limit $|{{\mathbf{Q}}_{\nu}}|\ll 1$. Other multiple-$\mathbf{Q}$ orderings can also be obtained for smaller values of the anisotropy [75, 76]. The following three ingredients are enough to obtain field-induced multiple-$\mathbf{Q}$ ordering: (1) C₆ symmetry, such as triangular and kagome lattices, (2) finite $|\mathbf{Q}|$ ordering due to competing interactions, and (3) easy-axis anisotropy. Fe_xNi_1−xBr₂ [79], Zn_xNi_1−xBr₂ [48], and an Fe monolayer on Ir(1 1 1) [80, 81] may be candidate materials to exhibit field-induced skyrmion crystal phase.

So far we have considered the effects of frustration in MIs, i.e. systems whose magnetic degrees of freedom are always localized. It is also interesting to consider the case of itinerant magnets. In this case we need to distinguish between two quite different situations: systems whose magnetic moments only arise from the spin or orbital degrees of freedom of the conduction electrons and systems that contain localized moments in addition to the magnetic degrees of freedom of the conduction electrons. While it is interesting to analyze the role of frustration in both systems, here we will mainly focus on the latter case, which describes the situation of many 4f and 3d-electron materials.

The small size of 4f-orbitals leads to the formation of localized magnetic moments when the energy of the f-orbitals is well below the Fermi level of the conduction electrons. Like in the case of MIs, double occupancy of the f-orbitals is suppressed by an intra-orbital Coulomb interaction, which is much larger than the hybridization between the f-orbital and the orbtials that generate the broad conduction band. The minimal model that describes this situation is the Kondo Lattice (KL) Hamiltonian:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{H}}_{\text{KL}}}=\underset{\mathbf{k},\sigma}{\sum}{{\varepsilon}_{\mathbf{k}}}c_{\mathbf{k}\sigma}^{\dagger}{{c}_{\mathbf{k}\sigma}}-J\underset{\mathbf{k},\mathbf{q},\sigma,{{\sigma}^{\prime}}}{\sum}c_{\mathbf{k}\sigma}^{\dagger}{{\sigma}_{\sigma {{\sigma}^{\prime}}}}{{c}_{\mathbf{k}+\mathbf{q}{{\sigma}^{\prime}}}}\cdot {{\mathbf{S}}_{\mathbf{q}}}. \end{array}\end{eqnarray} \tag{ 43 }$$

The operator $c_{\mathbf{k}\sigma}^{\dagger}\left({{c}_{\mathbf{k}\sigma}}\right)$ creates (annihilates) an itinerant electron with momentum $\mathbf{k}$ and spin σ. The bare dispersion relation of the conduction electrons is ${{\varepsilon}_{\mathbf{k}}}=\sum\nolimits_{jl}{{t}_{jl}}{{\text{e}}^{\text{i}\mathbf{k}\cdot {{\mathbf{r}}_{jl}}}}$, where t_jl are hopping integrals between sites j and l. The second term of ${{\mathcal{H}}_{\text{KL}}}$ corresponds to the exchange coupling (J) between the conduction electrons and the localized magnetic moments in the momentum space representation ${{\mathbf{S}}_{\mathbf{q}}}=\sum\nolimits_{\mathbf{r}}{{\text{e}}^{\text{i}\mathbf{q}\cdot \mathbf{r}}}{{\mathbf{S}}_{\mathbf{r}}}/\sqrt{N}$ (N is the number of lattice sites). For the case of 4f-electron systems, J corresponds to the so-called Kondo AFM interaction. This model also applies to 3d-electron transition metal oxides, like the manganites [82], in which the localized moments are provided by t_2g electrons, while the itinerant electrons occupy a band of e_g-orbitals. In this case, J corresponds to the FM Hund's coupling between electrons in the same atomic shell. The energy gap separating the localized t_2g-orbitals from the itinerant e_g-orbitals is provided by the strong crystal field produced by the oxygen octahedron that surrounds the transition metal.

Here we will only consider the case of classical magnetic moments ($|{{\mathbf{S}}_{i}}|=1$), i.e. the $S\to \infty$ limit. Although the Kondo effect disappears in this limit, we will still refer to ${{\mathcal{H}}_{\text{KL}}}$ as the KL Hamiltonian. The sign of J becomes irrelevant for classical moments because it can be absorbed by the transformation ${{\mathbf{S}}_{\mathbf{q}}}\to -{{\mathbf{S}}_{\mathbf{q}}}$. Therefore, without loss of generality, we will assume that the coupling is FM.

In the weak-coupling regime, $J\rho \left({{\varepsilon}_{\text{F}}}\right)\ll 1$, where $\rho \left({{\varepsilon}_{\text{F}}}\right)$ is the density of states at the Fermi level ${{\varepsilon}_{\text{F}}}$, it is possible to integrate out the conduction electrons to produce an effective spin Hamiltonian up to second order in J. This procedure leads to the RKKY Hamiltonian

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{H}_{\text{eff}}^{(2)}=-{{J}^{2}}\underset{\mathbf{q}}{\sum}{\chi}_{\mathbf{q}}^{0}{{\mathbf{S}}_{\mathbf{q}}}\cdot {{\mathbf{S}}_{-\mathbf{q}}}. \end{array}\end{eqnarray} \tag{ 44 }$$

Here, ${\chi}_{\mathbf{q}}^{0}=T\sum\nolimits_{\mathbf{k},{{\omega}_{n}}}G_{\mathbf{k},{{\omega}_{n}}}^{0}G_{\mathbf{k}+\mathbf{q},{{\omega}_{n}}}^{0}$ is the bare susceptibility of the conduction electrons, ${{\omega}_{n}}=(2n+1)\pi T$ are the Matsubara frequencies, and $G_{\mathbf{k},{{\omega}_{n}}}^{0}={{\left\{i{{\omega}_{n}}-\left[{{\varepsilon}_{\mathbf{k}}}-\mu \right]\right\}}^{-1}}$ is the bare Green function (μ is the chemical potential). We note that ${{\mathbf{S}}_{\mathbf{q}}}\cdot {{\mathbf{S}}_{-\mathbf{q}}}\geqslant 0$ and $\sum\nolimits_{\mathbf{q}}{{\mathbf{S}}_{\mathbf{q}}}\cdot {{\mathbf{S}}_{-\mathbf{q}}}=N$ because of the local constraint $|{{\mathbf{S}}_{i}}|\,\,=1$. According to equation (44), any single-$\mathbf{Q}$ helical state whose ordering vector maximizes ${\chi}_{\mathbf{q}}^{0}$ is a ground state of $\mathcal{H}_{\text{eff}}^{(2)}$. This is the origin of the helical spin arrangements in rare-earth and other itinerant magnets that were identified towards the end of fifties and the beginning of sixties [8–10].

When ${\chi}_{\mathbf{q}}^{0}$ has a single global maximum, $\mathcal{H}_{\text{eff}}^{(2)}$ has a unique ground state up to global spin rotations. In contrast, the existence of multiple wave-vectors ${{\mathbf{Q}}_{\nu}}$ that maximize ${\chi}_{\mathbf{q}}^{0}$ implies that $\mathcal{H}_{\text{eff}}^{(2)}$ has multiple ground states. In this case, we will say that the itinerant magnet is frustrated because there are many competing orderings. For instance, the circular Fermi surface (FS) of a 2D electron gas produces a highly-frustrated magnetic interaction: ${{\chi}^{0}}\left(\mathbf{q}\right)$ has a flat global maximum for $0\leqslant q\leqslant 2{{k}_{\text{F}}}$ and it decreases for $q>2{{k}_{\text{F}}}$ (see figure 15(a)), implying that any $\mathbf{q}$-spiral ordering with $q\leqslant 2{{k}_{\text{F}}}$ is a ground state of $\mathcal{H}_{\text{eff}}^{(2)}$.

Zoom In Zoom Out Reset image size

The ground state degeneracy of $\mathcal{H}_{\text{eff}}^{(2)}$ is removed by contributions of order J ⁴ and higher, which were neglected in going from ${{\mathcal{H}}_{\text{KLM}}}$ to $\mathcal{H}_{\text{eff}}^{(2)}$. As we will see below, these contributions often select exotic magnetic orderings which are qualitatively different from a simple single- Q spirals. More importantly, the resulting multiple-$\mathbf{Q}$ orderings are in general non-coplanar, i.e. they have a non-zero scalar spin chirality: $\langle {{\chi}_{jkl}}\rangle \ne 0$. This order parameter couples to the orbital degrees of freedom of the conduction electrons. The exchange interaction, J, tends to align the spins of the conduction electrons to the underlying local moment texture, inducing a Berry phase in the electron wave-function, which is indistinguishable from a 'real' magnetic flux (see figure 16). When conduction electrons propagate through such a spin texture, they may exhibit spontaneous Hall effect in the absence of any externally applied magnetic field [84]. Moreover, a spontaneous quantum Hall insulator can emerge from the chiral ordering if it gaps out the full Fermi surface [85–87].

Zoom In Zoom Out Reset image size

Here we will follow [88] to illustrate the role of frustration in the instability of single-$\mathbf{Q}$ spirals toward multiple-$\mathbf{Q}$ non-coplanar orderings. We start by considering the simple case of a square KL Hamiltonian (43), whose magnetic susceptibility has four global maxima connected by the C₄ symmetry operations of the square lattice. This is a quite general situation for arbitrary band filling fractions, which typically leads to susceptibility maxima at low-symmetry wave-vectors. For the case under consideration, these wave-vectors are $\pm {{\mathbf{Q}}_{\nu}}$, with $\nu =1,2$. We will also assume that these ordering wave-vectors connect opposite points on the FS.

We argue that the relevant variational ansatz is a double-Q vortex crystal state,

$$\begin{eqnarray}&&\mathbf{S}_{i}^{\text{vc}}(b)=\left(\begin{array}{c} \sqrt{\left(1-{{b}^{2}}\right)+{{b}^{2}}{{\cos}^{2}}\left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right)}\cos \left({{\mathbf{Q}}_{1}}\cdot {{\mathbf{r}}_{i}}\right) \\ \sqrt{\left(1-{{b}^{2}}\right)+{{b}^{2}}{{\cos}^{2}}\left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right)}\sin \left({{\mathbf{Q}}_{1}}\cdot {{\mathbf{r}}_{i}}\right) \\ b\sin \left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right) \end{array}\right),\end{eqnarray} \tag{ 45 }$$

where ${{\mathbf{r}}_{i}}$ is the position vector of the site i. For b  =  0, equation (45) describes an ordinary single-${{\mathbf{Q}}_{1}}$ spiral in the xy spin-plane. A nonzero b introduces an out-of-plane (z component) modulation of amplitude b and wave vector ${{\mathbf{Q}}_{2}}$, which deforms the single-${{\mathbf{Q}}_{1}}$ spiral into a vortex crystal state, as it is shown in figure 17. This state also exhibits a one-dimensional modulation of the scalar chirality (chirality stripes) along the ${{\mathbf{Q}}_{2}}$ direction.

Zoom In Zoom Out Reset image size

The basic motivation for proposing this particular multiple-$\mathbf{Q}$ ordering is that the vector amplitude of the ${{\mathbf{Q}}_{2}}$ component is orthogonal to the vector amplitude of the ${{\mathbf{Q}}_{1}}$-spiral. This condition guarantees that the amplitude of the higher harmonics required to enforce the constraint $\mathbf{S}_{i}^{2}=1$ is of order b², instead of the order b that would be obtained if the two vector amplitudes were not orthogonal to each other. As a consequence, this particular choice of multiple-$\mathbf{Q}$ ordering has an RKKY energy cost of order ${{J}^{2}}{{b}^{4}}$ instead of $\mathcal{O}\left({{J}^{2}}{{b}^{2}}\right)$.

The spin configuration of the vortex crystal in equation (45) has higher harmonics, which arise from an expansion of the square root prefactors introduced to satisfy the normalization condition $\mathbf{S}_{i}^{2}=1$. This is true in general for any multiple-$\mathbf{Q}$ solution, as long as ${{\mathbf{Q}}_{1}}$ and ${{\mathbf{Q}}_{2}}$ are not high-symmetry wave vectors. This observation implies that multiple-$\mathbf{Q}$ orderings have higher energy than the single-Q helical ordering at the RKKY level. However, as we will see below, the energy balance can change when higher-order terms are included.

The fourth-order contribution to the perturbative expansion in powers of J gives

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{H}_{\text{eff}}^{(4)}= & -{{J}^{4}}\underset{{{\mathbf{q}}_{j}},j=1,4}{\sum}{{\Pi}_{{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}},{{\mathbf{q}}_{3}},{{\mathbf{q}}_{4}}}}\left({{\mathbf{S}}_{{{\mathbf{q}}_{1}}}}\cdot {{\mathbf{S}}_{{{\mathbf{q}}_{2}}}}\right)\left({{\mathbf{S}}_{{{\mathbf{q}}_{3}}}}\cdot {{\mathbf{S}}_{{{\mathbf{q}}_{4}}}}\right) \\ {} & \times \,\delta \left({{\mathbf{q}}_{1}}+{{\mathbf{q}}_{2}}+{{\mathbf{q}}_{3}}+{{\mathbf{q}}_{4}}-\mathbf{G}\right), \end{array}\end{eqnarray} \tag{ 46 }$$

where ${{\Pi}_{{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}},{{\mathbf{q}}_{3}},{{\mathbf{q}}_{4}}}}$ is the coefficient for the four-spin scattering process with momenta ${{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}},{{\mathbf{q}}_{3}}$, and ${{\mathbf{q}}_{4}}$ and $\mathbf{G}$ is a reciprocal lattice vector. Because we are interested in scattering processes involving the lowest energy modes ${{q}_{1\leqslant j\leqslant 4}}\in \left\{\pm {{\mathbf{Q}}_{\nu}}\right\}$, we only need to keep the contributions ${{A}_{\mathbf{q}}}={{\Pi}_{\mathbf{q},-\mathbf{q},\mathbf{q},-\mathbf{q}}}$, ${{B}_{\mathbf{q},{{\mathbf{q}}^{\prime}}}}={{\Pi}_{\mathbf{q},-\mathbf{q},{{\mathbf{q}}^{\prime}},-{{\mathbf{q}}^{\prime}}}}$, and ${{W}_{\mathbf{q},{{\mathbf{q}}^{\prime}}}}={{\Pi}_{\mathbf{q},{{\mathbf{q}}^{\prime}},-\mathbf{q},-{{\mathbf{q}}^{\prime}}}}$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{A}_{\boldsymbol{q}}} & = & T\underset{\boldsymbol{k}\in 1\text{BZ},{{\omega}_{n}}}{\sum}{{(G_{\boldsymbol{k},{{\omega}_{n}}}^{0})}^{2}}{{(G_{\boldsymbol{k}+\boldsymbol{q},{{\omega}_{n}}}^{0})}^{2}}, \\ {{B}_{\boldsymbol{q},\boldsymbol{q}^{{\prime}}}} & = & T\underset{\boldsymbol{k}\in 1\text{BZ},{{\omega}_{n}}}{\sum}G_{\boldsymbol{k},{{\omega}_{n}}}^{0}{{(G_{\boldsymbol{k}+\boldsymbol{q},{{\omega}_{n}}}^{0})}^{2}}G_{\boldsymbol{k}+\boldsymbol{q}+\boldsymbol{q}^{{\prime}},{{\omega}_{n}}}^{0}, \\ {{W}_{\boldsymbol{q},\boldsymbol{q}^{{\prime}}}} & = & T\underset{\boldsymbol{k}\in 1\text{BZ},{{\omega}_{n}}}{\sum}G_{\boldsymbol{k},{{\omega}_{n}}}^{0}G_{\boldsymbol{k}+\boldsymbol{q},{{\omega}_{n}}}^{0}G_{\boldsymbol{k}+\boldsymbol{q}^{{\prime}},{{\omega}_{n}}}^{0}G_{\boldsymbol{k}+\boldsymbol{q}+\boldsymbol{q}^{{\prime}},{{\omega}_{n}}}^{0}, \end{array}\end{eqnarray} \tag{ 47 }$$

which correspond to the Feynman diagrams depicted in figure 18.

Zoom In Zoom Out Reset image size

In the following we will assume that $b\ll 1$. By Fourier transforming the spin configuration in equation (45) and substituting into equation (46), we obtain the energy up to $\mathcal{O}({{J}^{4}})$:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}E_{\text{eff}}^{(4)}(b)=-{{J}^{2}}\left[\left(1-\frac{{{b}^{4}}}{32}\right){\chi}_{{{\mathbf{Q}}_{1}}}^{0}+\frac{{{b}^{4}}}{32}{\chi}_{{{\mathbf{Q}}_{1}}+2{{\mathbf{Q}}_{2}}}^{0}\right]\\ -\frac{{{J}^{4}}}{2}\left[\left(1-{{b}^{2}}\right){{A}_{{{\mathbf{Q}}_{1}}}}+{{b}^{2}}{{B}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}+{{b}^{2}}{{W}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}\right].\end{array}\end{eqnarray} \tag{ 48 }$$

The energy difference between the proposed multiple-$\mathbf{Q}$ ordering and the single-$\mathbf{Q}$ spiral, $\Delta E_{\text{eff}}^{(4)}(b)\equiv E_{\text{eff}}^{(4)}(b)-E_{\text{eff}}^{(4)}(b=0)$, is

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\Delta E_{\text{eff}}^{(4)}(b)=\frac{{{J}^{2}}{{b}^{4}}}{32}\left[{\chi}_{{{\mathbf{Q}}_{1}}}^{0}-{\chi}_{{{\mathbf{Q}}_{1}}+2{{\mathbf{Q}}_{2}}}^{0}\right]\\ -\frac{{{J}^{4}}{{b}^{2}}}{2}\left[{{A}_{{{\mathbf{Q}}_{1}}}}+{{B}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}+{{W}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}\right]\equiv \alpha {{J}^{2}}{{b}^{4}}-\beta {{J}^{4}}{{b}^{2}}.\end{array}\end{eqnarray} \tag{ 49 }$$

The coefficient α is always positive because ${\chi}_{\mathbf{q}}^{0}$ is maximized at $\mathbf{q}=\pm {{\mathbf{Q}}_{1}},\pm {{\mathbf{Q}}_{2}}$. If β also turns out to be positive, the optimal value of b is ${{b}_{\text{opt}}}=\sqrt{\beta /\left(2\alpha \right)}J$. Indeed, explicit evaluation of the coefficients ${{A}_{{{\mathbf{Q}}_{1}}}}$, ${{B}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}$, and ${{W}_{{{\mathbf{Q}}_{2}},{{\mathbf{Q}}_{2}}}}$ gives $\beta >0$ at a low enough T, as long as ${{\mathbf{Q}}_{\nu}}$ connects opposite points on the FS. Moreover, ${{A}_{{{\mathbf{Q}}_{1}}}}$, ${{B}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}$, and ${{W}_{{{\mathbf{Q}}_{1}},{{\mathbf{Q}}_{2}}}}$ diverge for $T\to 0$. While this observation implies that the perturbative expansion is not valid in the $T\to 0$ limit, it suggests that b can become of order one even for very small J.

This hypothesis can be verified by introducing a regular (nondivergent) perturbative treatment. This approach requires to change the local reference frame for the spin components of the itinerant electrons. The local quantization x-axis is now aligned with the local moments of the single-${{\mathbf{Q}}_{1}}$ helical state. The local moments of the double-Q non-coplanar vortex crystal become

$$\begin{eqnarray}\begin{array}{*{20}{l}}{{\mathbf{S}}_{i}}\to {{\tilde{{\mathbf{S}}}}_{i}}=&\left(\begin{array}{c} \sqrt{1-{{b}^{2}}+{{b}^{2}}{{\cos}^{2}}\left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right)} \\ 0 \\ b\sin \left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right) \end{array}\right)\\ &\sim \left(\begin{array}{c} 1-\frac{1}{4}{{b}^{2}}+\frac{1}{4}{{b}^{2}}\cos \left(2{{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right) \\ 0 \\ b\sin \left({{\mathbf{Q}}_{2}}\cdot {{\mathbf{r}}_{i}}\right) \end{array}\right)\end{array}\end{eqnarray} \tag{ 50 }$$

in the new reference frame. By applying the same transformation to the creation and annihilation operators of conduction electrons, we obtain:

$$\begin{eqnarray}\begin{array}{*{35}{l}} c_{\mathbf{k}\sigma}^{\dagger} & \to \tilde{c}_{\mathbf{k}\sigma}^{\dagger}=c_{\left[\mathbf{k}+\sigma {{\mathbf{Q}}_{1}}/2\right]\sigma}^{\dagger}, \end{array}\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{c}_{\mathbf{k}\sigma}} & \to {{{\tilde{c}}}_{\mathbf{k}\sigma}}={{c}_{\left[\mathbf{k}+\sigma {{\mathbf{Q}}_{1}}/2\right]\sigma}}. \end{array}\end{eqnarray} \tag{ 52 }$$

Up to the quadratic order in b (we assume that $b\ll 1$), we obtain the following expression for $\mathcal{H}$ with the local moment configuration given by equation (45):

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{H}_{\text{eff}}^{\text{vc}}(b)= & \underset{\mathbf{k},\sigma}{\sum}{{\tilde{{\varepsilon}}}_{\mathbf{k},\sigma}}\tilde{{c}}_{\mathbf{k}\sigma}^{\dagger}{{\tilde{{c}}}_{\mathbf{k}\sigma}}-J\left(1-\frac{{{b}^{2}}}{4}\right)\tilde{{c}}_{\mathbf{k}\sigma}^{\dagger}{{\tilde{{c}}}_{\mathbf{k}-\sigma}}. \\ {} & -\,\text{i}\frac{Jb}{2}\left(\tilde{{c}}_{\mathbf{k}+{{\mathbf{Q}}_{2}}\sigma}^{\dagger}{{\tilde{{c}}}_{\mathbf{k}\sigma}}-\tilde{{c}}_{\mathbf{k}\sigma}^{\dagger}{{\tilde{{c}}}_{\mathbf{k}+{{\mathbf{Q}}_{2}}\sigma}}\right) \\ {} & +\,\frac{J{{b}^{2}}}{8}\left(\tilde{{c}}_{\mathbf{k}+2{{\mathbf{Q}}_{2}}\sigma}^{\dagger}{{\tilde{{c}}}_{\mathbf{k}\mathbf{-}\sigma}}+\tilde{{c}}_{\mathbf{k}\sigma}^{\dagger}{{\tilde{{c}}}_{\mathbf{k}+2{{\mathbf{Q}}_{2}}-\sigma}}\right), \end{array}\end{eqnarray} \tag{ 53 }$$

with ${{\tilde{\varepsilon}}_{\mathbf{k},\sigma}}={{\varepsilon}_{\mathbf{k}+\sigma {{\mathbf{Q}}_{1}}/2}}$. This is a more adequate basis for applying perturbation theory because electronic scattering with the ${{\mathbf{Q}}_{1}}$ component of the spin configuration is accounted for exactly via direct diagonalization of the first two terms of $\mathcal{H}_{\text{eff}}^{\text{vc}}(b)$. The divergence of the fourth-order contributions to the perturbative treatment in the original basis is avoided in this new basis because they are incorporated in second order contributions in the amplitude of the last two terms of $\mathcal{H}_{\text{eff}}^{\text{vc}}(b)$ [88].

The most significant second-order contributions to the total energy arise from the regions around the two pairs of FS points (hot spots) connected by the ordering vectors ${{\mathbf{Q}}_{1}}$ and ${{\mathbf{Q}}_{2}}$ (see figures 19(a) and (b)). The dispersion around these points can be parametrized in terms of the Fermi velocity, v_F, and the radius of curvature R₀ of the FS. Here we are assuming that ${{R}_{0}}\gg 1$, i.e. that the Fermi segments connected by ${{\mathbf{Q}}_{\nu}}$ are close to be nested.

Zoom In Zoom Out Reset image size

The change in energy induced by a nonzero b value (out-of-plane component) has a positive contribution (energy cost), $\Delta {{E}_{1}}$, arising from the first two terms of equation (53) and a negative contribution (energy gain), $\Delta {{E}_{2}}$, that arises from the last two terms of the same equation. By following the procedure described in [88], we obtain

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Delta {{E}_{1}} & =-\frac{\Delta k}{4{{\pi}^{2}}}\frac{{{J}^{2}}{{b}^{2}}}{16{{v}_{\text{F}}}}\ln \left[\frac{x+\sqrt{{{(\frac{J}{2{{v}_{\text{F}}}})}^{2}}+{{x}^{2}}}}{{{k}_{0}}+x+\sqrt{{{(\frac{J}{2{{v}_{\text{F}}}})}^{2}}+{{({{k}_{0}}+x)}^{2}}}}\right], \\ \Delta {{E}_{2}} & =\frac{\Delta k}{4{{\pi}^{2}}}\frac{{{J}^{2}}{{b}^{2}}}{16{{v}_{\text{F}}}}\ln \left[\frac{x+\sqrt{{{(\frac{Jb}{4{{v}_{\text{F}}}})}^{2}}+{{x}^{2}}}}{{{k}_{0}}+x+\sqrt{{{(\frac{Jb}{4{{v}_{\text{F}}}})}^{2}}+{{({{k}_{0}}+x)}^{2}}}}\right], \end{array}\end{eqnarray} \tag{ 54 }$$

where $\Delta k$ and k₀ define the circular rectangle of integration around the hot spots (see figure 19(a)), and $x=\Delta {{k}^{2}}/8{{R}_{0}}\ll 1$. By adding the energy cost and the energy gain, and expanding in small J and b, we obtain that the leading order contribution is

$$\begin{eqnarray}&&\Delta {{E}_{1}}+\Delta {{E}_{2}}\simeq -\frac{\Delta k}{4{{\pi}^{2}}}\frac{{{J}^{4}}{{b}^{2}}}{16v_{F}^{2}}\frac{1}{{{k}_{0}}+x}\left(\frac{1}{x}-x\right)\ln \left(\frac{{{k}_{0}}+x}{x}\right).\end{eqnarray} \tag{ 55 }$$

It is clear that $\Delta {{E}_{1}}+\Delta {{E}_{2}}<0$ in the region where the perturbative approach is valid ($b\ll 1$). Moreover, we have found that the energy difference between the vortex crystal state and the single-$\mathbf{Q}$ ordering is of order ${{J}^{4}}{{b}^{2}}$ for small J and b. This result agrees with our naive expectations based on equations (48) and (49). However, we have arrived to this result without encountering any divergence on the way.

Equation (54) reveals the origin of the non-coplanar magnetic ordering: the single-${{\mathbf{Q}}_{1}}$ state still has a high susceptibility for developing an additional ${{\mathbf{Q}}_{2}}$ modulation. In particular, this susceptibility is logarithmically divergent in the limit of perfect nesting (${{R}_{0}}\to \infty$), i.e. the energy gain $\Delta {{E}_{2}}$ becomes proportional to ${{b}^{2}}\ln b$, while the energy cost $\Delta {{E}_{1}}$ remains proportional to b².

Figure 19(c) shows the E(b) curve obtained from our perturbative approach after integrating over the $1\text{BZ}$ (not just around the hot spots). For small b values, E(b) is in good agreement with the exact evaluation for the double-Q state in equation (45), which can be computed by direct diagonalization of $\mathcal{H}$ when the magnetic unit cell associated with the wave-vectors ${{\mathbf{Q}}_{\nu}}$ is not too large. A deviation $\mathcal{O}\left({{10}^{-8}}\right)=\mathcal{O}$(${{J}^{2}}{{b}^{6}}$, ${{J}^{4}}{{b}^{4}}$, ${{J}^{6}}{{b}^{2}}$) for b  =  0.1 is consistent with the expected accuracy of the perturbative expansion. The stabilization of the non-coplanar vortex crystal described by equation (45) has been confirmed with unbiased Langevin dynamics simulations [88].

The above analysis shows that contributions which are not included at the RKKY level can radically change the nature of the magnetic ordering of frustrated itinerant magnets. In particular, these terms favor multiple-$\mathbf{Q}$ non-coplanar magnetic orderings which, in turn, produce an effective spin–orbit coupling via the Berry phase mechanism that was discussed at the beginning of this section. In particular, the chirality stripes of the non-coplanar vortex crystal described by equation (45) produce staggered electronic currents in the direction parallel to the stripes. The purpose of the next section is to show that the stabilization mechanism discussed here leads to other interesting non-coplanar magnetic orderings for special filling fractions.

The FS of the half-filled square lattice with nearest-neighbor hopping is a perfect square inscribed into the bigger square boundary of the first Brillouin zone (see figure 20(a)). This implies that the Fermi surface is perfectly nested with a single nesting wave-vector $\mathbf{Q}=\left(\pi,\pi \right)$. The corners of the square correspond to saddle points of the dispersion relation ${{\varepsilon}_{\mathbf{k}}}$, implying that the density of states diverges logarithmically at these points (Van-Hove singularity). The combination of perfect nesting and the Van-Hove singularity leads to magnetic susceptibility ${{\chi}_{\mathbf{q}}}$ that diverges as $-{{\ln}^{2}}\,|\mathbf{Q}-\mathbf{q}|$. Naturally, this divergence is stronger than the $-\ln q$ divergence of the uniform susceptibility, which is proportional to the density of states. Therefore, as it is shown in figure 20(a), the dominant magnetic ordering for the weak-coupling regime of the half-filled square KL model (KLM) is a plain collinear antiferromagnet with ordering wave-vector $\mathbf{Q}=\left(\pi,\pi \right)$. Frustration is absent in this example because the magnetic susceptibility has a single global maximum.

Zoom In Zoom Out Reset image size

We will consider now the case of the three-quarter filled triangular lattice with nearest-neighbor hopping [86]. Similarly to the case of the half-filled triangular lattice, the FS consists of a prefect regular hexagon that is inscribed in the bigger hexagonal boundary of the first Brillouin zone (see figure 20(b)). However, in contrast to the square lattice case, the hexagonal FS has three nesting wave-vectors ${{\mathbf{Q}}_{\nu}}$ with $\nu =1,2,3$. Like in the previous case, ${{\mathbf{Q}}_{\nu}}$ are half of reciprocal lattice vectors and the corners of the FS are saddle points of ${{\varepsilon}_{\mathbf{k}}}$. Consequently, the magnetic susceptibility diverges at each of these wave-vectors as $-{{\ln}^{2}}\,|{{\mathbf{Q}}_{\nu}}-\mathbf{q}|$, i.e. it has three degenerate global maxima. Given the frustrated nature of the magnetic susceptibility, it is natural to ask what is the optimal magnetic ordering.

To answer this question we will first consider the RKKY Hamiltonian $\mathcal{H}_{\text{eff}}^{(2)}$ given in (44). It is clear that a ground state of $\mathcal{H}_{\text{eff}}^{(2)}$ cannot have $\mathbf{q}$-components different from ${{\mathbf{Q}}_{\nu}}$. Consequently, the ground-space of $\mathcal{H}_{\text{eff}}^{(2)}$ is parametrized by the following expression:

$$\begin{eqnarray}&&{{\mathbf{S}}_{\mathbf{r}}}=\underset{\nu =1,3}{\sum}{{\mathbf{\Delta}}_{\nu}}\cos {{\mathbf{Q}}_{\nu}}\cdot \mathbf{r},\end{eqnarray} \tag{ 56 }$$

where ${{\mathbf{\Delta}}_{\nu}}$ are the vector amplitudes of the three wave-vectors ${{\mathbf{Q}}_{\nu}}$. The normalization condition of the local classical moments, $\mathbf{S}_{\mathbf{r}}^{2}=1$, implies that these vector amplitudes cannot take arbitrary values. Given that the wave vectors ${{\mathbf{Q}}_{\nu}}$ are half of reciprocal lattice vectors, the normalization condition implies that the three vector amplitudes, ${{\mathbf{\Delta}}_{\nu}}$, have to be orthogonal to each other:

$$\begin{eqnarray}&&{{\mathbf{\Delta}}_{1}}\bot {{\mathbf{\Delta}}_{2}}\bot {{\mathbf{\Delta}}_{3}}.\end{eqnarray} \tag{ 57 }$$

This condition leaves us with three possible orderings: (a) collinear ordering if two vector amplitudes are zero; (b) coplanar ordering if one vector amplitude is zero; and (c) non-coplanar ordering if the three vector amplitudes are non-zero. As it can be anticipated from the discussion in the previous section, direct evaluation of the exact energy of ${{\mathcal{H}}_{\text{KL}}}$ for these three states shows that the non-coplanar ordering with $|{{\mathbf{\Delta}}_{1}}|\,=\,|{{\mathbf{\Delta}}_{2}}|\,=\,|{{\mathbf{\Delta}}_{3}}|$ is the actual ground state [86].

Figure 21(b) shows that the real space version of the 3-$\mathbf{Q}$ non-coplanar magnetic ordering is a four-sublattice structure in which the spins of each sublattice point towards each of the four corners of a regular tetrahedron. In particular, the solid angle subtended by three spins from the same triangle jkl is one quarter of the full solid angle of the sphere, i.e. ${{\Omega}_{jkl}}=\pi$, as it is shown in figure 21(c). Correspondingly, the Berry phase that electrons pick up when they move around each triangle is ${{\Omega}_{jkl}}/2=\pi /2$. In other words, each triangle acquires a uniform spontaneous net 'magnetic flux' ${{\phi}_{jkl}}$ equal to one quarter of the flux quantum ${{\phi}_{0}}$. The real external magnetic field that would be required to produce such a flux is 10⁴T if we assume that the lattice parameter is of the order of a few Angstroms. Interestingly enough, Akagi and Motome found the same magnetic ordering for the intermediate-coupling regime of the quarter-filled (n  =  1/4) KLM [89, 90].

Zoom In Zoom Out Reset image size

Mermin–Wagner's theorem [78] precludes long-range magnetic ordering at finite T because it breaks the continuous SU(2) symmetry of ${{\mathcal{H}}_{\text{KL}}}$. In contrast, the scalar spin chirality $\langle {{\chi}_{jkl}}\rangle$ breaks only discrete symmetries, such as time-reversal and mirror symmetries. Therefore, in contrast to the magnetic ordering, long range chiral ordering can persist up to a finite critical temperature. MC and Langevin Dynamics simulations of ${{\mathcal{H}}_{\text{KL}}}$ for filling fractions n  =  1/4 (intermediate coupling) and n  =  3/4 (weak-coupling) show that this is indeed the case, and that the chiral ordering disappears through a first order thermodynamic phase transition [87, 91]. Like in the case of Mott insulators, the local physical observable associated with the scalar spin chirality is the orbital contribution to the magnetization, i.e. the component associated with orbital motion of the conduction electrons under the presence of the local effective flux ${{\phi}_{jkl}}$. In particular, given that the triple-$\mathbf{Q}$ state that we are considering has a uniform scalar spin chirality, the global order parameter is the uniform orbital magnetization, i.e. the chiral spin liquid is just an orbital ferromagnet.

The transverse or Hall conductivity, ${{\sigma}_{xy}}$, is a non-local physical property that can also be used to characterize the chiral liquid state of an itinerant magnet. The uniform effective flux should naturally induce a finite ${{\sigma}_{xy}}$. This phenomenon is known as anomalous Hall effect (AHE) and it has been studied for many years [92]. However, most of the original proposals for finding AHE in real materials relied on the relativistic spin–orbit interaction, as a key mechanism for stabilizing the non-coplanar magnetic orderings required to produce the emergent gauge field ${{\phi}_{jkl}}\propto \langle {{\chi}_{jkl}}\rangle \ne 0$. As we are showing in this section, magnetic frustration is an alternative route to induce AHE even in absence of relativistic spin–orbit coupling.

For the case under consideration, the uniform effective flux is induced by a non-coplanar ordering that fully gaps out the original Fermi surface. This implies that the Hall conductance must be quantized [93]. The value of ${{\sigma}_{xy}}$ can be anticipated by using the following trick. We will assume that the local moments remain in the 'tetrahedral' $3-\mathbf{Q}$ ordering depicted in figure 21, while we adiabatically increase the coupling constant J all the way to the strong coupling regime $JS/t\gg 1$. It is clear that the gap only increases along this process, leading to two disconnected sectors in the strong-coupling limit: the low (high) energy subspace of states in which conduction electrons hop between neighboring sites remaining parallel (anti-parallel) to the local moment ${{\mathbf{S}}_{i}}$.

The overall energy splitting between these two sectors is 2JS. Within each sector, the conduction electron spin orientation is now enslaved to the orientation of the local moment. This implies that the effective Hamiltonian that results from projecting ${{\mathcal{H}}_{\text{KL}}}$ onto the low (high)-energy subspaces consists of spinless fermions with an effective hopping

$$\begin{eqnarray}&&{{\tilde{t}}_{jk}}={{t}_{jk}}\cos \left({{\theta}_{jk}}/2\right){{\text{e}}^{\text{i}\sigma {{\Omega}_{jk\mathbf{n}}}/2}}\end{eqnarray} \tag{ 58 }$$

where ${{\Omega}_{jk\mathbf{n}}}$ is the solid angle subtended by the local moments on the sites j and l, and a unit vector ${{\mathbf{n}}_{jkl}}$ associated to each triangle jkl. Here we only need to consider half of the triangles (let's say the triangles that point upwards) because each bond belongs to two triangles. The variable σ denotes the relative orientation between the spin of the itinerant electron and the local moment, i.e. $\sigma =1$ (−1) for the low (high) energy subspaces. The arbitrary choice of the unit vectors ${{\mathbf{n}}_{jkl}}$ corresponds to the gauge freedom of this effective low-energy theory known as double-exchange model [94]. We note that the circulation of the the phase, ${{\Omega}_{jk}}+{{\Omega}_{kl}}+{{\Omega}_{lj}}={{\Omega}_{jkl}}$, is equal to twice the gauge invariant flux of the triangle jkl, which turns out to be equal to the solid angle subtended by the local moments ${{\mathbf{S}}_{j}}$, ${{\mathbf{S}}_{k}}$ and ${{\mathbf{S}}_{l}}$, as we already anticipated in the previous section (see figure 16). As we anticipated in figure 16, this phase indistinguishable from the Aharonov–Bohm phase that an electron would pick up in an externally applied magnetic field [95].

For the specific case of 'tetrahedral' ordering that we are considering (see figure 21), we have

$$\begin{eqnarray}&&{{\tilde{t}}_{jk}}=\frac{t}{\sqrt{3}}{{\text{e}}^{\text{i}\sigma \pi /6}},\end{eqnarray} \tag{ 59 }$$

when we circulate clockwise (from j to k) around a triangle jkl that is pointing upwards. Here we are using a particular gauge choice in which the unit vectors ${{\mathbf{n}}_{jkl}}$ are parallel to the sum of the three moments j, k and l.

For the specific case of the 'tetrahedral' ordering shown in figure 21, the solid angle subtended by local moments in corners of each elementary triangular plaquette is clearly $\Omega =4\pi /4$. Thus, in the limit of strong coupling, the problem is equivalent to one of spinless electrons on a triangular lattice in a uniform spin-dependent magnetic field that produces a flux ${{\phi}_{\sigma}}=\pi /2$ per triangular plaquette. The energy spectrum in the strong-coupling limit is

$$\begin{eqnarray}&&\epsilon _{\mathbf{k}}^{\pm}=\pm JS\pm \frac{2t}{\sqrt{3}}{{\left({{\cos}^{2}}\mathbf{k}{{\mathbf{a}}_{1}}+{{\cos}^{2}}\mathbf{k}{{\mathbf{a}}_{2}}+{{\cos}^{2}}\mathbf{k}{{\mathbf{a}}_{3}}\right)}^{1/2}},\end{eqnarray} \tag{ 60 }$$

where ${{\mathbf{a}}_{1}}$ and ${{\mathbf{a}}_{2}}$ are real space lattice vectors of the triangular lattice and ${{\mathbf{a}}_{3}}={{\mathbf{a}}_{1}}-{{\mathbf{a}}_{2}}$. Since the real space (magnetic) unit cell has doubled, this dispersion is defined on half of the original Brillouin zone. There are total of four energy-split bands. Thus, the system is an insulator for fillings that are integer multiples of 1/4. Filling n  =  1/2 corresponds to full filling of the spin-up sector, and therefore it is equivalent to a band insulator with ${{\sigma}_{xy}}=0$. From the explicit treatment of spinless electrons on triangular lattice [96], one finds that $|{{\sigma}_{xy}}|\,={{e}^{2}}/h$ for the filling fraction n  =  3/4 that we are considering here. This result is not surprising considering that four triangles contain one flux quantum for the tetrahedral ordering shown in figure 21, and that n  =  3/4 corresponds to one hole per two sites, i.e. one hole per flux quantum. As we mentioned before, the same magnetic ordering is obtained at n  =  1/4 for intermediate-coupling. In this case we just need to occupy the lowest energy band of spin-up sector. Clearly we still have $|{{\sigma}_{xy}}|\,={{e}^{2}}/h$ because in this case we have one electron per-flux quantum. However, ${{\sigma}_{xy}}$ has an opposite sign in this case because spin up and down electrons experience opposite fluxes.

Because the gap does not close when we reduce J all the way from the strong to the weak-coupling limit, relevant for n  =  3/4, $|{{\sigma}_{xy}}|$ remains equal to e²/h. This can be verified by explicit calculation of the Kubo formula directly in the weak-coupling regime [86]. In other words, the insulator that is induced by the tetrahedral ordering depicted in figure 21 is a Chern insulator. This phenomenon is known as quantum anomalous Hall effect (QAHE).

Our simple example of 3-$\mathbf{Q}$ in the triangular lattice clearly shows that spontaneous non-coplanar magnetic orderings can be exploited to modify the topological properties of the reconstructed bands and for instance induce QAH effect at ambient temperature [87]. As we explained in the previous section, magnetic frustration is reflected in the multiple wave-vectors that maximize the electronic susceptibility. By now there are several examples of non-coplanar orderings that produce QAH effect for different lattices and filling fractions [85, 89, 97–107]. It has also been shown that a QAH state can naturally arise form KL Hamiltonians of itinerant electrons interacting with chiral ice states of local moments, such as Kagome ice [103, 108]. The robustness of the QAH state against thermal fluctuations that do not destroy the chiral ordering of the local moments is a direct consequence of the robustness of the integer quantum Hall state against disorder. It is interesting to note that the QAH state can be suppressed by application of an external magnetic field. The Zeeman coupling between the field and the local moments favors a fully polarized state, implying that the charge gap must close at critical value of the field [87].

The non-coplanar four-sublattice ordering that we discussed in this section (see figure 21) has been proposed to describe the magnetic ordering in Mn monolayers on Cu(1 1 1) surfaces [109]. In general, non-coplanar magnetic orderings have been recently observed in 2D layers of 3d metals deposited on a nonmagnetic metallic surface [81, 110]. Moreover, first-principles calculations for these systems [109–110] indicate that these non-coplanar orderings arise from rather strong effective four-spin interactions. The results that we discussed in this section unveil the generic origin of these effective interactions and explain why they are so ubiquitous in frustrated itinerant magnets.

Equation (9) shows the functional form of the effective charge operator ${{\tilde{n}}_{j}}={{P}_{\phi}}{{\text{e}}^{\mathcal{S}}}{{n}_{j}}{{\text{e}}^{-\mathcal{S}}}{{P}_{\phi}}$ to $O\left({{t}^{3}}/{{U}^{3}}\right)$ for the strong-coupling limit of the hall-filled Hubbard model. As we discussed in the main text, this contribution from the smallest loop of frustrated lattices that is a triangle (jkl) (see figure A1(a)). Below we derive the proportionality constant ${{\beta}_{jkl}}=8{{t}_{jk}}{{t}_{kl}}{{t}_{lj}}/{{U}^{3}}+\mathcal{O}\left({{t}^{5}}/{{u}^{5}}\right)$ in this expression. To this end, it is sufficient to evaluate a particular nonzero matrix element on both sides of ${{\tilde{n}}_{j}}={{P}_{\phi}}{{\text{e}}^{\mathcal{S}}}{{n}_{j}}{{\text{e}}^{-\mathcal{S}}}{{P}_{\phi}}$ by choosing a convenient basis set. Here we consider the subspace characterized by (1) the total spin S  =  1/2, (2) S^z  =  1/2, and (3) even parity under exchanging k and l, i.e. ${{c}_{k,\sigma}}\to {{c}_{l,\sigma}}$, and ${{c}_{l,\sigma}}\to {{c}_{k,\sigma}}$. The corresponding basis states are (see figure A1(b))

$$\begin{eqnarray}\begin{array}{*{35}{l}} |{{\phi}_{1}}\rangle & =\frac{1}{\sqrt{2}}c_{j\uparrow}^{\dagger}(c_{k\uparrow}^{\dagger}c_{l\downarrow}^{\dagger}-c_{k\downarrow}^{\dagger}c_{l\uparrow}^{\dagger})|0\rangle, \\ |{{\phi}_{2}}\rangle & =\frac{1}{\sqrt{2}}c_{j\uparrow}^{\dagger}(c_{k\uparrow}^{\dagger}c_{k\downarrow}^{\dagger}+c_{l\uparrow}^{\dagger}c_{l\downarrow}^{\dagger})|0\rangle, \\ |{{\phi}_{3}}\rangle & =\frac{1}{\sqrt{2}}c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}(c_{k\uparrow}^{\dagger}+c_{l\uparrow}^{\dagger})|0\rangle, \\ |{{\phi}_{4}}\rangle & =\frac{1}{\sqrt{2}}(c_{k\uparrow}^{\dagger}c_{k\downarrow}^{\dagger}c_{l\uparrow}^{\dagger}+c_{k\uparrow}^{\dagger}c_{l\uparrow}^{\dagger}c_{l\downarrow}^{\dagger})|0\rangle, \end{array}\end{eqnarray} \tag{ A.1 }$$

where only $|{{\phi}_{1}}\rangle$ is a pure spin state with

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{\phi}_{1}}|{{\mathbf{S}}_{j}}\cdot ({{\mathbf{S}}_{k}}+{{\mathbf{S}}_{l}})-2{{\mathbf{S}}_{k}}\cdot {{\mathbf{S}}_{l}}|{{\phi}_{1}}\rangle =\frac{3}{2}. \end{array}\end{eqnarray} \tag{ A.2 }$$

Zoom In Zoom Out Reset image size

For the moment we assume ${{t}_{jk}}={{t}_{kl}}={{t}_{lj}}=t$ and denote the kinetic term in ${{\mathcal{H}}_{\text{Hub}}}$ as $\hat{T}$. The matrix elements of $\hat{T}$ and $\delta {{n}_{j}}={{n}_{j}}-1$ are represented, respectively, as

$$\begin{eqnarray}\langle {{\phi}_{\mu}}|\hat{T}|{{\phi}_{\nu}}\rangle =\left(\begin{array}{c} 0 & -2t & t & t \\ -2t & 0 & t & -t \\ t & t & -t & 0 \\ t & -t & 0 & t \end{array}\right),\\ \langle {{\phi}_{\mu}}|\delta {{n}_{j}}|{{\phi}_{\nu}}\rangle =\left(\begin{array}{c} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 \end{array}\right).\end{eqnarray} \tag{ A.3 }$$

We denote the unperturbed Hamiltonian as ${{H}_{0}}=(U/2)$ $\sum\nolimits_{j}{{\left({{n}_{j}}-1\right)}^{2}}$ where a chemical potential term to ensure the half-filling has been included. By denoting ${{G}_{0}}(E)={{\left(E-{{H}_{0}}\right)}^{-1}}$ and $Q=1-{{P}_{\phi}}$, we compute $|{{\psi}_{1}}\rangle ={{\text{e}}^{-\mathcal{S}}}|{{\phi}_{1}}\rangle$ to $\mathcal{O}\left({{t}^{3}}/{{U}^{3}}\right)$ as follows by using perturbation theory:

$$\begin{eqnarray}\begin{array}{*{35}{l}} |{{\psi}_{1}}\rangle = & |{{\phi}_{1}}\rangle +{{G}_{0}}(E)Q\hat{T}|{{\phi}_{1}}\rangle +{{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}|{{\phi}_{1}}\rangle \\ {} & +\,{{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}|{{\phi}_{1}}\rangle +\mathcal{O}\left({{t}^{4}}/{{U}^{4}}\right), \end{array}\end{eqnarray} \tag{ A.4 }$$

where

$$\begin{eqnarray}\begin{array}{r l} {{G}_{0}}(E)Q\hat{T}|{{\phi}_{1}}\rangle = & \frac{t}{E-U}\left(-2|{{\phi}_{2}}\rangle +|{{\phi}_{3}}\rangle +|{{\phi}_{4}}\rangle \right), \\ {{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}|{{\phi}_{1}}\rangle = & {{\left(\frac{t}{E-U}\right)}^{2}}\left(-3|{{\phi}_{3}}\rangle +3|{{\phi}_{4}}\rangle \right), \\ {{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}|{{\phi}_{1}}\rangle = & {{\left(\frac{t}{E-U}\right)}^{3}}\left(-6|{{\phi}_{2}}\rangle +3|{{\phi}_{3}}\rangle +3|{{\phi}_{4}}\rangle \right). \end{array}\end{eqnarray} \tag{ A.5 }$$

In fact we will not use the last piece because $\delta {{n}_{j}}$ is diagonal in the present basis set. We obtain

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle {{\psi}_{1}}|\delta {{n}_{j}}|{{\psi}_{1}}\rangle & =12{{\left(\frac{t}{U}\right)}^{3}}+\mathcal{O}\left({{t}^{5}}/{{U}^{5}}\right), \end{array}\end{eqnarray} \tag{ A.6 }$$

which we should compare with $\langle {{\phi}_{1}}|\left({{\tilde{n}}_{j}}-1\right)|{{\phi}_{1}}\rangle =3\beta /2$ (see equations (9) and (A.2)). Thus we conclude

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\beta}_{jkl}} & =8{{\left(\frac{t}{U}\right)}^{3}}+\mathcal{O}\left(\frac{{{t}^{5}}}{{{U}^{5}}}\right)\to 8\frac{{{t}_{jk}}{{t}_{kl}}{{t}_{lj}}}{{{U}^{3}}}+\mathcal{O}\left(\frac{{{t}^{5}}}{{{U}^{5}}}\right), \end{array}\end{eqnarray} \tag{ A.7 }$$

where in the last line we have recovered the bond dependence of the hopping integrals, the correctness of which can be suggested by the fact that we have computed the lowest order $\mathcal{O}\left({{t}^{3}}/{{U}^{3}}\right)$ term arising from a triangular loop.

Next we discuss how to determine the prefactor ${{\gamma}_{jkl}}$ that appears in the effective current operator (equation (8)). For this purpose, the above basis set $|{{\phi}_{1}}\rangle$ – $|{{\phi}_{4}}\rangle$ is not appropriate because $\langle {{\phi}_{1}}|{{\mathbf{S}}_{j}}\cdot \left({{\mathbf{S}}_{k}}\times {{\mathbf{S}}_{l}}\right)|{{\phi}_{1}}\rangle =0$. Instead, we use a different basis generating a chiral subspace characterized by (1) S  =  1/2, (2) S^z  =  1/2, and (3) the eigenvalue $\omega ={{\text{e}}^{2\pi \text{i}/3}}$ under the C₃ rotation (see figure A1(c)):

$$\begin{eqnarray}\begin{array}{*{35}{l}} |\phi _{1}^{\prime}\rangle & =\frac{1}{\sqrt{3}}(c_{j\uparrow}^{\dagger}c_{k\uparrow}^{\dagger}c_{l\downarrow}^{\dagger}+\omega c_{j\downarrow}^{\dagger}c_{k\uparrow}^{\dagger}c_{l\uparrow}^{\dagger}+{{\omega}^{2}}c_{j\uparrow}^{\dagger}c_{k\downarrow}^{\dagger}c_{l\uparrow}^{\dagger})|0\rangle, \\ |\phi _{2}^{\prime}\rangle & =\frac{1}{\sqrt{3}}(c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}c_{k\uparrow}^{\dagger}+\omega c_{k\uparrow}^{\dagger}c_{k\downarrow}^{\dagger}c_{l\uparrow}^{\dagger}+{{\omega}^{2}}c_{j\uparrow}^{\dagger}c_{l\uparrow}^{\dagger}c_{l\downarrow}^{\dagger})|0\rangle, \\ |\phi _{3}^{\prime}\rangle & =\frac{1}{\sqrt{3}}(c_{j\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}c_{l\uparrow}^{\dagger}+\omega c_{j\uparrow}^{\dagger}c_{k\uparrow}^{\dagger}c_{k\downarrow}^{\dagger}+{{\omega}^{2}}c_{k\uparrow}^{\dagger}c_{l\uparrow}^{\dagger}c_{l\downarrow}^{\dagger})|0\rangle, \end{array}\end{eqnarray} \tag{ A.8 }$$

where only $|\phi _{1}^{\prime}\rangle$ is a pure spin state and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \phi _{1}^{\prime}|{{\mathbf{S}}_{j}}\cdot ({{\mathbf{S}}_{k}}\times {{\mathbf{S}}_{l}})|\phi _{1}^{\prime}\rangle =\frac{\sqrt{3}}{4}. \end{array}\end{eqnarray} \tag{ A.9 }$$

By using the matrix elements of $\hat{T}$ and $\hslash {{I}_{jk}}=$ $-\text{i}{{t}_{jk}}\sum\nolimits_{\sigma}(c_{k\sigma}^{\dagger}{{c}_{j\sigma}}-c_{j\sigma}^{\dagger}{{c}_{k\sigma}})$, namely,

$$\begin{eqnarray}\begin{array}{*{20}{l}}\langle \phi _{\mu}^{\prime}|\hat{T}|\phi _{\nu}^{\prime}\rangle &=\left(\begin{array}{c} 0 & t\left(1-{{\omega}^{2}}\right) & t\left({{\omega}^{2}}-\omega \right) \\ t\left(1-\omega \right) & 0 & t\left(\omega -1\right) \\ t\left(\omega -{{\omega}^{2}}\right) & t\left({{\omega}^{2}}-1\right) & 0 \end{array}\right),\\ \hslash \langle \phi _{\mu}^{\prime}|{{I}_{jk}}|\phi _{\nu}^{\prime}\rangle & =\left(\begin{array}{c} 0 & \text{i}t\left({{\omega}^{2}}-1\right)/3 & \text{i}t\left({{\omega}^{2}}-\omega \right)/3 \\ -\text{i}t\left(\omega -1\right)/3 & 0 & \text{i}t\left(1-\omega \right)/3 \\ -\text{i}t\left(\omega -{{\omega}^{2}}\right)/3 & -\text{i}t\left(1-{{\omega}^{2}}\right)/3 & 0 \end{array}\right),\end{array}\end{eqnarray} \tag{ A.10 }$$

respectively, we evaluate $|\psi _{1}^{\prime}\rangle ={{\text{e}}^{-S}}|\phi _{1}^{\prime}\rangle$ perturbatively in the same way as in equation (A.4). We obtain

$$\begin{eqnarray}\begin{array}{r l} {{G}_{0}}(E)Q\hat{T}|\phi _{1}^{\prime}\rangle & =\frac{t\left(1-\omega \right)}{E-U}\left(|\phi _{2}^{\prime}\rangle +\omega |\phi _{3}^{\prime}\rangle \right), \\ {{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}|\phi _{1}^{\prime}\rangle & ={{\left(\frac{t\left(1-\omega \right)}{E-U}\right)}^{2}}[-\omega |\phi _{2}^{\prime}\rangle\\ &\quad -\left(1+\omega \right)|\phi _{3}^{\prime}\rangle ], \\ {{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}{{G}_{0}}(E)Q\hat{T}|\phi _{1}^{\prime}\rangle & ={{\left(\frac{t\left(1-\omega \right)}{E-U}\right)}^{3}}[\left(1+\omega \right)|\phi _{2}^{\prime}\rangle\\ &\quad -|\phi _{3}^{\prime}\rangle ], \end{array}\end{eqnarray} \tag{ A.11 }$$

and thus,

$$\begin{eqnarray}\begin{array}{*{35}{l}} \hslash \langle \psi _{1}^{\prime}|{{I}_{jk}}|\psi _{1}^{\prime}\rangle & =-\frac{6\sqrt{3}{{t}^{3}}}{{{U}^{2}}}+\mathcal{O}\left({{t}^{5}}/{{U}^{4}}\right). \end{array}\end{eqnarray} \tag{ A.12 }$$

By comparing this with equation (A.9), we conclude

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\gamma}_{jkl}} & =-\frac{24{{t}^{3}}}{{{U}^{2}}}+\mathcal{O}\left({{t}^{5}}/{{U}^{4}}\right)\to -24\frac{{{t}_{jk}}{{t}_{kl}}{{t}_{lj}}}{{{U}^{2}}}+\mathcal{O}\left({{t}^{5}}/{{U}^{4}}\right) \end{array}\end{eqnarray} \tag{ A.13 }$$

where we have recovered the bond dependence on the hopping integrals.