RKKY-dipolar Interactions and 3D Spin Supersolid on Stacked Triangular Lattice

Abstract

Inspired by the recent discovery of metallic spin supersolidity and its giant magnetocaloric effect in the rare-earth alloy EuCo$_2$Al$_9$ [Nature 651, 61 (2026)], we perform a combined study through electronic structure analysis, effective spin model, and Monte Carlo simulations on a stacked triangular lattice, and reveal a novel mechanism for the emergence of 3D spin supersolid in a metallic antiferromagnet. From first-principles inputs, we derive a minimal spin model on a stacked triangular lattice (STL), which arises from the interplay between Ruderman-Kittel-Kasuya-Yosida (RKKY) and dipolar interactions and accurately reproduces the experimental thermodynamics. Based on the STL model, we identify a ground state that simultaneously breaks discrete lattice translational symmetry and continuous spin-rotational symmetry -- the hallmark of a spin supersolid. Furthermore, we present the field-temperature phase diagram of the 3D STL model and discuss the various magnetic phases and associated phase transitions. Under zero field, the spin supersolid Y order establishes in two steps: an upper transition at $T_{N1}$, where an emergent U(1) symmetry appears and the system enters a fluctuating collinear regime, followed by a lower transition at $T_{N2}$ into the spin supersolid Y phase. In contrast, the supersolid V phase undergoes a single phase transition at $T_N^V$. Our results not only provide a comprehensive theoretical understanding of the metallic spin supersolid reported for EuCo$_2$Al$_9$ but also pave the way for further experimental investigations into its supersolid transitions and universality class.

Figures (5)

• Figure 1: $|~$ Electronic structure and localized magnetic moments in EuCo$_2$Al$_9$.(a), Stacked triangular lattice of Eu$^{2+}$ ions in EuCo$_2$Al$_9$ (ECA) with interlayer distance ($r_c = 3.901$ Å) much shorter than the intralayer one ($r_a = 7.871$ Å). The red arrows refer to the spins on site $i$ and $j$, and the unit vector $\hat{\mathbf{r}}_{ij}$ is also indicated. (b), Projected density of states (PDOS) of Eu, Co, and Al. (c-d), Fermi surfaces of ECA at (c)$\Gamma$-$K$-$A$ surface and (d)$\Gamma$-$K$-$M$ surface. The estimated nodes based on the free-electron approximation along the $k_c$ direction ($\Gamma$–$A$) and $k_a$ direction ($\Gamma$–$K$) are marked in light gray and light blue, respectively. (e), Schematic Doniach phase diagram of the Kondo lattice. The blue region where $T_{\text{RKKY}} > T_{\text{Kondo}}$ indicates the magnetic-ordered regime, and the red region where $T_{\text{RKKY}} < T_{\text{Kondo}}$ represents the heavy-fermion regime. The boundaries of the two regimes are given by $|T_{\text{RKKY}} - T_{\text{Kondo}}|$. The ECA is placed deeply at the $T_{\text{RKKY}} > T_{\text{Kondo}}$ regime. The inset shows the RKKY interaction profile based on the free-electron approximation.
• Figure 2: $|~$Parameter fitting and comparison between calculated and experimental results. (a) We fit $J_R(r_a)$ and $J_R(r_c)$ using a loss function involving the out-of-plane magnetization data (see main text), obtaining optimal parameters $J_R(r_c) = -0.34 \pm 0.04$ K and $J_R(r_a) = 0.152 \pm 0.003$ K. (b) Box plot of optimal parameters. The upper and lower error bars indicate the range of parameters with loss below 0.01, while the top and bottom edges of the box denote the standard deviation of the parameter distribution. The central line within each box represents the mean. Results are shown for the STL model parameters (left) and the RKKY and dipolar interaction parameters (right) for comparison. (c-h) Comparison between calculated (d, f, h) and experimental (c, e, g) results: (c-d)$c$-axis magnetization, (e-f)$a$-axis magnetization, and (g-h) specific heat for $B\parallel c$. The inset in (h) compares the experimental and calculated upper transition temperature $T_{N1}$, determined from the prominent specific heat peak, as a function of magnetic field $B$.
• Figure 3: $|~$ Calculated phase diagram and order parameters. (a) The simulated phase diagram of the minimal STL model. The blue and khaki open circles denote the phase boundaries, determined from the peak positions of $\langle |O_z|^2 \rangle - \langle |O_z| \rangle^2$ and $\langle |O_{xy}|^2 \rangle - \langle |O_{xy}| \rangle^2$, respectively. (b,d) Temperature dependence of the order parameters $|O_z|$ (blue dots) and $|O_{xy}|$ (orange dots) at magnetic fields (b)$B = 0$ T and (d)$B = 2.2$ T. (c,e) Log-log plot of $|O_z|$ and $|O_{xy}|$ as functions of $T - T_N$ at (c)$B = 0$ T and (e)$B = 2.2$ T. The red line in panels (b–e) represents the critical behavior of the 3D XY universality class, $(T - T_N)^\beta$, with $\beta \simeq 0.3487$. In panel (e), the purple solid line corresponds to a power-law fit with $\beta = 0.78$. (f, g) Distribution of the order parameter $O_z$ in the CLO$^*$ phase ($T = 2.4$ K) and near $T_{\text{N}1}$ ($T = 3.4$ K) at (f)$B = 0$ T and (g)$B = 0.05$ T. The calculations are performed with parameters $J_R(r_c) = -0.345$ K, $J_R(r_a) = 0.150$ K, and on $24 \times 24 \times 24$ lattice.
• Figure 4: $|~$Calculated results of the full STL model with all long-range dipolar interaction included.(a-b) Calculated distribution of the (a) in-plane magnetization $m_{xy}$ and (b) out-of-plane magnetic order parameter $O_z$ for the Y state within the full model at $B=0$ T and $T=0.5$ K. (c) Corresponding distribution of the phase angles $\theta$ of $O_z$ and $m_{xy}$. (d) Calculated magnetization curves for the minimal model and the full model at $T=1$ K . (e-h) Schematic diagrams of four magnetic configurations observed in ECA with a magnetic field $B$ applied along the $c$-axis: (e) Y state, (f) UUD state, (g) V state, and (h) paramagnetic (PM) state.
• Figure 5: $|~$ Field dependence of the order parameters at a fixed temperature $T = 0.5$ K. The phase transitions from Y to UUD, UUD to V, and V to PM occur at transition fields $B_{c1} \simeq 0.28$ T, $B_{c2} \simeq~2.12$ T, and $B_{c3}~\simeq 2.46$ T, respectively. A clear jump in the both $O_{z}$ and $O_{xy}$ occurs at the UUD-V phase transition around $B_{c2}$.