Antiferromagnetic skyrmion crystals in the Rashba Hund’s insulator on triangular lattice

Abstract

Motivated by the importance of antiferromagnetic skyrmions as building blocks of next-generation data storage and processing devices, we report theoretical and computational analysis of a model for a spin-orbit coupled correlated Hund’s insulator magnet on a triangular lattice. We find that two distinct antiferromagnetic skyrmion crystal (AF-SkX) states can be stabilized at low temperatures in the presence of external magnetic field. The results are obtained via Monte Carlo simulations on an effective magnetic model derived from the microscopic electronic Hamiltonian consisting of Rashba spin-orbit coupling, as well as strong Hund’s coupling of electrons to classical spins at half-filling. The two AF-SkX phases are understood to originate from a classical spin liquid state that exists at low but finite temperatures. These AF-SkX states can be easily distinguished from each other in experiments as they are characterized by peaks at distinct momenta in the spin structure factor which is directly measured in neutron scattering experiments. We also discuss examples of materials where the model as well as the two AF-SkX states can be realized.

Introduction

Strongly correlated electrons residing on geometrically frustrated lattices lead to intriguing ordered as well as disordered phases^1,2,3,4,5. While such systems are extremely challenging to study, suitably motivated approximate treatments not only lead to predictions of remarkable new phases of electronic matter, but also provide new paradigms for understanding the observed electronic properties of solids. Some examples of such phases are quantum and classical spin liquids^{6,7,8,9,10,11}, non-coplanar, non-collinear magnetic states^{12,13,14,15,16} and partially ordered states^17,18. Furthermore, the competition between spin-orbit coupling (SOC) and electronic correlations has emerged as one of the most interesting area of fundamental research in recent years^{19,20,21,22,23,24,25}. In particular, the possibility to tune Rashba SOC via a suitable material growth and design in terms of thin-film multilayers or interfaces has allowed for a realization of SOC induced effects in real materials. One important consequence, with potential applications in data storage and processing technologies, is the observation of skyrmions, antiskyrmions and antiferromagnetic skyrmions in various metals and insulators^{26,27,28,29,30}. Indeed, such topological spin textures are considered as building blocks of information storage in the race-track memory devices^31,32,33,34. Antiferromagnetic skyrmions are considered superior to skyrmions as the former do not exhibit skyrmion Hall effect which affects the device performance in case of skyrmions^35,36,37.

In this report, we present the results of our investigations of a prototype model that combines three of the most interesting aspects of electronic problems, namely, geometrical frustration, strong correlations and Rashba SOC. Our main motivation is to understand the physics of antiferromagnetic skyrmion formation in a microscopic electronic Hamiltonian. Most theoretical investigations of formation of skyrmion-like quasiparticles use suitable spin Hamiltonians as a starting point^38,39,40,41. Instead, here we begin with a microscopic model with itinerant electrons coupled to localized magnetic moments via Hund’s rule coupling in the presence of Rashba SOC. Such a model can be realized in thin films or interfaces of transition metal or heavy fermion compounds^{42,43,44,45,46,47,48}. We explicitly derive a low-energy magnetic Hamiltonian for the triangular lattice for the half-filled insulating case. Given the complex and competing nature of different terms in the Hamiltonian we investigate the low-temperature phases with varying external field via unbiased Monte Carlo simulation technique. In addition to the expected 120° state and the single-Q spiral states, we identify three non-trivial magnetic phases in the model: (1) a classical spin liquid (CSL) characterized via a diffuse ring pattern centered at the K points of the first Brillouin zone (BZ), (2) a 3Q AF-SkX1 characterized by hexagonal peak pattern in spin structure factor (SSF), and (3) a qualitatively distinct 6Q AF-SkX2 phase that has never been reported before. The degeneracy-induced CSL state can be understood as the parent of both the AF-SkX phases. Our study reports, that not only the AF-SkX states can be described within a microscopic electronic model, but also two distinct AF-SkX phases exist in the triangular lattice model. These two phases can be easily detected in experiments as they lead to qualitatively different peak structure in the neutron scattering data.

Derivation of the low energy spin hamiltonian

The Rashba Hund’s model (RHM) describing Rashba electrons coupled to localized magnetic moments residing on a triangular lattice is described by the lattice Hamilonian⁴⁹,

