Multiple Incommensurate Magnetic States in the Kagome Antiferromagnet Na₂Mn₃Cl₈

A polycrystalline sample (mass $2.1$ g) of Na₂Mn₃Cl₈ was prepared by sealing a stoichiometric mixture MnCl₂ and NaCl in SiO₂ after heating at 250 ^∘C under dynamic vacuum overnight. The mixture was heated to 750 ^∘C for several hours and quenched by removing from the furnace. The sample was ground, sealed with $\frac{1}{4}$-atm argon, and annealed at $350$ ^∘C for at total of $\approx 260$ h with an additional intermediate grinding. All handling of this very air-sensitive sample was conducted in an inert-atmosphere glovebox, and the samples were kept under inert atmosphere when they were transferred to the vacuum lines for sealing of the silica tubes.

Magnetization measurements were performed using Quantum Design magnetometers with data below $1.8$ K collected using a ³He insert. The samples were loaded into measurement straws in an inert-atmosphere glovebox with grease to protect the powders from air during the rapid loading process.

Neutron-diffraction measurements were performed using the HB-2A powder diffractometer at the High Flux Isotope Reactor of Oak Ridge National Laboratory. The incident neutron wavelength $\lambda=2.4109\,\text{\AA}$. Our powder sample of mass $2.1$ g was loaded into a 4-mm-diameter cylindrical Cu container in a He glovebox. The sample was cooled using a cryostat with a ³He insert, affording a base temperature of $\approx 0.3$ K. Counting times were $\approx 3$ hr at 0.3, 0.8, 2.0, 5.0, and 40 K, and $\approx 0.5$ hr at other temperatures below $T_{N1}$. The data were corrected for neutron absorption by the sample (Hewat_1979, ).

Magnetic diffuse-scattering refinements were performed using the Spinteract program to refine the values of the exchange interactions (Paddison_2022a, ). The spin Hamiltonian included Heisenberg exchange interactions and the magnetic dipolar interaction (see Section III.5). The input data were collected at $2$ K and $5$ K and were placed in absolute intensity units (barn sr^-1 Mn^-1) by normalization to the nuclear Bragg profile. A high-temperature ($40$ K) data set was subtracted from these data.

The magnetic diffuse scattering $I(Q)$ and bulk susceptibility $\chi T$ were calculated using Onsager reaction-field theory (Paddison_2022a, ; Logan_1995, ; Brout_1967, ), and a $40$ K calculation was subtracted from the calculated $I(Q)$. In this approach, the Fourier transform of the magnetic interactions is calculated as $J_{ij}^{\alpha\beta}(\mathbf{q})=\sum_{\mathbf{r}}J_{ij}^{\alpha\beta}(\mathbf{r})\exp(-\mathrm{i}\mathbf{q}\cdot\mathbf{r})$, where $\alpha,\beta$ denote Cartesian spin components, $i,j\in\{1,3\}$ denote sites within the primitive unit cell, and $\mathbf{r}$ is the vector connecting unit cells containing sites $i$ and $j$. The interaction matrix formed by the $J_{ij}^{\alpha\beta}(\mathbf{q})$ is diagonalized on a grid of up to $50^{3}$ points in the first Brillouin zone to determine its eigenvalues $\lambda_{\mu}(\mathbf{q})$ and eigenvector components $U_{i\mu}^{\alpha}(\mathbf{q})$,

where $\mu\in\{1,3\}$ indexes the normal modes. The long-range dipolar interaction is included using Ewald summation (Enjalran_2004, ). Within a reciprocal-space mean-field approximation, the magnetic propagation vector of the first ordered state is the wavevector at which $\lambda_{\mu}$ reaches its maximal value. The $I(Q)$ and $\chi T$ are given in terms of the $\lambda_{\mu}(\mathbf{q})$ and $U_{i\mu}^{\alpha}(\mathbf{q})$, as described in Ref. (Paddison_2022a, ).

During the refinements, we minimized the function

where subscript “expt” and “calc” indicate measured and calculated diffuse scattering patterns, respectively, $\sigma$ is an experimental uncertainty, and $s$ is a refined overall scale factor common to the neutron-scattering data and the magnetic susceptibility $\chi T$. The minimization was performed using the Minuit program (James_1975, ; James_1994, ). To identify local minima in $\chi^{2}$, we performed $25$ refinements for each model, with different randomly-chosen initial parameter values in each case.

Rietveld refinements were performed using the Fullprof software (Rodriguez-Carvajal_1993, ; Rodriguez-Carvajal_1993a, ). A crystal-structure refinement was first performed at $T=2$ K ($>T_{N1}$). In addition to the crystallographic parameters given in Table 1, we refined the intensity scale factor, $2\theta$ zero-offset, peak-shape, and background parameters. The peak shape was modeled using a pseudo-Voigt function initialized with the instrument resolution parameters, with $U$, $V$, and $W$ parameters subsequently refined. The background was fitted using Chebychev polynomials.

