Stripe-yz magnetic order in the triangular-lattice antiferromagnet KCeS2

Yb- and Ce-based delafossites were recently identified as effective spin-1/2 antiferromagnets on the triangular lattice. Several Yb-based systems, such as NaYbO2, NaYbS2, and NaYbSe2, exhibit no long-range order down to the lowest measured temperatures and therefore serve as putative candidates for the realization of a quantum spin liquid. However, their isostructural Ce-based counterpart KCeS2 exhibits magnetic order below TN = 400 mK, which was so far identified only in thermodynamic measurements. Here we reveal the magnetic structure of this long-range ordered phase using magnetic neutron diffraction. We show that it represents the so-called"stripe-yz"type of antiferromagnetic order with spins lying approximately in the triangular-lattice planes orthogonal to the nearest-neighbor Ce-Ce bonds. No structural lattice distortions are revealed below TN, indicating that the triangular lattice of Ce3+ ions remains geometrically perfect down to the lowest temperatures. We propose an effective Hamiltonian for KCeS2, based on a fit to the results of ab initio calculations, and demonstrate that its magnetic ground state matches the experimental spin structure.


Introduction
In low-dimensional quantum magnets, competing exchange interactions may lead to a strong frustration accompanied by enhanced quantum fluctuations. Ultimately this can prevent the systems from long-range order, and the ground state is supposed to be a quantum spin liquid (QSL) [1][2][3][4]. The theoretically predicted QSL state remains experimentally elusive, and a major ongoing research effort is directed towards the identification of promising QSL candidate materials [5].
A short while ago, a whole new family of compounds with the undistorted triangular lattice of rare-earth ions with little or no structural disorder has been identified [39][40][41][42]. These are ternary compounds with the general chemical formula A + R 3+ X 2− 2 , where A is an alkali metal (e.g. Na, K, Rb, Cs), R is a rare-earth ion (in our case, Ce or Yb), and X is a chalcogen atom (O, S, or Se). They crystallize in the layered delafossite structure depicted in Fig. 1 that derives from the rock-salt structure upon cation ordering, realizing ABAB-type stacking of magnetic triangular-lattice layers [43,44]. Low-temperature thermodynamic measurements reveal no signatures of magnetic ordering down to subkelvin temperatures in NaYbO 2 [39,[45][46][47][48], NaYbS 2 [39,40,49], and NaYbSe 2 [39,[50][51][52], indicating that these Yb-delafossite compounds may possess QSL ground states. Long-range magnetic order in these systems can be achieved in a moderate magnetic field [47,50,53]. On the other hand, their isostructural Ce-based counterpart KCeS 2 develops magnetic order below ∼400 mK already in zero field, according to the specific-heat data [54].
Because Ce 3+ and Yb 3+ ions have one electron or hole in the f shell, respectively, and the CEF schemes of all the mentioned compounds have ground-state doublets that are well separated from the excited CEF states [40,46,48,51,54,55], one would generally expect the spin-1/2 approximation to be equally applicable to these systems. However, they may still differ in terms of the magnetocrystalline anisotropy and have different eigenfunctions of the CEF ground state with respect to their decomposition in the |J, J z basis, which would place them in different regions of the parameter space of a generic effective Hamiltonian proposed for triangular-lattice antiferromagnets [56,57]. To understand these differences, we have investigated the magnetic structure of the low-temperature ordered phase in KCeS 2 and its low-temperature lattice structure using powder neutron diffraction.

Sample preparation
The powder sample of KCeS 2 was prepared through the following synthesis route. Potassium carbonate and cerium dioxide were mixed in a 20:1 molar ratio and thoroughly ground in a porcelain mortar. A glassy carbon crucible was filled with this mixture and placed in a tube furnace. Before heating up to the target temperature, the whole apparatus including a reservoir for carbon disulfide was flushed with argon for 30 min. The mixture was heated up to 1050 • C within 3 hours under an unloaded stream of argon (5 L/h). While dwelling one hour, a stream of argon of 0.5 L/h was used to carry CS 2 into the hot zone to enable the sulfidization. Finally, the apparatus was allowed to cool down to 600 • C within 6 hours and without further control down to ambient temperature under a slight argon stream. The ingot that formed in the glassy carbon crucible was then dissolved in water to set the insoluble KCeS 2 free. Filtering via paper and funnel and washing with water and ethanol produced the sample, mainly composed of intergrown crystals of the target compound and some Ce 2 O 2 S as a minority phase. The formation of oxidic by-products during the formation of rare-earth metal chalcogenides (and other compounds of rare-earth metals) is well known and can be attributed to the high oxygen susceptibility of the rare-earth metal ions. It can only be prevented in a completely oxygen-free setup. We also observed that mechanical manipulations (such as pressing, grinding, etc.) can foster decomposition of KCeS 2 with the formation of the Ce 2 O 2 S phase, therefore we avoided grinding the sample and used the as-synthesized powder in neutron-diffraction measurements.

Powder neutron diffraction
To reveal the spin structure in the low-temperature AFM state of KCeS 2 , we measured neutron powder diffraction at the D1B high-intensity two-axis powder diffractometer at ILL, France. Neutrons with a calibrated wavelength λ = 2.5286 Å were selected with a highly oriented pyrolytic graphite (HOPG) (002) monochromator. The contribution of the instrument to the peaks broadening was determined from the refinement of Na 2 Ca 3 Al 2 F 14 standard sample, while wavelengths were refined using a Si standard. Parasitic diffraction peaks arising from the sample environment were eliminated by a radial oscillating collimator. Our delafossite sample represented rough powder consisting of small crystallites with lateral sizes up to 0.5 mm. This sample was placed in a Cu can to ensure good thermal contact at subkelvin temperatures and was cooled down using a dilution refrigerator. The measurements were carried out while the sample was rotated about the vertical axis within 180 • with a step of 1 • (for the 72 and 600 mK datasets) or 3 • (other temperatures), with subsequent averaging of the datasets to compensate for the preferred orientations of individual crystallites. The data were collected at several temperatures between 72 and 600 mK in zero magnetic field, allowing for sufficient time for temperature stabilization before every data collection. The measurements at base temperature and at 600 mK were measured with higher statistics, corresponding to approximately 13 h acquisition time per temperature for the 72 and 600 mK datasets.
The 600 mK dataset, collected above T N , is plotted in Fig. 2 (a). Apart from the structural reflections from the main KCeS 2 phase, it contains reflections from the paramagnetic impurity phase Ce 2 O 2 S [58,59], which constitutes about 15% of our sample judging from the results of structural refinement but exhibits no magnetic Bragg reflections. There are also two intense peaks from metallic copper that originate from the sample container, which we included as the third phase. Our Rietveld refinement of the neutron data was done in the assumption of the rhombohedral R3m space group, previously reported for KCeS 2 from x-ray structure refinement [54,60]. The agreement between the lattice structures obtained at room temperature and at 600 mK (apart from the changes in lattice constants due to thermal expansion) demonstrates the absence of any symmetry-lowering transitions down to T N . The low-temperature lattice parameters resulting from the refinement of our neutron data are summarized in Table 1.
The magnetic peaks, shown in the 72 mK dataset in Fig. 2 (b) with black arrows, set in below 400 mK -the transition temperature revealed earlier in specific-heat measurements [54]. The difference of the low-and high-temperature datasets, plotted in Fig. 2 (c), contains only magnetic Bragg peaks, whereas the structural scattering from KCeS 2 and the oxysulphide impurity phase is subtracted to zero, which Table 1. Low-temperature unit cell parameters, χ 2 values, and R factors for our refinement of neutron diffraction data on the majority phase KCeS 2 and minority phase Ce 2 O 2 S at T = 72 and 600 mK. Here V is the unit cell volume, z(S) is the z position of S on the 6c Wyckoff site in KCeS 2 , and z(Ce) and z(O) are the positions of Ce and O on the 4h Wyckoff site in Ce 2 O 2 S, respectively.  presents strong evidence for the absence of any structural phase transition associated with the magnetic ordering. The commensurate AFM propagation vector (0 1 2 1 2 ) was determined by the K_Search algorithm in the FullProf Suite based on difference of the two datasets. Representation analysis carried out with SARAh [61] for the R3m space group with the aforementioned propagation vector suggests two possible common irreducible representations: Γ 2 and Γ 4 , which we have checked against the temperature-subtracted I(0.072 K) − I(0.6 K) dataset. Our final result of the magnetic structure refinement, which was done using FullProf Suite [62], is shown in Fig. 2 (c) with a black line. The magnetic signal can be perfectly described only using the Γ 2 irreducible representation. Note that the two peaks marked with red arrows are of nonmagnetic origin, as they result from the imperfect subtraction of strong structural reflections from Cu. The magnetic structure of KCeS 2 is collinear, with spins lying approximately in the ab plane and pointing orthogonally to the nearest-neighbor Ce-Ce bonds, as shown in the inset to Fig. 2 (d). The refinement suggests a statistically insignificant out-of-plane canting of Ce 3+ spins by (5.5 ± 5) • , yet fixing this angle to zero has no significant effect on the quality of the fit, so also the in-plane spin arrangement would be consistent with our data. This type of AFM structure on the triangular lattice is known as "stripe-yz" order among theorists [56,57,63,64]. The stacking of the AFM layers along the out-of-plane direction is illustrated in Fig. 2 (e). Note that all magnetic Bragg peaks in Fig. 2 (c) are well described with our model for KCeS 2 , suggesting that the impurity phase Ce 2 O 2 S develops no magnetic order down to 20 mK and could be therefore itself considered a promising spin-liquid candidate.

Parameters
According to our combined refinement of the magnetic and structural data, the ordered magnetic moment on Ce 3+ at base temperature is approximately 0.32(1)µ B . It decreases monotonically with increasing temperature, as shown in Fig. 2 (d), following order-parameter behavior. Fitting this dependence with the empirical function µ(T ) = A tanh π 2 (T N /T − 1) α , shown with a dotted line, gives a transition temperature T N ≈ 430 mK, which is slightly higher than the value previously observed in thermodynamic measurements [54]. The ordered moment in our data saturates below 300 mK, which may be a consequence of poor temperature stabilization in the powder sample close to the base temperature of the dilution refrigerator. Therefore, it cannot be excluded that we somewhat underestimate the value of the ordered moment at the base temperature. From the structural refinement of neutron diffraction data, we also obtained the temperature dependence of lattice constants, as it is shown in Fig. 3 with solid circles. For comparison, we also show our own and previously published cell parameters [60,65] obtained using x-ray diffraction. Both lattice parameters a and c show conventional (positive) thermal expansion above 100 K, whereas below this temperature the thermal expansion is rather small and not measurable within the experimental error. We compare this result with the c-axis thermal expansion measured by capacitive dilatometry (solid line), carried out using our in-house dilatometer with a sensitivity to relative length changes of about 10 −7 [66]. The specimen consisted of several flakelike single crystals of KCeS 2 stacked along the c axis to achieve a minimum required sample thickness. Therefore, only thermal expansion along the c axis but not in the ab plane could be measured. In good agreement with the neutron diffraction data, the capacitive dilatometry confirms a really small thermal expansion below 120 K, where the relative length change is only 3.5 × 10 −4 . The linear expansion coefficient α stays less than 6 × 10 −6 K −1 in this temperature range and only increases noticeably at higher temperatures. Additionally, nearly no magnetostriction could be detected during field sweeps (not shown), where the relative length change is only about −1 × 10 −5 up to 8 T at T = 2 K. The crossover in the behavior of the thermal expansion, observed consistently in diffraction and dilatometry, occurs more than two orders of magnitude higher in temperature than the magnetic ordering transition, which suggests that magnetostrictive effects become effective far above the Néel temperature due to the strong frustration in the system. Moreover, we observe no magnetostriction anomaly at the magnetic ordering temperature (dashed vertical line) in our neutron data. This is another indication that the AFM state is accompanied with no structural distortions, and the delafossite lattice structure with geometrically perfect triangular-lattice layers remains stable even after the stripe-yz AFM order sets in.

DFT/PBE/PAW level
We start the theoretical investigation at the density-functional theory (DFT/PBE/PAW) level using VASP suit [67][68][69]. Admittedly any single-determinant approximation would be a weak one for the 4f system without special arrangements for f -shell. To deal with the complication due to the 4f shell, at first the f shell in core potential for Ce is considered [69]. The Γ -centered k grid of 4 × 4 × 2 was used for system optimization and charge density convergence with plane-wave basis cutoff of 400 eV. The original rhombohedral structure has been optimized down to 10 −4 eV/Å level with no restrictions on symmetry. The optimized structure preserves the original R3m space group but acquires a slightly expanded unit cell compared to the experimental data, a th = 4.23 Å vs. a exp = 4.22 Å and c th = 22.00 Å vs. c exp = 21.83 Å. While the lattice parameter a remains within the experimental uncertainty from the measured value, the c lattice constant experiences a significant (0.78%) expansion along the c axis, which accounts for ∼0.026 Å expansion of the CeS 6 trigonal antiprism along the C 3 axis. Nevertheless, the electronic structure is only slightly affected by these differences. Figure 4 shows the band structure profiles along the high symmetry lines of KCeS 2 with 4f shell in the core and explicit 4f -shell treatment. Expectedly, with 4f electrons in the core, KCeS 2 represents an indirect-gap semiconductor with large 2.7 and 2.3 eV direct and indirect gaps, respectively, as seen in Fig. 4 (a). Explicit 4f -shell treatment produces a localized state around the Fermi level with a small dispersion driven by d-f hybridization. At the Γ point, the d-f hybridization is minimal with a formal gap of ∼0.05 eV (within the f shell) and approximately 1.0 eV difference between  Next, we estimate the magnetization energy, using the same DFT level as before, by optimizing a single-determinant wave function for different total magnetization densities: (ρ α (r) − ρ β (r))dr = 4 for the FM state and (ρ α (r) − ρ β (r))dr = 0 for the AFM state. Here, a proper consideration would require a supercell at least accommodating one propagation vector. However, the system of this size will suffer a lot more from well-known DFT shortcomings than a smaller system, and the question of applicability alone would need a separate study. Therefore, here we restrict the consideration to a fragment of a single isolated CeS 2 layer in a 2 × 2 supercell, as shown in Fig. 5. We find that thermodynamically the AFM state is more stable than the FM state by 2.0 meV/Ce, as indicated in Figs. 5 (a,c). Similar to the complete system, the f shells form condensed and localized bands around the Fermi level for both spin constraints [see Fig. 5 (b,d)]. Although the f shell in the partial density of states (PDOS) has significant rearrangements, the spin density is well localized on the Ce sites and looks very similar before and after the spin flip. Further investigation of the moment orientation will require noncollinear consideration with a spin-orbit coupling Hamiltonian in place. This consideration would, even more, stretch the reliability of the DFT model in its application to f -shell systems modeling. Moreover, 4f systems have a strictly multiconfigurational character of the ground state |J, m J -multiplet (with H SO > H LS ). This cannot be easily accounted for with any DFT model. Multiconfigurational methods such as the complete active space self-consistent field (CASSCF) are generally inapplicable to a periodic system and limited to cluster approximations.

CASSCF(1,7)/RASSI/SO modeling
Luckily, we find our f shell well isolated in the system. From band structures in Fig. 4, we find that around the Γ point, the 4f states remain factorized from unoccupied states (lack of dor p-projections and small d contributions along the k path). Also, the noticeable contribution of the f shell to the valance bands is likely an artifact of d+p and f -shell projection procedure.
Thus, to understand the local properties, we can restrict our consideration to the local polyhedron (minimal cluster of CeS 6 ) only. Moreover, as it was shown previously, the local properties of Ce centers in ESR spectra can be rationalized from very simple approximations [54]. Here, we elaborate more on this simple theoretical picture using the full strength of 2 F 5/2 multiplet modeling based on Stevens parameters predicted through either point-charge models (using MCPHASE [70] program and in-house Python scripts) or the CASSCF(1,7)/ANO-RCC-VDZP/RASSI/SO modeling (with OPENMOLCAS code [71,72]). We find that apart from some small scaling, both models provide virtually identical results and, for all intents and purposes, can be used interchangeably. We also apply the same level of theory to the impurity phase Ce 2 O 2 S for comparison and better understanding of the experimental results. Figures 6 (b, d) show two different projections of the CeS 6 and CeS 3 O 4 , respectively, based on the experimentally refined lattice parameters. In the case of CeS 3 O 4 with an easy-axis anisotropy, the direction of the moment coincides with the direction of the main magnetic axis along c. For the "easy-plane" situation in the CeS 6 cluster, the anisotropy ellipsoid and the radii are given by the g tensor provided in Table 2. Still, it is possible to evaluate total angular momentum components along the main magnetic With no external perturbation, this moment has Larmor precession in a plane that is tilted 10 • with respect to the ab plane, which is consistent with the moment canting suggested by the refinement of the magnetic structure from neutron diffraction in section 3.  To gain more confidence in the employed theoretical framework, we first studied the impact of the compression/elongation along the C 3 -axis of local Ce polyhedra (CeS 6 for KCeS 2 and CeS 3 O 4 for Ce 2 O 2 S). Figure 6 shows the multiplet transformation (state energies and g tensor of the ground state) as a function of the elongation δ z defined relative to the experimentally refined crystal structure. Interestingly, drastic changes in the state g-energy diagram are observed within a small deviation of only ±0.1 Å from the experimental value. For example, for KCeS 2 (CeS 6 cluster) the energies for the second and third states can change almost twofold and threefold, respectively, in this range. At the same time for Ce 2 O 2 S the situation is much more complex. Apart from the pronounced energy change, the system undergoes two local phase transitions from the easy-axis to easy-plane and back to easy-axis as δ z changes from 0.0 to 0.1. The detailed analysis of this phenomenon is outside the scope of this work, so we leave it for follow-up studies. In the following, we compare relative state energies around the experimental geometries and their compositions ( Table 2).
In our previously published inelastic neutron scattering (INS) data [54], we observed one additional CEF excitation beyond those allowed for the Ce 3+ ion in KCeS 2 . Such additional CEF lines may arise from an impurity phase [40], from modified local environment due to lattice defects, such as stacking faults [48,73], or from vibron quasibound states that result from strong magnetoelastic coupling [74,75]. It is therefore important to verify whether the Ce 2 O 2 S impurity phase identified in our present work can account for this additional excitation. In the framework of the point-charge model, we have optimized the S-ion charge at q S = −1.5e, so that the first excited state of the 2 F 5/2 multiplet matches the most intense experimental CEF transition at 46 meV, while using the experimental geometries for polyhedral CeS 6 and CeS 3 O 4 . In the case of CeS 3 O 4 , the O-ion point charge is assumed equal to that of the S ion, scaled by the ratio of Hirshfeld population (P O /P S ) predicted by the DFT model: q O = q S (P O /P S ) = −1.5e · 5.0/3.7 ≈ −2.0e. The calculated intensities of the INS cross section for CEF transitions in both compounds are calculated using the state energies and the wave functions and are given in Table 2. Taking the linear combination of the resulting spectra in the proportion KCeS 2 : Ce 2 O 2 S = 85 : 15 (corresponding to the results of structural refinement given in Table 1), we can reproduce both the energies and intensities of the experimental INS spectrum with remarkable accuracy (see Fig. 7). This modeling suggests that the highest-energy peak at 74 meV originates from the KCeS 2 majority phase, despite its low intensity, whereas the additional peak at 62 meV can be assigned to the oxysulphide phase. As long as the charge of the ions changes consistently for both phases, the Ce 2 O 2 S excitation is always and Ce 2 O 2 S (P 2 ) phases with experimental local polyhedron geometries (δ z = 0) and an average spectrum weighted with experimental proportion between P 1 and P 2 . located between the two peaks originating from KCeS 2 . This leaves no ambiguity in the assignment of the experimental CEF lines.

The effective Hamiltonian parameterization
To understand the origin of magnetic ordering in KCeS 2 and be able to compare it with rare-earth delafossites with a quantum-disordered ground state, we need to estimate effective interactions between local moments. The chosen DFT level predicts 2.0 meV/Ce, but this method is notoriously ill-equipped for such tasks. Thus, ab initio verifications are needed.
For ab initio estimations, however, one has to consider at least two interacting Ce 3+ ions at the same quantum-mechanical level. Fortunately, a pair of Ce 3+ atoms with 4f 1 configuration each and an AFM ground state (M = 2S + 1 = 1) can be described by CASSCF (2,14)/ANO-RCC-VDZP/RASSI/SO problem with just 105 roots (Slater determinants) [71]. This model produces the six low-lying spin orbitals within 60 meV. They comprise a single ground-state singlet, followed by a quasi-doublet at 5.1 meV, an excited singlet at 10.3 meV, and an excited quasi-doublet with an energy of 53.3 meV. For the single Ce 3+ -ion, the ab-initio ground-state multiplet structure can be reproduced by the minimal |J, m J basis using the derived ab-initio Stevens parameters, B k q , with 99.98% efficiency [72]. Assuming no significant electronic charge rearrangements (Ce 3+ sites are well-localized) and the Lines model as a valid approximation [76], we can project the two-site ab initio solution on the following effective Hamiltonian in the |J 1 , m J 1 ; J 2 , m J 2 basiŝ in the attempt to retrieve J αβ couplings. Using the Hamiltonian in Eq. (1), we fit the free parameters J αβ to ab initio states energies predicted at CASSCF(2,14)/RASSI/SO. We find using the PHI program [77] that the three lowest-energy states can be reproduced exactly with the optimized parameter values given in Table 3. Above 15.0 meV, the deviation in energy exceeds tens of meV, presumably due to contributions of the higher-order (multipolar) interactions for which the anisotropic Lines model, with dipolar interactions only, does not account. Moreover, the extracted Hamiltonian is based on the minimal two-site model which lacks the exact crystal symmetry and may overestimate local anisotropy contributions. Nonetheless, this model can describe the four lowlying states including the ground state. Thus, we conclude that this center symmetry model [ Fig. 8 (a)] is able to predict noticeable anisotropy in interactions between local moments and sufficient to control total magnetic order in cluster models. To prove that, these interactions were further implemented in the localized-moments ordering models using the self-consistent mean-field method (see MCPHASE [70] for more details) for clusters of 2, 3, 4, and 13 sites and for fully periodic systems with proper q-space sampling ( Figure 8). Having allowed for the q-sampling and a maximum unit cell up to 4 atoms in any direction, within the proposed interaction Hamiltonian, we obtain from this model the magnetic unit cell shown in Fig. 8 (f). This is consistent with the experimentally observed stripe-yz order and the ordered moment of 0.33 µ B /Ce, which also agrees with the experimental value.

Discussion and conclusions
In this paper, we have presented the results of magnetic structure refinement for the low-temperature AFM state, which was recently revealed in the effective spin-1/2 triangular-lattice antiferromagnet KCeS 2 . It represents collinear stripe-yz AFM order with spins lying orthogonal to the nearest-neighbor Ce-Ce bonds in the ab plane, possibly with a small (∼10 • ) out-of-plane canting of the magnetic moments that is also expected from theory. A similar stripe order was recently proposed for the closely related compound CeCd 3 As 3 [64]. In addition, we have shown that there are no lattice distortions or symmetry-lowering structural transitions associated with the magnetic ordering to within ∼10 −4 . The thermal expansion remains very small (α c < 6×10 −6 K −1 ) below 120 K, which we confirmed for the c lattice constant using capacitive dilatometry. Our experimental results also indicate that cerium oxysulphide, Ce 2 O 2 S, which was present in our sample as a minority phase, does not order magnetically down to 20 mK and may therefore represent a promising spin-liquid candidate deserving a separate study.
Following the experimental determination of the magnetic structure, we could explain it using firstprinciples calculations that predict exchange-coupling anisotropy in the pair interaction between two neighboring Ce 3+ ions. This resulting coupling tensor can be regarded as ab initio because it was based on a direct fit of low-energy spectra of the exact RASSI/SO problem with no further approximations. We find that only a limited number of spin-orbitals can be used in the fit, which is direct evidence that at higher energy, coupled states are mediated by higher-order couplings beyond dipolar. We also find that, surprisingly, plain-vanilla DFT model at PBE/PAW level with explicit 4f 1 treatment correctly predicts the dominance of AFM phase over FM with the reasonable magnetization density and the difference in energy between AFM and FM configurations of 2 meV/Ce, which is remarkably close to the mean isotropic coupling, (J xx + J yy + J zz )/3 = 2.1 meV/Ce. Furthermore, comparing the ratio of J i j components in our model with the recent group-theoretical study with a pseudospin-1/2 moment model on a similar lattice [57,63], we find that with J yz /J xx ∼ 1 the KCeS 2 system should be deep in the stripe-yz region of the phase diagram, in consistency with our experimental spin structure.