$$\begin{aligned} H_{\text {RHM}}= & \, - t \sum _{i,\gamma ,\sigma } (c^\dagger _{i, \sigma } c^{}_{i+\gamma , \sigma } + \text {H.c.}) - \text{ i }\lambda \sum _{i,\gamma , \sigma \sigma '} c_{i\sigma }^{\dagger } [{\hat{z}}\cdot ({\varvec{\tau }} \times {\hat{\varvec{\gamma }}})]_{\sigma \sigma '} c_{j\sigma '} \nonumber \\&- J_\text {H} \sum _{i} {\mathbf{S}}_i \cdot {\mathbf{s}}_i - h_z \sum _{i} S_i^z. \end{aligned}$$

(1)

Operator $c_{i\sigma }^\dagger$ ($c_{i\sigma }$) creates (annihilates) an electron at site i with spin $\sigma \in \{\uparrow , \downarrow \}$. $\varvec{\tau }$ is a vector operator with the three Pauli matrices as components. ${\mathbf{S}}_i$ denote the localized spins which we assume to be classical vectors with $|{\mathbf{S}}_i| \equiv 1$. t, $\lambda$ and $J_{\text {H}}$ denote the strengths of hopping amplitude, Rashba coupling and Hund’s rule coupling, in that order. Assuming the lattice constant to be unity, ${\hat{\gamma }} \in \{ {\mathbf{a}}_1, {\mathbf{a}}_2, {\mathbf{a}}_3 \}$ denotes the basis unit vector of the triangular Bravais lattice with ${\mathbf{a}}_1$=(1,0), ${\mathbf{a}}_2$=(1/2,$\sqrt{3}/2$) and ${\mathbf{a}}_3$=(-1/2,$\sqrt{3}/2$). ${\mathbf{s}}_i$ is the electronic spin operator. $h_z$ is the strength of Zeeman field applied along z axis, and $\text {i} = \sqrt{-1}$.

The Hamiltonian specified in Eq. (1) above can be realized in a variety of magnetic compounds comprising of transition metal or rare-earth ions where more than one bands are partially filled. In addition, the existence of Rashba SOC requires inversion symmetry breaking which can be realized in thin films or at interfaces^{42,43,44,45,46,50,51,52,53,54,55}. We are interested in a situation where charge degree of freedom is completely frozen due to strong correlations and the low energy physics is described in terms of an effective magnetic Hamiltonian. Such a condition is met in Mott insulators where strong Hubbard term disfavours any transfer of charge. In the Hund’s model, a similar scenario is realized for large $J_{\text {H}}$ at half filling. For large $J_{\text {H}}$, it is useful to work in a site-dependent spin-quantization basis achieved via local SU(2) rotations, given by,

$$\begin{aligned} \begin{bmatrix} c_{i\uparrow } \\ c_{i\downarrow } \end{bmatrix} = \begin{bmatrix} \cos (\frac{\theta _i}{2}) &\, -\sin (\frac{\theta _i}{2}) e^{-\text {i} \phi _i} \\ \sin (\frac{\theta _i}{2}) e^{\text {i} \phi _i} &\, \cos (\frac{\theta _i}{2}) \end{bmatrix} \begin{bmatrix} d_{ip} \\ d_{ia} \end{bmatrix}. \end{aligned}$$

Here, $d_{ip} (d_{ia})$ annihilates an electron at site i with spin parallel (anti-parallel) to the localized spin. The polar and azimuthal angle pair {$\theta _i, \phi _i$} specifies the orientation in three dimensions of the local moment ${\mathbf{S}}_i$.

The transformed Hamiltonian takes the form,

$$\begin{aligned} H_{\text{RHM}} = \sum _{i, \gamma , \sigma \sigma '} \big ( g_{i, \gamma }^{\sigma \sigma '} d_{i, \sigma }^{\dagger } d_{i+\gamma , \sigma '} + \text {H.c.} \big ) - \frac{J_\text {H}}{2} \sum _i \big (n_{ip} - n_{ia} \big ), \end{aligned}$$

(2)

where $\sigma \in \{p,a\}$. The transformation to local basis puts the interaction term in a diagonal form. However, the spin dependence now resides in the hopping parameters. The projected hopping amplitudes $g_{i,\gamma }^{\sigma \sigma '} = t_{i,\gamma }^{\sigma \sigma '} + \lambda _{i,\gamma }^{\sigma \sigma '}$ have contributions from standard tight-binding hopping integral t and the Rashba SOC $\lambda$. Following the 2nd order perturbation approach applied to an isolated pair of sites, we derive the classical super-exchange (CSE) model for the triangular lattice⁵⁶.

The parallel to anti-parallel hopping contributions $g_{i,\gamma }^{pa}$, are given by,

$$\begin{aligned} t_{i,\gamma }^{pa}= & \, -t \Big [\sin (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) e^{-\text {i} \phi _i} -\cos (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{-\text {i} \phi _j} \Big ], \nonumber \\ \lambda _{i,{\mathbf{a}}_1}^{pa}= & \, \lambda \Big [\cos (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) + \sin (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{-\text {i} (\phi _i+\phi _j)} \Big ], \nonumber \\ \lambda _{i,{\mathbf{a}}_2}^{pa}= & \, \frac{\lambda }{2} \Big [(1-\text {i}\sqrt{3})\cos (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) + (1+\text {i}\sqrt{3})\sin (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{-\text {i} (\phi _i+\phi _j)} \Big ], \nonumber \\ \lambda _{i,{\mathbf{a}}_3}^{pa}= & \, -\frac{\lambda }{2} \Big [(1+\text {i}\sqrt{3})\cos (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) + (1-\text {i}\sqrt{3})\sin (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{-\text {i} (\phi _i+\phi _j)} \Big ]. \end{aligned}$$

(3)

The anti-parallel to parallel contributions, $g_{ij}^{ap}$, are given by,

$$\begin{aligned} t_{i,\gamma }^{ap}= & \, -t \Big [\cos (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{\text {i} \phi _j} - \sin (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) e^{\text {i} \phi _i} \Big ], \nonumber \\ \lambda _{i,{\mathbf{a}}_1}^{ap}= & \, -\lambda \Big [\cos (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) + \sin (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{\text {i} (\phi _i+\phi _j)} \Big ], \nonumber \\ \lambda _{i,{\mathbf{a}}_2}^{ap}= & \, -\frac{\lambda }{2} \Big [(1+\text {i}\sqrt{3}) \cos (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) + (1-\text {i}\sqrt{3}) \sin (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{\text {i} (\phi _i +\phi _j)} \Big ], \nonumber \\ \lambda _{i,{\mathbf{a}}_3}^{ap}= & \, \frac{\lambda }{2} \Big [(1-\text {i}\sqrt{3})\cos (\frac{\theta _i}{2}) \cos (\frac{\theta _j}{2}) + (1+\text {i}\sqrt{3})\sin (\frac{\theta _i}{2}) \sin (\frac{\theta _j}{2}) e^{\text {i} (\phi _i+\phi _j)} \Big ]. \end{aligned}$$

(4)

In the above equations, $j = i+\gamma$. The general expression for the second order perturbative energy correction,

$$\begin{aligned} \Delta E^{(2)}_{\mathbf {\{S\}}} = \sum\limits_k \frac{ |\left\langle \psi _k|H^{\prime }|\psi _0\right\rangle |^2 }{ E_0-E_k } = -\Bigg [ \frac{|g_{ij}^{pa}|^2}{J_\text {H}} + \frac{|g_{ij}^{ap}|^2}{J_\text {H}} \Bigg ], \end{aligned}$$

(5)

involves modulus squares of the hopping amplitudes which are obtained from Eqs. (3) and (4) in the following form:

$$\begin{aligned} |g^{pa}_{i,{\mathbf{a}}_1}|^2 = |g^{ap}_{i,{\mathbf{a}}_1}|^2= & \, \bigg [ \frac{t^2}{2}(1-{\mathbf{S}}_i \cdot {\mathbf{S}}_j) + \frac{\lambda ^2}{2} (1+{\mathbf{S}}_i \cdot {\mathbf{S}}_j-2S_i^{y}S_j^{y}) - t\lambda (S_i^zS_j^x - S_i^xS_j^z) \bigg ], \nonumber \\ |g^{pa}_{i,{\mathbf{a}}_2}|^2 = |g^{ap}_{i,{\mathbf{a}}_2}|^2= & \, \bigg [\frac{t^2}{2}(1-{\mathbf{S}}_i \cdot {\mathbf{S}}_j) + \frac{\lambda ^2}{4} \{ 2(1+S_i^zS_j^z) - S_i^xS_j^x + S_i^yS_j^y + \sqrt{3}S_i^yS_j^x + \sqrt{3}S_i^xS_j^y \} \nonumber \\&+ \frac{t\lambda }{2} \{ (S_i^xS_j^z - S_i^zS_j^x) + \sqrt{3}(S_i^yS_j^z - S_i^zS_j^y) \} \bigg ], \nonumber \\ |g^{pa}_{i,{\mathbf{a}}_3}|^2 = |g^{ap}_{i,{\mathbf{a}}_3}|^2= & \, \bigg [t^2(1-{\mathbf{S}}_i \cdot {\mathbf{S}}_j) + \frac{\lambda ^2}{2}\{ 2(1+S_i^zS_j^z) - S_i^xS_j^x + S_i^yS_j^y - \sqrt{3}S_i^yS_j^x - \sqrt{3}S_i^xS_j^y \} \nonumber \\&- \frac{t\lambda }{2} \{ (S_i^xS_j^z - S_i^zS_j^x) - \sqrt{3}(S_i^yS_j^z - S_i^zS_j^y) \} \bigg ]. \end{aligned}$$

(6)

Substituting the expression from Eq. (6) into Eq. (5) and taking the sum over all nn pairs, we arrive at the classical super-exchange (CSE) Hamiltonian on a triangular lattice,

$$\begin{aligned} H_{\text {CSE}}= & \, -1/J_{\text {H}} \sum _{i,\gamma } \big [ t^2(1-{\mathbf{S}}_i \cdot {\mathbf{S}}_j) + 2t\lambda \hat{\gamma '} \cdot ({\mathbf{S}}_i \times {\mathbf{S}}_j) \nonumber \\&+ \lambda ^2(1+{\mathbf{S}}_i \cdot {\mathbf{S}}_j-2 (\hat{\gamma '} \cdot {\mathbf{S}}_i)(\hat{\gamma '} \cdot {\mathbf{S}}_j)) \big ] - h_z \sum _{i} S_i^z, \end{aligned}$$

(7)

with $\hat{\gamma '} = {\hat{z}} \times {\hat{\gamma }}$. We also note that the Hamiltonian is written in a general form which is valid for any Bravais lattice.

Figure 1

$1/J_{\text {H}}$ dependence of energy per site calculated via exact diagonalization of the Rashba Hund’s Hamilonian (red circles) and that calculated from the low-energy spin Hamilonian (blue triangles) for different magnetic states and for different $\alpha$: (a) 120° state, (b) random configuration representative of a paramagnetic state, (c) single-Q spiral state, and (d) ferromagnetic state. All energies are measured in units of the bare hopping parameter $t_0 = 1$.

We parameterize by $\alpha$ the relative strength of the Rashba coupling as compared to hopping parameter as $t = (1-\alpha ) t_0$ and $\lambda = \alpha t_0$, where $t_0=1$ sets the reference energy scale. The resulting model consists of antiferromagnetic coupling terms along with anisotropic terms resembling Dzyaloshinskii–Moriya (DM) and Kitaev-like interactions^56,57,58 . The ground state of the Hamiltonian Eq. (7) for $\alpha =0$ is the well known three-sublattice 120° state which stabilizes due to geometrical frustration. In absence of external magnetic field ($h_z = 0$), increasing $\alpha$ favours non-collinear spin arrangement due to the presence of DM terms. However, the frustrating nature of the DM and Kiteav-like terms leads to a large degeneracy of states, as discussed by us for the case of square lattice⁵⁶. One consequence of this large degeneracy is the presence of entropically stabilized filamentary domain states at low temperatures.

Figure 2

Energy difference, $\Delta E$, between different pairs of magnetic states, calculated within exact RHM and the derived classical super-exchange model (CSE), for various $\alpha$ values. The magnetic phases considered to calculate the energy differences are antiferromagnetic (AFM), ferromagnetic (FM), single-Q (SQ), and antiferromagnetic skyrmion states (AF-SkX1, AF-SkX2) states. All energies are measured in units of the bare hopping parameter $t_0 = 1$.

Before proceeding with the investigations of the magnetic properties of the electronic Hamiltonian Eq. (1) in terms of the effective low energy CSE Hamiltonian Eq. (7), we explicitly check the validity of CSE Hamiltonian by comparing energies of different magnetic states obtained within the exact and approximate Hamiltonians. The energies calculated via the CSE Hamiltonian match very well with the exact values, provided $1/J_{\text {H}} < 0.1$ (see Fig. 1). Note that the large $J_{\text {H}}$ expansion is only valid when the parallel and antiparallel bands are split and the chemical potential resides in the gap. For the triangular lattice tight-binding bands, such splitting will occur for $J_{\text {H}} = 9$ for the largest bandwidth corresponding to ferromagnetic background. Therefore, the comparison in Fig. 1 shows that the effective Hamilonian is quantitatively accurate as long as the gap-opening condition is satisfied. In order to further demonstrate the validity of the derived model, we compare the energy differences between various pairs of magnetic states calculated within exact and approximate models (see Fig. 2). Given that the CSE model is explicitly derived via perturbation theory, it is not surprising that this serves as an accurate model for describing magnetism of the $H_{\text {RHM}}$ in the large $J_{\text {H}}$ limit.

Results and discussion

Phases in the absence of external magnetic field

In this section, we present the Monte Carlo simulation results for the CSE Hamiltonian Eq. (7) (see “Methods” section). We begin by studying the model in the absence of external magnetic field. In the limit $\alpha \rightarrow 0$, the Hamilonian reduces to a simple antiferromagnetic Heisenberg model and the lowest energy is obtained for a three-sublattice 120° arrangement of spins. We track the temperature dependence of the spin structure factor (SSF), as obtained in simulations, at relevant value of ${\mathbf{q}}$. The 120° state is characterized by a SSF peak at ${\mathbf{q}}_0 = (2\pi /3, 2\pi /\sqrt{3})$, and the symmetry related points (see Fig. 3f). For small values of $\alpha$, the SSF at ${\mathbf{q}}_0$ exhibits an order parameter like rise upon lowering temperature (see Fig. 3a). The sharp upturn point is identified as the ordering temperature, which is further verified via specific heat calculations (filled symbols in Fig. 3a). Beyond a critical value of $\alpha$, the 120° ground state is destabilized in favour of a single-Q (SQ) spiral state with SSF peak located at $\pm {\mathbf{q}}$ for one specific ${\mathbf{q}}$ (see Fig. 3h). In contrast to the 120° state, the SQ state lifts the three-fold degeneracy as one pair of K points is spontaneously selected from three equivalent choices. The ordering temperature for the SQ state, as inferred from the temperature dependence of SSF at relevant ${\mathbf{q}}$, monotonically decreases upon increasing $\alpha$ (see Fig. 3b). For $\alpha > 0.43$, a specific SQ state which consists of ferromagnetic (FM) chains oriented antiferromagnetically, labelled as spin stripe (SS) state, is stable over a wide range of $\alpha$. This is characterized by a peak in the SSF at a pair of M points (see Fig. 3i). Similar to the SQ state, the SS state is also three-fold degenerate and the degeneracy is spontaneously lifted. Furthermore, the specific heat displays a sharp peak at the temperature corresponding to the on-set of SSF peak (see Fig. 3c) allowing us to reliably infer the ordering temperatures. Eventually a FM state, characterized by magnetization, becomes the ground state in the limit of strong Rashba coupling. As discussed above, the transitions from the high temperature paramagnetic (PM) state to any of the ordered states can be described with the help of the relevant peak in the SSF. These transitions are also identified as sharp peaks in the specific heat as shown in Fig. 3a,c. We find that for most of the $\alpha$ values, the sharp rise in the order parameter is accompanied by a peak in the specific heat. However, for $0.07< \alpha < 0.25$ we find a broad hump feature in the specific heat in addition to a sharp peak (see Fig. 3d). On a careful observation of the SSF, we find the presence of diffuse circular patterns centered about the K-points of the first BZ (see Fig. 3g). This allows us to identify this finite-temperature phase as a classical spin liquid, similar to the one reported in the square lattice⁵⁶. We summarize our results of Monte Carlo simulations in the absence of magnetic field as a $T-\alpha$ phase diagram (see Fig. 3e). While there is a similarity with the square lattice phase diagram, it is surprizing to note that the geometrical frustration inherent in the triangular geometry disfavors the zero-field skyrmion state found in the square lattice⁵⁶. The zero-field skyrmion state reported in the square-lattice model consists of skyrmions packed in a square geometry which is not compatible in a triangular lattice. Therefore, the SS state is preferred over the zero-field skyrmion crystal state for $0.43< \alpha < 0.67$.

Figure 3

(a)–(c) Temperature dependence of spin structure factor (open symbols) at different values of ${\mathbf{q}}$, and for different values of relative Rashba coupling strength $\alpha$. The right y-axis in (a) and (c) is for specific heat (filled symbols) corresponding to one of the $\alpha$ values. (d) Specific heat as a function of temperature showing a broad hump followed by a sharp peak for two values of $\alpha$. (e) Phase diagram obtained by tracking the features in the SSF and the specific heat with the different phases described as follows: the 120° denotes the well known three-sublattice order on triangular lattice. Single-Q (SQ) denotes a state where SSF displays peak only at a pair ${\pm {\mathbf{q}}}$ of momentum. Spin stripe (SS) phase consists of ferromagnetic lines oriented antiferromagnetically w.r.t. neighboring lines. Classical spin liquid (CSL) denotes a phase with short-range correlations but no long-range ordering. (f)–(i) SSF plots for the four non-trivial phases displayed in the phase diagram. Note that the CSL state is characterized by a diffuse ring-like pattern centered at the K point. A zoomed in view is shown in panel (g). Inset in (e) shows variation in the magnitude q of the relevant wave-vector ${\mathbf {q}}$ with $\alpha$ for SQ state.

Phases in the presence of external magnetic field

Having established the zero-field phase diagram for the triangular lattice, we now discuss the effect of Zeeman field on the magnetic states. In Fig. 4a–d, we show temperature dependence of SSF, $C_V$ and topological susceptibility (see “Methods” section) at finite $h_z$. We discuss the case of $\alpha = 0.2$ as a representative of finite Rashba SOC. For small $h_z$, the ground state remains a three-fold degenerate SQ spiral state discussed in the previous subsection (Fig. 4e displays another choice of the degenerate ${\mathbf{q}}$ point). The ordering temperature, as inferred from the SSF at the relevant ${\mathbf{q}}$ (see Fig. 4a), decreases with increasing $h_z$. Interestingly, over a moderate range of $h_z$ values ($0.07< h_z < 0.14$) the SSF does not display any prominant peak. Specific heat also shows only a broad hump and no sharp anomaly (see Fig. 4b). We confirm the existence of diffuse circular pattern in SSF at low temperatures in this range of $h_z$, similar to that shown in Fig. 3g (see Fig. 4f). Therefore, we conclude that the CSL state gets stabilized down to very low temperatures for finite Zeeman field. This stabilization of a short-range ordered spin-liquid state by application of external magnetic field is unusual and is analogous to melting of a solid under pressure as magnetic field in spin systems is analogous to external pressure in real solids. Upon increasing $h_z$ further, we find two exotic ordered states: a 3Q antiferromagnetic skyrmion crystal, henceforth labelled as AF-SkX1, (see Fig. 4g) and a novel antiferromagnetic skyrmion crystal, henceforth labelled as AF-SkX2, with SSF peaks on the boundaries (straight line joining nearest pairs of K and M points) of the first BZ (see Fig. 4h)^41,59 .

Figure 4

(a) SSF peak as a function of temperature for SQ state, (b) specific heat variation with temperature to identify CSL state transition, (c,d) variation of SSF peak (left axis) and topological susceptibility (right axis) for AF-SkX1 and AF-SkX2 states respectively. (e–h) SSF peak locations for four magnetic states displayed in the phase diagram (Fig. 7) SQ, CSL, AF-SkX1 and AF-SkX2 respectively.

We further characterize the two multi-Q states with the help of skyrmion density ($\langle {\mathscr {T}} \rangle$) and topological susceptibility ($\chi _{{\mathscr {T}}}$) (see “Methods” section). Indeed, the topological susceptibility peaks at the on-set temperature of the multi-Q order as inferred from the SSF (see Fig. 4c,d) in both the skyrmion states, AF-SkX1 and AF-SkX2. The peaks in $\chi _{{\mathscr {T}}}$ are clear indications of non-topological to topological phase transitions.

In Fig. 5 we show the evolution with increasing Zeeman field of low temperature magnetic states via representative spin configurations. The SQ state consists of spins spiraling in the xz plane with the ordering wavevector residing on the x axis in the reciprocal space (the corresponding SSF is shown in Fig. 4e). Upon increasing magnetic field the system enters a short-range ordered CSL phase. A typical spin configuration in the CSL state consists of filamentary domain segments (see Fig. 5b). The existence of such a disordered state relies on the presence of an unusual degeneracy of the SQ spirals that involved a simultaneous change of the spiral wavevector and the spin plane. This is discussed by us in recent papers^56,57. Indeed, the CSL state is similar to the antiferromagnetic string state discussed in⁵⁶, with the difference that the CSL state emerges in the background of 120° state on triangular lattice. It is interesting to note that some of the filaments existing in the CSL state are short, and hence acquire a skyrmion-like modulations of the spins. This is suggestive that the CSL state is unstable towards a state hosting skyrmions.

Figure 5

Low-temperature spin configurations for (a) single-Q spiral, (b) classical spin liquid, (c) AF-SkX1 and (d) AF-SkX2 states. The z-component is represented by the color bar and the planar components by the arrow lengths and directions. For clarity, we only display a $30 \times 30$ section of the simulated lattice in panels (a), (c), (d) and $60 \times 60$ for (b). The insets in (c) and (d) show the zoomed-in view of spins in the core of skyrmions. The inset of (a) denotes the three-sublattices (A, B, C) in a single triangular plaquette.

Figure 6

Sublattice resolved spin configuration plots for AF-SkX1 (upper panels) and AF-SkX2 (lower panels) states. All three-sublattices, denoted as (A), (B) and (C) on the vertices of a triangular plaquette, display perfect Néel skyrmion textures.

Indeed, this is confirmed as increasing magnetic field leads to the formation of AF-SkX1 state. A typical configuration of spins in the AF-SkX1 state is shown in Fig. 5c where a triangular arrangement of antiferromagnetic skyrmions is observed in the background of 120° state. With a further increase in the strength of Zeeman field, we find the AF-SkX2 as the ground state (see Fig. 5d). While it is difficult to distinguish between AF-SkX1 and AF-SkX2 looking at the real-space spin configurations, the SSF for AF-SkX2 is qualitatively different from that for AF-SkX1 (compare Fig. 4g,h). This can be interpreted as a superposition of two counter-rotated triangular arrangements of the skyrmions. In order to understand the underlying spin structure Fig. 5c,d, we plot the sublattice resolved spin configurations^41,60. The spin textures for AF-SkX1 and AF-SkX2 on three-sublattices (A, B, C) are shown in Fig. 6. The Néel type nature of skyrmion states on each sublattice is clear from Fig. 6. These plots also clarify that the difference between AF-SkX1 and AF-SkX2 is purely in terms of how the three-sublattices are oriented relative to one another. This is what generates a very different skyrmion density map for the two states (see Fig. 7b,c). Another interpretation is that the AF-SkX1 state is closer to SQ (SSF peaks on the $\Gamma$-K line) while the AF-SkX2 is closer to SS (SSF peaks on the K–M line). Therefore, AF-SkX1 and AF-SkX2 can be visualized as originating from the underlying classical spin liquid state with circular pattern in SSF (see Fig. 5f) by lifting the degeneracy in two different ways.

Our main findings are summarized in the form of T vs. $h_z$ phase diagram shown in Fig. 7a. We discover three non-trivial states in our study. A liquid-like short range ordered state existing between the PM and the SQ state at zero magnetic field becomes stable at low temperatures with increasing magnetic field. The AF-SkX1 becomes the ground state near $h_z = 0.14$ which then destabilizes in favour of AF-SkX2 near $h_z = 0.21$. The boundaries seperating different phases were inferred from a combination of temperature dependence of relevant components of SSF, specific heat and topological susceptibility as discussed before. The open circles display the variation of skyrmion density ($\langle {\mathscr {T}} \rangle$) as a function of applied field at low temperatures. Note that the presence of phase boundaries is clearly reflected in the sharp changes in the skyrmion density. The existence of a finite skyrmion density in the CSL state indicates the existence of a isolated skyrmions in this phase when the filaments lengths become of the same order as their width (see Fig. 5b). The sharp jump within AF-SkX2 state does not represent a phase transition as the SSF remains qualitatively identical on two sides of the discontinuity. Inset in Fig. 7a shows magnetization, $M_z$, (blue) and magnetic susceptibility, $\chi _M$, (red) as a function of applied field. Note that the phase changes affect the manner in which magnetization increases with applied field and this gets clearly reflected in the peak structure in $\chi _M$ that exactly matches the indicative phase transitions shown by the skyrmion density variation.

Figure 7

(a) The temperature versus magnetic field phase diagram of the Hamilonian Eq. 7. The different phase boundaries are inferred from the order parameter plots shown in Fig. 4. The right axis shows the variation of skyrmion density $\langle {\mathscr {T}} \rangle$ with external magnetic field $h_z$. Field dependence of magnetization (blue) and that of its derivative (red) are shown in inset. The real-space maps of the skyrmion density for the two antiferromagnetic skyrmion phases (b) AF-SkX1, (c) AF-SkX2.

We observe that the total skyrmion density remains almost unchanged in the $h_z$ window corresponding to the AF-SkX1 state. This suggests that the AF-SkX1 state is highly incompressible, and is similar to the packed skyrmion phase discussed by us in a recent paper⁵⁸. In contrast, the skyrmion density in the AF-SkX2 state gradually decreases upon increasing magnetic field. The step wise reduction of $\langle {\mathscr {T}} \rangle$ in AF-SkX2 is a finite size effect, which can accommodate only particular number of skyrmions with the imposed periodic conditions. In continuum limit it is expected that $\langle {\mathscr {T}} \rangle$ should smoothly go to zero. The qualitative difference between the two antiferromagnetic skyrmion states, AF-SkX1 and AF-SkX2, is clearly observed in the corresponding skyrmion density map plots (see Fig. 7b,c) in terms of the opposite polarity of $\langle {\mathscr {T}} \rangle$ in the skyrmion cores.

Conclusion

Starting from a microscopic Rashba-Hund’s model on a triangular lattice in large Hund’s coupling limit, we derived an effective magnetic model in the insulating limit. A comprehensive Monte Carlo simulation study of the model uncovers a variety of intriguing magnetic phases. In particular, we find two distinct antiferromagnetic skyrmion crystals as the ground states of the model in the presence of external magnetic field. The sublattice resolved spin configuration analysis reveals that the antiferromagnetic skyrmion phases consist of three interpenetrating Néel skyrmion states of usual types observed in many multilayer systems^61,62,63. The effective magnetic model allows us to understand the origin of these two AF-SkX states. Existence of a short-range ordered state characterized by circular diffuse pattern in the SSF serves as the parent of the two SkX states. These two states correspond to two different ways of breaking the degeneracy present in the classical spin liquid. A realization of the model studied here can be achieved in the interfaces and heterostructures of transition-metal oxides (TMOs) along (111) direction. The transition-metal (TM) ions with large moment, such as Mn, Fe etc., are particularly relevant for validating our approximation of a classical spin. The neccessary features to realize our model can be found in many real materials^64,65,66,67. One candidate material is $\hbox {GdI}_2$, in which Gd ions in the $4f^75d^1$ state form a triangular lattice arrangement and electrons from partially filled d bands are coupled to localized f bands^68,69. In a recent study by Chakhalian et al., the complex oxides $\hbox {A}_2$ $\hbox {B}_2$ $\hbox {O}_7$ in [111] directional growth opens up a new route where triangular arrangement of high atomic number transition metal ions induce strong spin-orbit coupling⁷⁰. Other potential candidate materials are $\hbox {LaFeO}_3$, $\hbox {LaMnO}_3$, $\hbox {LaFeO}_3$/$\hbox {LaCrO}_3$⁷¹, ($\hbox {LaMnO}_3$)$_2$/($\hbox {LaScO}_3$)$_4$^72,73 bilayers [111]. In $\hbox {LaMnO}_3$ half-filled $\hbox {t}_{2g}$ electrons and $\hbox {e}_g$ are coupled via Hund’s coupling and the bilayers of $\hbox {LaScO}_3$ provide significant spin-orbit coupling. Given that the low-energy magnetic Hamilonian for a Rashba coupled Mott insulator will have the identical form to what we derived here for the Hund’s model, our results regarding the existence of AF-SkX states are also relevant to Mott-insulators on triangular lattices. Since these two skyrmion states can be easily distinguished based on the spin structure factor, our results provide a clear prediction for their observation in neutron scattering experiments.

Methods

We simulate the spin Hamiltonian Eq. (7) via the Classical Monte Carlo technique based on conventional heat bath method⁷⁴. Periodic boundary conditions are implemented along each direction. Temperature parameter is reduced in small steps starting at high temperature to capture the phase transition from paramagnetic to ordered state. For a given value of T and $h_z$, single spin updates are performed by proposing a new spin configuration from a set of uniformly distributed points on the surface of a unit sphere. Note that the option for picking a completely new orientation for spin reduces the tendency of the system to get stuck in the metastable state. The new configuration is accepted based on the standard Metropolis algorithm^75,76. A Monte Carlo run at each magnetic field and temperature consist of $\sim 1 \times 10^5$ Monte Carlo steps (MCSs) for equilibration and twice the number for calculations of the desired observables. For detailed exploration of parameter space we used lattice size $N = 60 \times 60$, and the stability of results is ensured by simulating sizes up to $N = 120 \times 120$ for some selected parameter values. For simulations in the presence of external magnetic field, we use the field cooled protocol, where the temperature is lowered in the presence of finite external field.

The various phases, obtained via Monte Carlo simulations, can be distinguished from the corresponding real-space spin textures (see Fig. 5). Additionally, we have calculated various physical observables to precisely identify the phase transitions. We calculate the magnetization (M), magnetic susceptibility ($\chi _M$), specific heat ($C_V$) and the topological susceptibility ($\chi _{\mathscr {T}}$)⁷⁷, defined as,

$$\begin{aligned} \begin{aligned} M&= \frac{1}{N} \bigg \langle \sum _i S_i^z \bigg \rangle , \\ \chi _M&= \frac{dM}{dh_z}, \\ C_V&= \frac{d \langle E \rangle }{dT}, \\ \chi _{\mathscr {T}}&= \frac{\langle {\mathscr {T}}^2 \rangle - \langle {\mathscr {T}} \rangle ^2}{NT}. \end{aligned} \end{aligned}$$

(8)

Figure 8

Schematic diagram showing locations of nn sites of a central site in the triangular lattice.

The angular brackets denote the Monte-Carlo average of the quantity, $\langle E \rangle = \frac{1}{N} \langle H_{\text {CSE}} \rangle$, and ${\mathscr {T}}$ denotes the discretized skyrmion density, given as⁴¹,

$$\begin{aligned}&{\mathscr {T}} = \frac{1}{4\pi } \Bigg [ \sum _i A_i^{(12)} \text {sgn} [{\mathscr {L}}^{(12)}_{i}] + A_i^{(45)} \text {sgn} [{\mathscr {L}}^{(45)}_{i}] \Bigg ], \end{aligned}$$

(9)

$$\begin{aligned}&{\mathscr {L}} = \frac{1}{8\pi } \Bigg \langle \sum _i {\mathscr {L}}_i^{(12)} + {\mathscr {L}}_i^{(45)} \Bigg \rangle , \end{aligned}$$

(10)

where, $A_i^{(ab)} = ||(\mathbf{S}_{i_a} - \mathbf{S}_{i}) \times (\mathbf{S}_{i_b} - \mathbf{S}_{i})||/2$ is the local area of the surface spanned by three spins on every elementary triangular plaquette $\mathbf{r}_i,\mathbf{r}_a,\mathbf{r}_b$. Here ${\mathscr {L}}_{i}^{(ab)} = \mathbf{S}_i.(\mathbf{S}_{i_a} \times \mathbf{S}_{i_b})$ is the so-called local chirality and $\mathbf{r}_i, \mathbf{r}_1 - \mathbf{r}_5$ (see Fig. 8) are the sites involved in the calculation of $\langle {\mathscr {T}} \rangle$.

Most importantly, we also compute the component resolved spin structure factor (SSF) to characterize the conventional ordered magnetic phases. The SSF is given by,

$$\begin{aligned} S_f(\mathbf{q})= & \, S^{x}_f(\mathbf{q}) + S^{y}_f(\mathbf{q}) + S^{z}_f(\mathbf{q}), \nonumber \\ S^{\mu }_{f}(\mathbf{q})= & \, \frac{1}{N^2} \bigg \langle \sum _{ij} S^{\mu }_i S^{\mu }_j~ e^{-\mathrm{i}{} \mathbf{q} \cdot (\mathbf{r}_i - \mathbf{r}_j)} \bigg \rangle . \end{aligned}$$

(11)

with $\mu = x,y,z$.