Magnetic Rietveld refinements were performed against $0.3$ K and $0.8$ K data from which the $2$ K data had been subtracted. This subtraction isolates the magnetic Bragg signal by subtracting the nuclear and background contributions, which are essentially unchanged between $0.3$ and $2$ K. In the magnetic refinements, asymmetry, Chebychev background, and magnetic-structure parameters were refined, as described in Section III.4; all other parameters were fixed at the values obtained from the $2$ K refinement. Magnetic-structure figures were prepared using the Vesta program (Momma_2008, ).

Density functional theory calculations were performed using the all-electron-density functional code Wien2K (Sjostedt_2000, ; Blaha_2001, ). The linearized augmented plane wave method (Singh_2006, ) and the generalized-gradient approximation of Perdew, Burke, and Ernzerhof (Perdew_1996, ) were utilized. The $RK_{\mathrm{max}}$ generated by the smallest linearized augmented plane wave sphere radius ($R$) and the interstitial plane-wave cutoff ($K_{\mathrm{max}}$) was set as 7.0 for good convergence. The muffin-tin radii of Na, Cl, and Mn atoms were $2.47$ a.u., $2.14$ a.u., and $2.49$ a.u., respectively. The number of $\mathbf{q}$-points in the full Brillouin zone was $200$. Lattice parameters of Na₂Mn₃Cl₈ were fixed to the experimental values of $a=b=7.423$ Å and $c=19.497$ Å. Then, the internal atomic coordinates were relaxed until forces on all of the atoms were less than $1$ mRy/bohr, with non-magnetic, ferromagnetic, and interlayer antiferromagnetic states. It turns out the relaxed crystal structure with a ferromagnetic state is highly similar to the experimental crystal structure. However, in the non-magnetic state, the atomic coordination changes significantly, as Cl atoms moves towards the Mn layers. The relaxed antiferromagnetic crystal structure is same as the ferromagnetic crystal structure, but exhibits lower energy and smaller forces. The energy difference between ferromagnetic and antiferromagnetic states is $1.19$ meV/f.u.. We therefore used the relaxed antiferromagnetic structure to calculate the intralayer and interlayer magnetic couplings.

The crystal structure of Na₂Mn₃Cl₈ is shown in Figure 1(a), and comprises of triangular Na⁺ layers separating kagome Mn²⁺ layers (Loon_1975, ; Devlin_2022, ). We performed Rietveld refinements against our $2$ K and $40$ K neutron-diffraction data to investigate the possibility of a crystallographic distortion from the published structure (space group $R\bar{3}m$). Good agreement was obtained with the published structural model (Loon_1975, ) at both temperatures, except for two very weak peaks at $1.64$ and $2.16\,\text{\AA}^{-1}$ that were not accounted for, and were unchanged between $0.3$ and $40$ K. Since these peaks could not be explained by simple multiples of the crystallographic unit cell, or by possible impurity phases (NaCl, MnCl₂, or NaMnCl₃), we concluded that the sample or its environment contained a small fraction of unknown impurity. Our results do not show evidence for any structural phase transition between $2$ K and $40$ K, indicating that the broad $\sim$$6$ K specific-heat anomaly reported previously (Devlin_2022, ) is probably due to magnetic ordering of a minor NaMnCl₃ impurity phase (a possibility noted in Ref. (Devlin_2022, )).

Each nearest-neighbor Mn–Mn bond is bridged by a Cl1 ion and a Cl2 ion, which provide the nearest-neighbor superexchange pathways. The Mn–Cl1–Mn and Mn–Cl2–Mn bond angles are $92.42^{\circ}$ and $94.03^{\circ}$, respectively. Since these values are close to $90^{\circ}$, the Goodenough-Kanamori rules predict weak ferromagnetic nearest-neighbor exchange interactions. Further-neighbor interactions have more complicated pathways and, consequently, are difficult to predict. In particular, there are two inequivalent third-neighbor exchange pathways with the same interatomic separation [Figure 1(c)].

Our high-temperature bulk magnetic susceptibility measurements and Curie-Weiss fits are shown in Figure 2(a). They reveal an effective magnetic moment of $5.99$ $\mu_{\mathrm{B}}$, close to the spin-only value of $5.92$ $\mu_{\mathrm{B}}$ for Mn²⁺, and a Weiss temperature of $\theta_{\mathrm{W}}=-4.6(1)$ K, indicating net antiferromagnetic interactions. An anomaly is observed at $T_{N1}\approx 1.6$ K in our low-temperature magnetic susceptibility data, which is suppressed to lower temperature with increasing applied field, consistent with a magnetic ordering transition [Figure 2(b)]. The “frustration parameter”, $f=\theta_{\mathrm{W}}/T_{N1}\approx 3$, indicates a relatively small degree of frustration. Since the nearest-neighbor kagome antiferromagnet is highly frustrated, this result hints at the presence of significant further-neighbor couplings or anisotropies; however, the nature of these couplings cannot be determined from bulk characterization data alone. Interestingly, and consistent with Ref. (Devlin_2022, ), we also observe a second magnetic-susceptibility anomaly at $T_{N2}\approx 0.6$ K [Fig. 2(b)], suggesting a multi-stage magnetic ordering process. Such behavior is unusual and hints that, despite the relatively small value of $f$, the frustrated topology of the kagome lattice may cause several magnetic structures to be nearly degenerate. Figure. 2(c) shows the low temperature field dependence of the magnetization, which does not follow the Brillouin function, in qualitative agreement for theoretical predictions for the kagome lattice (Nakano_2015, ). Figure 1(a) shows the field derivative of the magnetization, $dM/dH$. Several anomalies are observed in $dM/dH$ below $T_{N2}$ at small applied fields, suggesting that the magnetic ground state is fragile to external perturbations.

