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

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.


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 49 , Operator c † iσ ( c iσ ) creates (annihilates) an electron at site i with spin σ ∈ {↑, ↓} . τ is a vector operator with the three Pauli matrices as components. S i denote the localized spins which we assume to be classical vectors with |S i | ≡ 1 . t, and J H denote the strengths of hopping amplitude, Rashba coupling and Hund's rule coupling, in that order. Assuming the lattice constant to be unity, γ ∈ {a 1 , a 2 , a 3 } denotes the basis unit vector of the triangular Bravais lattice with a 1 =(1,0), a 2 =(1/2, √ 3/2 ) and a 3 =(-1/2, √ 3/2 ). s i is the electronic spin operator. h z is the strength of Zeeman field applied along z axis, and i = √ −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-46,50-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 H at half filling. For large J H , it is useful to work in a site-dependent spin-quantization basis achieved via local SU(2) rotations, given by, 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 { θ i , φ i } specifies the orientation in three dimensions of the local moment S i .
The transformed Hamiltonian takes the form, where σ ∈ {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,γ = t σ σ ′ i,γ + σ σ ′ i,γ have contributions from standard tight-binding hopping integral t and the Rashba SOC . 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 56 . The parallel to anti-parallel hopping contributions g pa i,γ , are given by, The anti-parallel to parallel contributions, g ap ij , are given by, In the above equations, j = i + γ . The general expression for the second order perturbative energy correction, with γ ′ =ẑ ×γ . We also note that the Hamiltonian is written in a general form which is valid for any Bravais lattice. We parameterize by α the relative strength of the Rashba coupling as compared to hopping parameter as t = (1 − α)t 0 and = α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 α = 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 α 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 56 . One consequence of this large degeneracy is the presence of entropically stabilized filamentary domain states at low temperatures.
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 H < 0.1 (see Fig. 1). Note that the large J 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 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 RHM in the large J 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 α → 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 q . The 120° state is characterized by a SSF peak at q 0 = (2π/3, 2π/ √ 3) , and the symmetry related points (see Fig. 3f). For small values of α , the SSF at 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 α , the 120° ground state is destabilized in favour of a single-Q (SQ) spiral state with SSF peak located at ±q for one specific 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 q , monotonically decreases upon increasing α (see Fig. 3b). For α > 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 α . This is characterized by a peak in the SSF at a pair of M points (see Fig. 3i). www.nature.com/scientificreports/ 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 α values, the sharp rise in the order parameter is accompanied by a peak in the specific heat. However, for 0.07 < α < 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 56 . We summarize our results of Monte Carlo simulations in the absence of magnetic field as a T − α 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 56 . 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 < α < 0.67.
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 α = 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 q point). The ordering temperature, as inferred from the SSF at the relevant 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 www.nature.com/scientificreports/ 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 . We further characterize the two multi-Q states with the help of skyrmion density ( T ) and topological susceptibility ( χ 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 χ 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 56 , 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.
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 Ŵ -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 www.nature.com/scientificreports/ 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 ( T ) 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, χ 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 χ M that exactly matches the indicative phase transitions shown by the skyrmion density variation.
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 58 . In contrast, the skyrmion density in the AF-SkX2 state gradually decreases upon increasing magnetic field. The step wise reduction of T 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 T 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 T 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    Fig. 4. The right axis shows the variation of skyrmion density T 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. . In LaMnO 3 half-filled t 2g electrons and e g are coupled via Hund's coupling and the bilayers of 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 74 . 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 ∼ 1 × 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 × 60 , and the stability of results is ensured by simulating sizes up to N = 120 × 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 ( χ M ), specific heat ( C V ) and the topological susceptibility ( χ T ) 77 , defined as, The angular brackets denote the Monte-Carlo average of the quantity, �E� = 1 N �H CSE � , and T denotes the discretized skyrmion density, given as 41 , (9) T = 1 4π i A (12) i sgn[L (12) i ] + A (45) i sgn[L (45) i ] ,  = ||(S i a − S i ) × (S i b − S i )||/2 is the local area of the surface spanned by three spins on every elementary triangular plaquette r i , r a , r b . Here L (ab) i = S i .(S i a × S i b ) is the so-called local chirality and r i , r 1 − r 5 (see Fig. 8) are the sites involved in the calculation of T .
Most importantly, we also compute the component resolved spin structure factor (SSF) to characterize the conventional ordered magnetic phases. The SSF is given by, with µ = x, y, z.