Interleaved bond and magnetic frustration in triangular lattice LnCd₃P₃

Geometrical frustration on triangular lattice networks can lead to number of unconventional states arising from the competition between electronic degrees of freedom. Magnetic interactions, for example, across the triangular lattice network are predicted to stabilize states such as the quantum spin liquids or native quantum disordered ground states [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. This frustration is, however, not limited to spin interactions, and the charge degree of freedom can also be impacted via forms of orbital or bond frustration [13, 14, 15, 16, 17, 18].

Interfacing the charge degree of freedom with highly frustrated magnetic states is an experimental challenge in the field of quantum materials. Creating such an interface is key to testing many of predictions of unconventional superconductivity, quantum criticality, and other phenomena that are predicted to arise when itinerant charge carriers interact within a fluctuating background of entangled spins [19, 20, 21, 22, 23, 24, 25, 26]. However, many of the leading materials candidates hosting a highly frustrated magnetic lattice are built from highly localized lanthanide (Ln) ions whose charge degree of freedom is strongly gapped out and are not readily dopable [1]. Conversely, many dopable magnetic materials are often built from transition metal ions that possess extended exchange interactions or other interactions that lower crystallographic symmetries and lift magnetic frustration [27, 28, 29]. Bridging this divide and establishing a cleanly tunable material platform that allows for conduction electrons to coexist within a frustrated, local moment lattice network remains an important goal.

There has been recent progress in stabilizing quantum disordered or quantum spin liquid states across triangular lattice networks of $S_{\mathrm{eff}}=1/2$ magnetic moments derived from lanthanide (Ln) ions. A number of quantum disordered or nearly disordered states have been reported in layered ALnX₂ compounds for instance [3, 30, 31, 32, 33, 34, 35, 36, 37]; however doping itinerant charge carriers into these large gap insulators is a challenge. In the pursuit of an alternative platform that can host a similar triangular lattice network of Ln moments with a smaller, dopable semiconductor charge gap, the Ln(Cd,Zn)₃(P,As)₃ family of compounds presents a number of interesting candidates [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56].

Ln(Cd,Zn)₃(P,As)₃ compounds, based on the prototypical ScAl₃C₃ structure type, host an ideal triangular lattice of rare-earth moments whose anisotropies and ground states can be tuned via Ln composition. There is also evidence that they are tunable semiconductors that can be driven from insulators into metals via steric effects [40] or via slight self-doping effects during crystal growth [42]. In contrast to the delafossite-variant ALnX₂ compounds, strong covalency between the transition metal and pnictide ions leads to a highly disperse valence band, likely enabling p-type doping [56, 39, 47]. An important open question is whether the carriers in these compounds, derived mainly from the transition metal-pnictide interstitial layers, experience comparable frustration effects and whether additional charge or bond order instabilities can be interfaced with the frustrated Ln network.

In this paper, we demonstrate that the LnCd₃P₃ family of compounds hosts a highly frustrated lattice instability via the formation of local bond order across a honeycomb network rooted in distorted trigonal-planar CdP₃ layers within the unit cell. A combination of single-crystal diffuse X-ray scattering and powder pair distribution function analyses identify a local shortening of 2/3 of the trigonal Cd-P bonds, resulting in the formation of one-dimensional CdP chains which are prevented from forming long-range order via geometric frustration in a manner analogous to kagome-ice correlations [57]. Interactions between these distortions can be mapped to a frustrated tiling of dimers across a honeycomb network [58, 59], leading to fluctuating orthorhombic domains which locally break the in-plane rotational and translational symmetry.

Our results establish the LnCd₃P₃ family of triangular lattice antiferromagnets as a rare material platform where frustration of the charge degree of freedom can be interfaced with frustrated magnetism within the same lattice framework, and we further suggest that the unidentified order resolved in previous studies of metallic LnCd₃P₃ [42] and LnCd₃As₃ crystals [49] is a partial freezing of this bond order.

A representative unit cell for LnCd₃P₃ is shown in Figure 1(a). The structure is composed of alternating layers of LnP₆ octahedra, CdP₄ tetrahedra, and CdP₃ trigonal planar units [Fig. 1(b)]. The latter form a honeycomb pattern, with Cd and P alternating at the vertices, reminiscent of hexagonal BN and of the Zn(As/P) layer in hexagonal KZn(As/P)[64, 65] [Figure 1(c)]. Here we denote the two separate Cd crystallographic sites as Cd_trig and Cd_tet; the P sites are denoted in a similar fashion.

In first examining the structure across the LnCd₃P₃ compound series, refinement of average structural models using powder X-ray diffraction measurements reveal large anisotropic displacement parameters for the Cd_trig and P_trig sites in all four Ln = La, Ce, Pr, and Nd compounds. Data and refinements are shown in Figure 2. The values for the structural refinements are reported in the supplementary information (Table S1 [66]) and are consistent with a previous report of PrCd₃P₃ [38]. In all four compounds, Cd_trig cations exhibit large displacements along the c-axis at room temperature, displacing toward the large voids above and below the site. P_trig anions, on the other hand, have larger in-plane displacements at this temperature. Notably, upon cooling toward T = 5 K, both the Cd_trig cation and P_trig anion displacements evolve to be strongly in-plane, revealing a tendency toward two-dimensional disorder within the $ab$-plane.

To explore this further, X-ray PDF measurements of the local structure of polycrystalline samples of LnCd₃P₃ were performed at T = 80 K [Fig. 3(a), Fig. S1 [66]] and 300 K (Fig. S2 [66]). Small-box modeling and fitting to the data reveal a preference for a distorted, orthorhombic Cmcm cell over the average hexagonal P6₃/mmc cell that is enhanced as the fitting range becomes more local [Fig. 3(c)]. The primary driver of the proposed distortion and resulting symmetry lowering is the breaking of rotational symmetry in the trigonal planar CdP₃ units within the CdP layer. This suggests a local contraction of two of the three equivalent Cd-P bond lengths within the Cd_trig layers and a potential local bond order that is frustrated within the plane [Fig. 3(d,e)].

To further study the presence of local symmetry lowering, single crystal synchrotron X-ray diffraction measurements were performed on crystals of LaCd₃P₃, CeCd₃P₃, PrCd₃P₃, and NdCd₃P₃. Because the data are qualitatively similar in reciprocal-space structure amongst all four compounds, a representative $(H,K)$-plane diffraction pattern for PrCd₃P₃ is presented in Figure 4(a). Corresponding data for all four compounds are shown in Fig. S3 [66]. Aside from sharp Bragg peaks at integer values of H and K, extensive diffuse scattering appears connecting Bragg reflections. The scattering is quasi-two dimensional and completely diffuse along L, generating structured planes of scattering perpendicular to $(1,1,0)$ and equivalent reciprocal lattice directions.

Because the diffuse scattering in PrCd₃P₃ shows the highest relative intensities, the stronger signal in PrCd₃P₃ was chosen for analysis as representative of the behavior in all four compounds. Line-cuts of the diffuse scattering are presented along directions perpendicular [Fig. 4(b)] and parallel [Fig. 4(c)] to the planes of diffuse scattering for PrCd₃P₃. Additional line cuts for the other LnCd₃P₃ are presented in Fig. S4 [66]. The line-cuts within the plane show that the diffuse plane is structured and peaked at half-integer positions, indicating a local breaking of translational symmetry.

The peak shapes were fit using a pseudo-Voigt profile, and the in-plane correlation lengths, extracted via the inverse of the full width at half maximum (FWHM) of the Lorentzian component of the profile, are shown in Figure 4(d,e). Correlation lengths as well as the intensity of the local correlations do not change significantly with temperature and are seemingly decoupled from thermally-induced atomic displacements, suggesting that the diffuse intensity arises from frustrated charge correlations within the trigonal planar CdP layer.

To further resolve the origin of the short-range charge correlations, 3D-$\Delta$PDF analysis of the X-ray scattering data [67, 68] was performed using a modified approach to the punch-and-fill method [69]. A representative $(H,K)$-plane of the isolated diffuse scattering with the average structure removed and the real part of the Fourier transform are shown for PrCd₃P₃ at T = 300 K in Figure 5(b,c).

The vector probabilities of atomic displacements in the experimentally-derived 3D-$\Delta$PDF are highest in intensity at integer lattice units in real space, corresponding to distances between like atoms in different unit cells. A quadrupole-like displacement feature appears at the first-nearest neighbor $\langle 1,0,0\rangle$ positions and another in-phase displacement correlation appears at the $\langle 2,0,0\rangle$ positions, indicative of strongly one-dimensional correlations along this direction. These correlations decay exponentially and have correlation lengths on the order of 5 unit cells, as demonstrated in Fig. S5 [66].

While the three Cd-P bond lengths in the CdP₃ layers are equal in the hexagonal cell, powder PDF results indicate one bond is shorter than the other two in the Cmcm cell. For instance, in the T = 80 K PrCd₃P₃ data, the shorter bonds are 2.42(11) $\mathrm{\AA}$, whereas the longer bond is 2.53(2) $\mathrm{\AA}$, representing an approximate 2% decrease/increase in the bond length with respect to the hexagonal cell. This distortion effectively reduces the CdP₃ trigonal planar units to one-dimensional CdP chains and maps the local distortion to a two-dimensional tiling of the honeycomb Cd-P network with “dimers”, where each dimer corresponds to the single long bond within each unit. This type of bond order is expected to be frustrated across the honeycomb network and is a likely candidate for the local correlations isolated in the 3D-$\Delta$PDF analysis.

To assess the microscopic nature of these correlations, we employed forward Monte Carlo modeling to generate a large-box model with short-range correlations between sites on the CdP honeycomb lattice, the expected diffuse scattering, and resulting $\Delta$PDF maps. The dimer covering of the honeycomb network is canonically defined such that each vertex of the honeycomb is reached by exactly one dimer. This problem can be mapped to a 2D Ising model on the (dual) triangular lattice [58, 59]. In this mapping, Ising spins decorate the center of each hexagon, and dimers are placed on those edges which sit between frustrated pairs of spins (up-up or down-down).

Within the dual model of Ising spins on an antiferromagnetic triangular lattice with couplings up to third nearest neighbor ($J_{1},J_{2},J_{3}<0$), in the limit of large nearest neighbor $J_{1}$, frustration effects manifest one-dimensional “stripe” textures for $J_{3}/J_{2}<0.5$ and “zig-zag” textures for $J_{3}/J_{2}>0.5$ [70]. The stripe pattern in the Ising model maps to the “staggered” configuration in the dimer model, and the zig-zag pattern in the Ising model maps to a “Herringbone”-like configuration of dimers on the honeycomb lattice (Fig. S6) [66, 71, 72]. The forward Monte Carlo model was built around the 2D Ising problem, and, after relaxing the spin configuration, the coordinates of each Cd_trig and P_trig site were modified to elongate one of the three Cd-P bonds for each CdP₃ unit. The magnitude of these displacements were parameterized based on the orthorhombic Cmcm solution from the aforementioned small box modeling of the powder PDF data, though the model is not constrained by this symmetry.

By fixing a relatively large $J_{1}$ and sweeping the ratio of $J_{3}/J_{2}$ in the forward Monte Carlo simulation, we are able to tune the short-range correlations in a 100$\times$100 supercell between the staggered and Herringbone configurations of the dimer model (Fig. S7 [66]). Examination of the calculated diffuse scattering for each case reveals that the diffuse scattering planes observed in the experiment are closely replicated for $J_{3}/J_{2}\geq 0.5$, and the experimentally observed half-integer maxima are recreated for $J_{3}/J_{2}>0.5$ due to the nature of the symmetry lowering in the Herringbone configuration of dimers.

Optimized Monte Carlo scattering results are presented in Figure 5, where $J_{3}/J_{2}=0.58$. The resulting dimer model configuration shown in Figure 5(a) features an arrangement of one-dimensional CdP chains which are correlated over short length scales ($\approx$ 3 to 5 unit cells), with a weak tendency toward breaking translational symmetry by alternating chain direction. Due to the three-fold degeneracy of distortion pathways, three domains form, each rotated by 120^∘ relative to the other two. The calculated diffuse scattering that results from this short-range ordered supercell is shown in Figure 5(b). The planes of diffuse scattering, including the modulations in the structure factor and the weak half-integer maxima, are well-captured by the model. For further comparison, the 2D-$\Delta$PDF transformation of this simulated scattering also agrees well with the corresponding 3D-$\Delta$PDF transformation of the experimental scattering data [Fig. 5(c)], replicating the key features at integer coordinates in the experimental 3D-$\Delta$PDF.

In order to further strengthen the interpretation presented in our model, reverse Monte Carlo modeling (RMC) was used to refine a large-box structural model against the experimental diffuse scattering. Because the scattering is completely diffuse along $L$, data from the $(H,K,\frac{5}{2})$ scattering plane was used as the target for the RMC model. An 75$\times$75$\times$1 supercell including only a trigonal Cd-P plane was initialized with a random distribution of Cd and P distortions along directions A, B, and C [as defined in Fig. 1(d)], and each RMC move involved a swap of the displacement type of two randomly selected sites of the same element. The termination criteria was defined as 10 times the number of atoms in the supercell. The refined supercell, diffuse scattering, and $\Delta$PDF after convergence are shown in Figure 5(d-f). The RMC generally agrees with the results of the FMC model, replicating a similar chain correlation length and exhibiting a weak tendency to form chains which alternating directions. Taken together, the FMC and RMC results confirm that the origin of the diffuse scattering is bond-dependent, short-range interactions between Cd-P displacive distortions which follow the energy landscape enforced by frustrated bond order and result in a lower symmetry at the local scale.

In $Ln$Cd₃P₃ compounds, the valence band maximum is composed primarily of phosphorus p-states, specifically from those P anions within the CdP layer [39]. The conduction band maximum is composed of cadmium s-states. Given that the Cd_trig cations are bonded solely to the three P anions within the CdP layer, the bonding interactions are therefore likely between Cd_trig s-states and phosphorus p_x and p_y atomic orbitals. These atomic orbitals likely hybridize to enable a trigonal planar coordination environment.

In the ionic limit, any degree of charge transfer from the P^3- anions, which possess filled p_x and p_y bands (valence band maximum), to the trigonal planar Cd²⁺ cations, which possess empty Cd²⁺ states (conduction band minimum), will result in a second-order Jahn-Teller instability in the phosphorus p-states. Similarly, in the covalent limit, strong hybridization between the Cd s- and P p-states would result in sp² hybridized orbitals, which would contain two holes and be subject to a first-order Jahn-Teller instability. In either scenario, a bonding instability is expected to arise from the unconventional trigonal planar CdP₃ bonding environment.

The resulting bond instability generates local CdP chains within the honeycomb network whose configurations are governed via chemical bonding considerations – avoiding over- and under-bonded P ions. The effective interactions of this local bond order on a honeycomb lattice frustrate the formation of long-range order and map to the dual triangular lattice with Ising-like interactions [58, 59]. An alternative way to envision things is to consider local 120^∘ displacements of Cd ions toward P ions in the honeycomb network. These displacement modes have a local Ising-like degree of freedom (toward or away) and can be mapped onto the midpoints of the CdP honeycomb bond network. This creates an emergent kagome network that maps to the kagome-ice problem of antiferromagnetic Ising spins on a kagome lattice, and whose solution in the classical limit can be attained via dimer tiling on a honeycomb network [57].

Within the dual triangular lattice Ising model, extended interactions within a $J_{1}-J_{2}-J_{3}$ model can generate a range of local patterns of order [70]. The “zig-zag” state accessed in the limit of appreciable antiferromagnetic $J_{3}$ in the triangular lattice corresponds to the Herringbone local order of the dimer model revealed in reverse Monte Carlo analysis of the diffuse scattering in LnCd₃P₃. The resulting steric frustration is therefore similar to the frustrated charge order reported in ScV₆Sn₆ [73, 74] and related kagome compounds [75, 76, 77] with an out-of-plane lattice chain instability. An unusual off-centering of the Cd_trig and As_trig ions was reported for similar LnCd₃As₃ compounds, suggesting a generic in-plane instability for this entire class of materials and warrants further structural investigations [56].

The extended interactions mediated via a strain field imply a potentially rich set of states accessible via external strain in the LnCd₃P₃ family of compounds. Pronounced strain effects and coupling of lattice fluctuations are known to couple to the charge order state and electronic properties of kagome materials with a similarly frustrated charge ordered state [78, 79, 80, 81]. Given the narrow band gap of the LnCd₃P₃ family, these compounds are likely dopable into a metallic state as well. Previous studies of LnCd₃P₃ crystals grown via self-flux techniques indeed reported a weakly metallic state with a small carrier density [42]. Crystals grown via this method are likely self-doped (hole-doped) via excess phosphorus within the CdP blocks, but, curiously, the onset of metallic behavior coincides with the emergence of an unknown phase transition at high-temperature [42]. We propose this weak phase transition is a partial freezing of the frustrated bond-order native to the lattice, one whose frustration is relieved via the introduction of local charge carriers.

More generally, the ability to dope the LnCd₃P₃ compounds opens a fascinating new frontier of exploring carrier tuning in a lattice hosting both strong magnetic and steric frustration. A quantum disordered magnetic ground state hosted within the local moments of the $S_{eff}=1/2$ Ln layer can in principle be interfaced with the frustrated bond order within the CdP₃ planes via charge carriers, potentially leading to a number of interesting effects. These could include giant magnetocaloric responses [82, 83] in the magnetic channel or emergent behavior, such as unconventional superconductivity [84, 85] within the charge channel. Future experiments on electronic symmetry breaking in carrier-doped LnCd₃P₃ are therefore highly desired.

In conclusion, we have studied the LnCd₃P₃ family of compounds and have identified a frustrated form of bond order akin to kagome-ice within its honeycomb, trigonal planar CdP layer. The bond frustration derives from a chain instability within the honeycomb network and therefore maps to frustrated Ising spins on the dual triangular lattice. Signatures of this frustration persist to room temperature and exist independently from the magnetic frustration within the triangular lattice Ln $S_{\mathrm{eff}}=1/2$ moments. Our findings motivate future experiments probing the coupling between the two (bond and moment) frustrated networks as well as the impact of carrier doping as a means of stabilizing emergent phase behavior in this fascinating family of compounds.

SDW acknowledges helpful discussions with Leon Balents and Jacob Ruff. The authors acknowledge various forms of support from Casandra Gomez Alvarado, Guang Wu, Matt Krogstad, Dibyata Rout, and Joe Paddison. This work was supported by the US Department of Energy (DOE), Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Grant No. DE-SC0017752. S.J.G.A. acknowledges the additional financial support from the National Science Foundation Graduate Research Fellowship under Grant No. 1650114. J.R.C. acknowledges support through the NSF MPS-Ascend Postdoctoral Fellowship (DMR-2137580). This research made use of the shared facilities of the NSF Materials Research Science and Engineering Center at UC Santa Barbara (DMR-2308708). Use was made of computational facilities purchased with funds from the National Science Foundation (CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF DMR-2308708) at UC Santa Barbara. G.P., B.R.O., and L.K. acknowledge support from the National Science Foundation (NSF) through Enabling Quantum Leap: Convergent Accelerated Discovery Foundries for Quantum Materials Science, Engineering and Information (Q-AMASE-i): Quantum Foundry at UC Santa Barbara (DMR-1906325). This research used resources of the Advanced Photon Source, a U.S. Department of Energy (DOE) Office of Science user facility operated for the DOE Office of Science by Argonne National Laboratory under Contract No. DE-AC02-06CH11357. Research conducted at the Center for High-Energy X-ray Sciences (CHEXS) is supported by the National Science Foundation (BIO, ENG and MPS Directorates) under award DMR-2342336.