We performed neutron-diffraction measurements to obtain microscopic insight into the magnetic interactions and structures of Na₂Mn₃Cl₈ (see Methods). An overview of the temperature dependence of our neutron data is shown in Figure 3(a). Several new Bragg peaks appear below $T_{N1}\approx 1.6$ K, most prominently at wavevectors of approximately $0.6$ and $1.5\,\text{\AA}^{-1}$. We identify these as magnetic Bragg peaks arising from the onset of long-range magnetic ordering, since they appear at the same temperature as the magnetic-susceptibility anomaly at $T_{N1}$. Interestingly, the positions of the magnetic Bragg peaks suddenly shift at $T_{N2}\approx 0.6$ K, revealing that the second phase transition involves a change in magnetic propagation vector. At temperatures above $T_{N1}$, broad magnetic diffuse scattering features can be seen, indicating the development of short-range magnetic correlations as $T_{N1}$ is approached from above. We discuss the Bragg and diffuse magnetic scattering in Section III.4 and III.5, respectively.

We first discuss possible ordered magnetic structures of Na₂Mn₃Cl₈, as determined by analyzing the magnetic Bragg profiles obtained at temperatures below $T_{N1}$.

We used the program KSearch of the Fullprof suite (Rodriguez-Carvajal_1993, ; Rodriguez-Carvajal_1993a, ) to identify possible propagation vectors at $0.3$ K ($T<T_{N2}$) and $0.8$ K ($T_{N2}<T<T_{N1}$). The positions of 10 magnetic peaks (at $0.8$ K) and 15 magnetic peaks (at $0.3$ K) were provided as input, and a systematic search of candidate propagation vectors $\mathbf{q}=q_{x}\text{{a}}^{\ast}+q_{y}\text{{b}}^{\ast}+q_{z}\text{{c}}^{\ast}$ was performed, starting with those that lie on a symmetry point, line, or plane of the Brillouin zone. However, none of the high-symmetry propagation vectors was compatible with the observed Bragg positions, at either temperature. The best-fit propagation vectors were instead of the form $(q+\delta,q-\delta,\frac{3}{2})$, with $\delta\ll q$. We obtain $(q,\delta)_{\mathrm{}}\approx(0.29,0.02)$ at $0.8$ K, and $(q,\delta)\approx(0.27,0.06)$ at $0.3$ K; precise values are given in Table 2. These propagation vectors lie on a general position, but they are close to the high-symmetry $(q,q,\frac{3}{2})$ plane.

Having determined possible propagation vectors, we used the program Sarah (Wills_2000, ) to identify symmetry-allowed magnetic structures. The primitive unit cell contains three Mn²⁺ sites, with fractional coordinates $\text{{r}}_{1}=(\frac{1}{2},0,\frac{1}{2})$, $\text{{r}}_{2}=(0,\frac{1}{2},\frac{1}{2})$, and $\text{{r}}_{3}=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ with respect to the conventional axes $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$. Each site has three magnetic degrees of freedom, which are not further constrained by symmetry. We choose these as basis-vector components along orthonormal axes $\mathbf{q}_{\parallel}$, $\mathbf{q}_{\perp}$, $\hat{\mathbf{c}}$, where $\hat{\mathbf{c}}$ is parallel to the $c$-axis, $\mathbf{q}_{\parallel}$ is parallel to the projection of $\mathbf{q}$ in the $ab$-plane, and $\mathbf{q}_{\perp}=\hat{\mathbf{c}}\times\mathbf{q}_{\parallel}$ is perpendicular to $\mathbf{q}_{\parallel}$ and c. These structures are amplitude-modulated spin-density waves (sine structures), with different spin magnitudes and orientations for each site,

where $\mathbf{q}$ denotes the propagation vector, $\mathbf{R}$ denotes a lattice vector, $j\in\{1,3\}$ labels sites within the unit cell, and $\mu_{\mathbf{q}_{\parallel}},\mu_{\mathbf{q}_{\perp}},\mu_{\mathbf{c}}$ are basis-vector components. Alternatively, it is possible to construct helical structures such as

where, in this case, the spin plane is perpendicular to the $c$-axis. The ordered magnetic-moment length can be identical on all sites in the crystal in a helical structure, for example if $\mu_{\mathbf{q}_{\parallel}}=\mu_{\mathbf{q}_{\perp}}$ in Eq. (3).

Due to the relatively large number of variable parameters and the limitations of powder data, we make two assumptions when testing candidate magnetic structures. First, we only consider structures that order with a single propagation vector (single-$\mathbf{q}$ structures). While multi-$\mathbf{q}$ structures are possible, they cannot generally be distinguished from single-$\mathbf{q}$ structures by powder diffraction Kouvel_1963 . Second, we initially assume that the basis vectors at sites $\mathbf{r}_{1}$, $\mathbf{r}_{2}$ and $\mathbf{r}_{3}$ are parallel; this assumption is reasonable because the interactions between nearest and next-nearest neighbors are ferromagnetic, as we will show in Section III.5. Magnetic-structure models were tested against the magnetic Bragg profile using Rietveld refinement (see Methods).

We first considered amplitude-modulated sine structures. The assumption of parallel basis vectors reduces the number of refined parameters from 9 to 3. Sine structures yield excellent agreement with our data at both 0.3 and 0.8 K, as shown in Figure 4(a) and (e), respectively. The magnetic moment is predominantly oriented along $\mathbf{q}_{\perp}$ for the corresponding structures, which are shown in Figure 4(c) and (g), respectively. The refined parameter values and goodness-of-fit metric $R_{\mathrm{wp}}$ are given in Table 2. To determine if sine structures are physically reasonable, we calculated the maximum value of the ordered magnetic moment, $\mathrm{max}(\mu_{\mathrm{ord}})$. For a spin-only ion, this value should not normally exceed $2S\thinspace\mu_{\mathrm{B}}$ ($=5.0\thinspace\mu_{\mathrm{B}}$ for Mn²⁺). This expectation is confirmed by the low-temperature magnetization of Na₂Mn₃Cl₈, which saturates to approximately $5\thinspace\mu_{\mathrm{B}}$ per Mn²⁺ [Figure 2(d)]. Unfortunately, we find $\mathrm{max}(\mu_{\mathrm{ord}})\gg 5.0\mu_{\mathrm{B}}$ for the refined sine structures: $\mathrm{max}(\mu_{\mathrm{ord}})=6.15(7)\thinspace\mu_{\mathrm{B}}$ at $0.3$ K, and $5.40(7)\thinspace\mu_{\mathrm{B}}$ at $0.8$ K. These values are physically unreasonable, suggesting that the correct structures of Na₂Mn₃Cl₈ are not single-$\mathbf{q}$ sine structures.

Circular helices are promising alternative structures, since all sites have equal magnetic moment lengths. Initially, we consider circular helices with magnetic moments in either the $\mathbf{q}_{\parallel}\mathbf{c}$ plane, the $\mathbf{q}_{\perp}\mathbf{c}$ plane, or the $\mathbf{ab}$ plane (equivalent to the $\mathbf{q}_{\perp}\mathbf{q}_{\parallel}$ plane). At both $0.3$ and $0.8$ K, the best fit is obtained for the $\mathbf{ab}$-helix, with slightly worse agreement for the $\mathbf{q}_{\perp}\mathbf{c}$-helix [Table 2]. The $\mathbf{q}_{\parallel}\mathbf{c}$-helix yields much worse agreement than the other structures, so we do not consider it further. The fits for $\mathbf{ab}$-helices at $0.3$ and $0.8$ K are shown in Figure 4(b) and (f), respectively, and the corresponding structures are shown in Figure 4(d) and (h). The agreement with the data is very good, although marginally worse than for the corresponding sine structures. Importantly, however, the refined values of $\mu_{\mathrm{ord}}$ are now physically reasonable, with a maximum value of $4.77(5)\thinspace\mu_{\mathrm{B}}$ at $0.3$ K. This result favors the helical structures.

We tested two variations of the helical structures in an effort to improve the fit quality. First, we considered the $\mathbf{ab}$-helix and relaxed our previous assumption of parallel basis vectors, by refining a clockwise rotation $\Delta\phi$ of the basis vector at position $\mathbf{r}_{3}=(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ about the $c$-axis. The optimal fit is obtained for relatively small $\Delta\phi\approx 20^{\circ}$ at $0.3$ K, and substantial $\Delta\phi\approx 90^{\circ}$ at $0.8$ K. Second, we maintain the assumption that the basis vectors are parallel, but vary the spin plane as $\mathbf{q}_{\perp}(\mathbf{q}_{\parallel}\cos\theta+\mathbf{c}\sin\theta)$. At $0.3$ K, a minimum in $R_{\mathrm{wp}}$ occurs for $\theta\approx 35^{\circ}$, whereas at $0.8$ K, fit quality is essentially unchanged for all $\theta\lesssim 30^{\circ}$. Each of these variations yields a similar or slightly improved fit compared to the simple $\mathbf{ab}$-helix (see Table 2).

Figure 5 shows the temperature evolution of the refined parameter values for the $\mathbf{ab}$-helix with collinear basis vectors. A discontinuity in the propagation vector is apparent at $T_{N2}$ [Figure 5(a)]. No such anomaly is apparent in the refined value of the ordered magnetic moment, which increases smoothly on cooling the sample below $T_{N1}$ [Figure 5(b)]. The temperature dependence of this order parameter at low temperatures ($T\leq 1.2$ K) is consistent with the phenomenological form $\mu_{\mathrm{ord}}\propto 1-cT^{2}$ for a three-dimensional magnet with half-integer spin Kobler_2003 . Its temperature dependence for $1.2\leq T\leq 1.6$ K is consistent with critical form for a three-dimensional Heisenberg magnet, $\mu_{\mathrm{ord}}\propto(T_{N1}-T)^{0.365}$, although the small number of data points precludes fitting the critical exponent.

In conclusion, our powder-diffraction data are well explained by circular helical magnetic structures. Basis vectors are close to the $\mathbf{ab}$ plane and nearly collinear at $0.3$ K, with a possibility of greater noncollinearity at $0.8$ K. We emphasize, however, that the possibility of multi-$\mathbf{q}$ structures cannot be ruled out, and we discuss this further in Section IV.2.

We seek to parametrize the spin Hamiltonian of Na₂Mn₃Cl₈ by analyzing the diffuse magnetic scattering measured above $T_{N1}$. This approach is an alternative to spin-wave analysis of inelastic neutron-scattering data, and has recently been applied to several frustrated antiferromagnets (Samarakoon_2020, ; Scheie_2021, ; Welch_2022, ). The magnetic diffuse scattering measured at $2$ K and $5$ K (with $40$ K data subtracted) is shown in Figure 6. Diffuse magnetic peaks are sharper at $2$ K than at $5$ K, consistent with an increase in the magnetic correlation length on cooling the sample. The bulk magnetic susceptibility expressed as $\chi T$ is also shown in Figure 6 and confirms the development of antiferromagnetic correlations.

To model these data, we consider a Hamiltonian that includes Heisenberg exchange interactions and the long-range magnetic dipolar interaction,

where $\mathbf{S}_{i}$ is modeled as a classical vector of magnitude $\sqrt{S(S+1)}$, $S=5/2$ is the spin quantum number of Mn²⁺, $\hat{\mathbf{r}}_{ij}=|\mathbf{r}_{j}-\mathbf{r}_{i}|/r_{ij}$ is a unit vector parallel to the separation of spins $i$ and $j$, and $r_{1}=3.7124(1)\,\text{\AA}$ is the nearest-neighbor distance. The exchange interactions include the nearest-neighbor exchange $J_{1}$, the inter-layer coupling $J_{c}$, and the further-neighbor couplings shown in Figure 1(c), so that $J_{ij}\in\{J_{1},J_{2},J_{3a},J_{3b},J_{c}\}$. The magnitude of the dipolar interaction at the nearest-neighbor distance, $D=\mu_{0}(g\mu_{\mathrm{B}})^{2}/4\pi r_{1}^{3}k_{\mathrm{B}}=0.0487$ K, is determined by the crystal structure. The values of the exchange interactions were optimized against our $I(Q)$ and $\chi T$ data (see Methods).

We tested interaction models against our data in order of increasing number of exchange parameters, as follows. Unless otherwise noted, the dipolar interaction was fixed at $D=0.0487$ K. For each exchange model, the best fit to $I(Q)$ and $\chi T$ data is shown in Figure 6. The refined values of the exchange interactions are shown in Table 3, along with the goodness of fit metric $R_{\mathrm{wp}}$, and two quantities estimated from the Onsager-reaction-field calculation that may indicate model quality: the predicted magnetic ordering temperature $T_{N}^{\mathrm{calc}}$ and the predicted magnetic propagation vector $\mathbf{q}_{\mathrm{calc}}$.

First, we considered a minimal model in which only $J_{1}$ and $J_{c}$ were refined [model (a)]. This model does not represent our neutron data well [Figure 6(a)], and the predicted propagation vector is commensurate, in contrast to the incommensurate propagation vector observed experimentally. Second, we refined $J_{1}$, $J_{2}$, and $J_{c}$ parameters [model (b)]. This model yields a substantially improved fit, but some misfit is still evident in the $I(Q)$ and, especially, the $\chi T$ data [Figure 6(b)]. The calculated propagation vector is now incommensurate, but different to the experimental one. Third, we refined $J_{1}$, $J_{2}$, $J_{3a}$, $J_{3b}$, and $J_{c}$ parameters [model (c)]. This model yields an excellent fit to both the $I(Q)$ and $\chi T$ data [Figure 6(c)]. Moreover, the calculated propagation vector, $(0.30,0.30,\frac{3}{2})$, is close to the experimental value of $(0.3102(4),0.2646(4),\frac{3}{2})$ in the first ordered state at $0.8$ K, and the calculated $T_{N}^{\mathrm{calc}}\approx 1.6$ K agrees with the measured value. This refinement was stable despite the relatively large number of free parameters; no large parameter covariances ($\sigma_{ij}\geq 80\%$) were noted, and initializing the refinement with different parameter values yielded only one possible local minimum, which had significantly worse $R_{\mathrm{wp}}^{\mathrm{neutron}}=23.0\%$ and $R_{\mathrm{wp}}^{\chi T}=4.1\%$.

Our results suggest that model (c) represents well the interactions of Na₂Mn₃Cl₈. This model has weak ferromagnetic $J_{1}$, consistent with the Goodenough-Kanamori rules. The inter-layer coupling $J_{c}$ is antiferromagnetic, consistent with the antiferromagnetic layer stacking observed below $T_{N1}$. The third-neighbor couplings $J_{3a}$ and $J_{3b}$ are antiferromagnetic and significantly larger than $J_{1}$. Hence, Na₂Mn₃Cl₈ is an unusual system where strong antiferromagnetic third-neighbor interactions compete with ferromagnetic nearest-neighbor interactions. To the best of our knowledge, a similar $J_{1}$-$J_{3}$ competition has been identified in only one other kagome material, vesignieite Boldrin_2018 . However, this material differs from Na₂Mn₃Cl₈ in its magnetic properties as well as its chemistry, as it has $S=1/2$ and shows commensurate magnetic ordering Boldrin_2018 .

Finally, we considered the relevance of the long-ranged dipolar interaction by performing a fourth refinement in which $J_{1}$, $J_{2}$, $J_{3a}$, $J_{3b}$, and $J_{c}$ parameters were varied, while $D$ was fixed at zero [model (d)]. This refinement yielded worse agreement with $I(Q)$ and $\chi T$ data, and significantly underestimates the value of $T_{N1}$ [Table 3]. This result shows that the dipolar interaction has a significant effect on the magnetic properties, as expected since $D$ is of comparable magnitude to the exchange interactions. However, the refined values of all parameters except $J_{1}$ are equivalent (within $1\sigma$) for models (c) and (d), suggesting that the effect of the dipolar term on these refinements is largely confined to nearest neighbors.

To gain insight into the exchange interactions, we performed first-principles calculations using density-functional theory (see Methods). The values of the interactions calculated using DFT are given in Table 4 for different values of the Hubbard $U$ between $0$ and $5.25$ eV. Based on other materials, we anticipate that $U$ is likely between $4$ and $5.25$ eV.

The first-principles exchange interactions show similarities with the experimentally-determined values, but also substantial differences. On the one hand, the first-principles values of $J_{1}$ and $J_{c}$ are ferromagnetic and antiferromagnetic, respectively, consistent with the values fitted to experimental data. The magnitudes of $J_{1}$ and $J_{c}$ for $U=5.25$ eV are also comparable to the experimentally-determined magnitudes, in contrast to a previous DFT study that reported interactions larger than $30$ K . On the other hand, the first-principles values of $J_{2}$, $J_{3a}$, and $J_{3b}$ are opposite to the experimentally-determined values; moreover, the calculated magnitudes of these interactions are very large compared to the other interactions.

We carefully checked whether the first-principles results could be consistent with our experimental data. Taking $U=5.25$ eV, we calculate the Weiss temperature as $\theta_{\mathrm{DFT}}=\frac{4}{3}S(S+1)[J_{1}+J_{2}+J_{c}+J_{3a}+J_{3b}/2]=3.0$ K. Hence, DFT predicts a ferromagnetic Weiss temperature, which is not consistent with the antiferromagnetic value ($\theta=-4.6(1)$ K) measured experimentally. We also estimate the magnetic ordering temperature to be $4.8$ K, which is much larger than the experimental value of $1.6$ K. Finally, we performed additional refinements to neutron and $\chi T$ data as described in Section III.5, except we constrained the signs of the exchange interactions to be the same as those from DFT, while allowing their magnitudes to refine freely. These refinements yielded $J_{3a}=J_{3b}\approx 0$, essentially reproducing the results of model (b) in Section III.5.

We therefore conclude that the DFT results are not fully consistent with our experimental data, making Na₂Mn₃Cl₈ a model material for benchmarking developments in first-principles calculations. The reason for the inaccuracy of the DFT exchange interactions beyond nearest-neighbors is not yet clear; an interesting possibility is that it may relate to the neglect of the Stoner coupling on the Cl ligand sites, as recently proposed in the related material NaMnCl₃ (Solovyev_2022, ).

In this section, we discuss the origin of the multiple incommensurate ordering transitions in Na₂Mn₃Cl₈, using a combination of field-theoretic and Monte Carlo simulations.

Incommensurate magnetic structures are relatively uncommon in kagome antiferromagnets. For example, to the best of our knowledge, all known jarosite minerals that exhibit long-range order have either $(0,0,0)$ or $(0,0,\frac{3}{2})$ propagation vectors (see (Mendels_2011, ) and references therein). Similarly, commensurate states are observed for many other insulating materials in which the kagome lattice is undistorted or slightly distorted; for example, MgFe₃(OH)₆Cl₂ with $\mathbf{q}=(0,0,\frac{3}{2})$ (Fujihala_2017, ), centennialite CaCu₃(OH)₆Cl${}_{2}\cdot$0.6H₂O (Iida_2020, ), CdCu₃(OH)₆(NO₃)${}_{2}\cdot$0.6H₂O (Ihara_2022, ), Nd₃Sb₃Mg₂O₁₄ (Scheie_2016, ), and Sr-vesignieite SrCu₃V₂O₈(OH)₂ with $\mathbf{q}=(0,0,0)$, $\alpha$-Cu₃Mg(OH)₆Br₂ (Wei_2019, ) and YCu₃(OH)₆Cl₃ with $\mathbf{q}=(0,0,\frac{1}{2})$ (Zorko_2019, ), and Ba-vesignieite BaCu₃V₂O₈(OH)₂ with $\mathbf{q}=(\frac{1}{2},0,0)$ (Boldrin_2018, ; Okamoto_2009, ). By contrast, the distorted-kagome material Ba₂Mn₃F₁₁ is one of the only insulating kagome materials with incommensurate magnetic ordering (Hayashida_2018, ). Incommensurate modulations are more frequently observed in metallic kagome systems, such as Tb₃Ru₄Al₁₂ (Rayaprol_2019, ) and YMn₆Sn₆, the latter of which undergoes an incommensurate-to-commensurate transition on cooling (Ghimire_2020, ).

To understand the preference for kagome magnets to form commensurate structures, and the conditions where incommensurate structures may appear, we use the reciprocal-space mean-field approximation introduced in Section II.3 to investigate the stability of different phases as a function of the interactions $J_{1}$, $J_{3a}$, $J_{3b}$, and $J_{c}$ [Figure 1(c)]. Throughout large regions of this interaction space, the classical ground state is one of the commensurate “regular magnetic orders” described in Ref. (Messio_2011, ). Of the models previously investigated theoretically, the most relevant one to Na₂Mn₃Cl₈ is the $J_{1}$-$J_{3a}$ Heisenberg model studied in Refs. (Grison_2020, ; Li_2022, ). The phase diagram for this model is shown in Figure 7(a), and contains five phases: ferromagnetic layers with antiferromagnetic stacking [$\mathbf{q}=(0,0,\frac{3}{2})$], $\mathbf{q=0}$ antiferromagnet, $\sqrt{3}\times\sqrt{3}$ antiferromagnet [$\mathbf{q}=(\frac{1}{3},\frac{1}{3},\frac{3}{2})$], three-sublattice antiferromagnet [$\mathbf{q}=(0,\frac{1}{2},\frac{1}{4})$], and an incommensurate region. This result reproduces the result of Ref. (Grison_2020, ) for isolated kagome planes, except that we include a small antiferromagnetic inter-layer coupling $J_{c}\rightarrow 0^{-}$ to stabilize three-dimensional ordering.

While the $J_{1}$-$J_{3a}$ phase diagram is relatively complicated, it is nevertheless simpler than our model for Na₂Mn₃Cl₈, which also includes significant $J_{c}$, $J_{3b}$, and dipolar couplings. We therefore extended the $J_{1}$-$J_{3a}$ phase diagram to consider the effects of these additional couplings, which are needed for a full description of our Na₂Mn₃Cl₈ data. Notably, for all models, antiferromagnetic $J_{3a}$ is necessary to stabilize incommensurate ordering with $\mathbf{q}=(q,q,\frac{3}{2})$. In Figure 7(b), we fix ferromagnetic $J_{1}=1$ and consider the phase diagram in the $J_{3a}$-$J_{3b}$ plane for antiferromagnetic $J_{c}\rightarrow 0^{-}$. Nonzero $J_{3b}$ has a dramatic effect on the phase diagram; in particular, including antiferromagnetic $J_{3b}$ extends the stability region of the incommensurate phases observed for antiferromagnetic $J_{3a}$. Figure 7(b)–(d) show the effect of increasing the magnitude of $J_{c}$, the antiferromagnetic interlayer coupling ($J_{c}\rightarrow 0^{-}$, $-0.15$, and $-0.65$, respectively, in the same units as $J_{1}$). The effect of increasing $|J_{c}|$ is to increase further the region of phase space in which incommensurate order is stable within the mean-field approximation. Finally, in Figure 7(e), we show the $J_{3a}$-$J_{3b}$ phase diagram including the long-range dipolar interaction $D/J_{1}\approx 0.55$ appropriate for Na₂Mn₃Cl₈. The inclusion of $D$ has a relatively small effect on the positions of the phase boundaries.

The reciprocal-space mean-field theory provides a useful overview of the phase space, but has several important limitations. First, for a non-Bravais lattice such as kagome, it only determines a lower bound on the energy of the ground state. As discussed in Ref. (Grison_2020, ), for the incommensurate region of the $J_{1}$-$J_{3a}$ phase diagram, a physical spin configuration could not be identified that reached this lower bound; hence, the actual magnetic ground state is uncertain in this region. Second, since this theory considers instabilities of the paramagnetic phase, it predicts only the propagation vector of the first ordered state that develops on cooling; it provides no information about the possibility of multiple phase transitions, as are observed experimentally in Na₂Mn₃Cl₈.

We performed classical Monte Carlo simulations to address these limitations. Since the periodicity of an incommensurate magnetic structure does not “fit” within any finite-sized configuration, finite-size artifacts are encountered, which can be reduced by studying relatively large system sizes. However, the long-ranged nature of the magnetic dipolar interaction makes large system sizes computationally expensive. We therefore consider first an approximation to the full Hamiltonian, Eq. (4), where we simulate the parameters that best describe our diffuse-scattering data [model (c) in Table 3], but truncate the dipolar interaction $D$ at the nearest-neighbor distance; we will call this the “nearest-neighbor dipolar model”. For comparison, we also simulated the same model (c) except with $D=0$. To identify finite-size effects, we considered different system dimensions from $10\times 10\times 4$ hexagonal unit cells ($3600$ spins) to $20\times 20\times 8$ hexagonal unit cells ($28800$ spins). For the $10\times 10\times 4$ and $20\times 20\times 8$ simulations only, we slightly adjusted the model (c) interaction parameters to stabilize $\mathbf{q}=(\frac{3}{10},\frac{3}{10},\frac{3}{2})$ ordering, which is commensurate with the system size; this was achieved by multiplying the best-fit values of $J_{3a}$ and $J_{3b}$ by $0.936$. To investigate the effect of a different system geometry, we defined a orthogonal unit cell with axes $\mathbf{a}_{\mathrm{o}}=\mathbf{a}$, $\mathbf{b}_{\mathrm{o}}=\mathbf{a}+2\mathbf{b}$, and $\mathbf{c}_{\mathrm{o}}=\mathbf{c}$, and performed simulations of $12\times 6\times 4$ and $18\times 9\times 6$ orthogonal unit cells ($5184$ and $17496$ spins, respectively). Simulations were run for up to $4.1\times 10^{6}$ moves per spin at low temperatures, where a single move involved one microcanonical (over-relaxation) update followed by a proposed spin rotation of a randomly-chosen spin, which was accepted or rejected according to the Metropolis criterion. Measurements of the autocorrelation function showed that these conditions allowed the system to decorrelate at all temperatures above $0.1$ K. Simulations including the long-ranged dipolar interaction, implemented using Ewald summation (Wang_2001, ), were also performed for a small system size of $10\times 10\times 4$ hexagonal unit cells, without over-relaxation updates.

Results of our Monte Carlo simulations are shown in Figure 8. For the model with $D=0$, a sharp anomaly indicating a single magnetic phase transition is observed at $\approx 0.9$ K; we do not consider the low-temperature state here. The nearest-neighbor dipolar model shows a more complex temperature evolution. In all our simulations, sharp specific heat anomaly is observed at $\approx 0.9$ K, with a second feature between $1.3$ and $1.6$ K that is resolved as either a single broadened peak or two peaks close in temperature, depending on system dimensions. Hence, unlike the Heisenberg model, the nearest-neighbor dipolar model shows at least two magnetic phase transitions, in qualitative agreement with the experimental data for Na₂Mn₃Cl₈. Properties of the magnetic phases obtained for a model of $18\times 9\times 6$ orthogonal unit cells are shown at $1.4$, $1.0$, and $0.3$ K, in Figure 8(c), (d), and (e) respectively. The phases observed at $1.0$ and $0.3$ K are resolved for all other system sizes and geometries. However, the $1.4$ K phase is not resolved in the $20\times 20\times 8$ simulation, suggesting its appearance for some other system sizes may be a finite-size artifact. The calculated magnetic powder diffraction patterns show remarkably good agreement with our experimental powder-diffraction data, especially at $1.4$ and $1.0$ K [Figure 8(c)–(e)]. Calculations of the single-crystal magnetic diffraction patterns reveal magnetic Bragg peaks corresponding to a single incommensurate wavevector at $1.4$ K, indicating a single-$\mathbf{q}$ magnetic structure at this temperature [Figure 8(c)]. Remarkably, however, the same calculation shows magnetic Bragg peaks corresponding to two wavevectors at $1.0$ and $0.3$ K. The intensity of each wavevector is approximately equal at $1.0$ K but significantly different at $0.3$ K [Figure 8(d) and (e)]. The same effect was observed across all our simulations at $1.0$ and $0.3$ K, suggesting this is likely not an artifact due to domain formation, but instead indicates the formation of double-$\mathbf{q}$ states in the Monte Carlo simulations. Our simulations of the long-ranged dipolar model also suggest a possible change in magnetic structure below approximately $1.0$ K, although a second transition is not clearly resolved in the heat capacity for this small simulation size [Figure 8(b)]. For this model, the magnetic structure is clearly 2-$\mathbf{q}$ only below $1.0$ K.

Our results suggest the enticing possibility that the ordered incommensurate states may, in fact, be multi-$\mathbf{q}$ structures rather than single-$\mathbf{q}$ helices. Given the good agreement of our microscopic model with powder-diffraction data and its correct prediction of multiple phase transitions, this scenario is certainly possible. Further theoretical studies including the long-ranged dipolar interaction would be useful to elucidate the relative stabilities of single-$\mathbf{q}$ and multi-$\mathbf{q}$ states, which may be close in energy.