Critical spin liquid versus valence-bond glass in a triangular-lattice organic antiferromagnet

In the quest for materials with unconventional quantum phases, the organic triangular-lattice antiferromagnet κ-(ET)2Cu2(CN)3 has been extensively discussed as a quantum spin liquid (QSL) candidate. The description of its low temperature properties has become, however, a particularly challenging task. Recently, an intriguing quantum critical behaviour was suggested from low-temperature magnetic torque experiments. Here we highlight significant deviations of the experimental observations from a quantum critical scenario by performing a microscopic analysis of all anisotropic contributions, including Dzyaloshinskii–Moriya and multi-spin scalar chiral interactions. Instead, we show that disorder-induced spin defects provide a comprehensive explanation of the low-temperature properties. These spins are attributed to valence bond defects that emerge spontaneously as the QSL enters a valence-bond glass phase at low temperature. This theoretical treatment is applicable to a general class of frustrated magnetic systems and has important implications for the interpretation of magnetic torque, nuclear magnetic resonance, thermal transport and thermodynamic experiments.

with H (1) including the Zeeman operators, H (2) including bilinear spin interactions, etc. We first present calculations of the interactions in the crystallographic coordinate system (H) up to fourth order, and then describe the effects of the local coordinate transformation H →H introduced in the main text, which removes the anisotropic spin interactions at lowest order.
We have previously 1 estimated the nearest neighbour bilinear couplings using a combination of hopping integrals obtained from ORCA at the PBE0/def2-VDZ level and exact diagonalization of small clusters of molecules. Using this approach, we have extended the calculations to estimate also longer-range couplings and higher order ring-exchange terms.
The Zeeman operator can be generally written: in terms of the external field H, local g-tensor G i and spin operator S i at dimer site i. Within each P 2 1 /c unit cell, there are two molecular dimers, which are related by 2 1 screw axes and c-glide planes, shown in Supplementary Figure 1(a). As mentioned in the main text, the g-tensors may differ for the two dimer sublattices, labelled A and B, by symmetry. For this reason, it is useful to divide the g-tensor into uniform G u and staggered G s components, with: In order to estimate G u and G s for κ-Cu, we performed density functional theory calculations on isolated dimers at the PBE0/IGLO-III level using the ORCA package 2,3 . The molecular geometry was taken from the 5 K crystal structure reported in Supplementary Ref. 4. The principle axes (p, q, r) of the local g-tensors for the A and B sublattices are illustrated in Supplementary Figure 1(b). The largest value g r = 2.010 corresponds to the long axis of the molecules, while the second largest g q = 2.008 lies along the axis connecting the two molecules within each dimer. The final value of g p = 2.002 was found for the third principal axis. In the crystallographic (a, b, c * ) coordinates, the uniform and staggered tensors are estimated as: and: At second order in the spin operators, the bilinear interactions can be generally written: where J ij describes the Heisenberg coupling, D ij is the Dzyaloshinskii-Moriya vector, and Γ ij is a traceless symmetric tensor describing the pseudo-dipolar interaction. We label the unique interactions according to Supplementary Figure 2(a). For example, the anisotropic triangular lattice of nearest neighbour bonds is composed of J and J interactions, while longer range second neighbour couplings are labelled J and J .
(g p , g q , g r ) Within the P 2 1 /c space group, the presence of a crystallographic inversion center forbids a DM interaction between sites i, j belonging to the same sublattice (i.e. D ij = 0). Between different sublattices, D ij is finite. By symmetry, the signs of the a and c * axes components D a and D c * are staggered with k = (π, π) periodicity with respect to the square lattice bonds, as shown in Supplementary Figure 2 Here, I 2×2 is the 2 × 2 identity matrix, t ij is the spin-diagonal hopping, while off-diagonal hopping λ ij arises as a result of spin-orbit coupling. The above restriction on the sum of DM-vectors around any closed loop is equivalent to restricting the product of hopping matrices around any closed path to be a multiple of the identity matrix: We then consider making site-dependent spin rotations, which transform the operators as: in terms of some arbitrary vector v i . The transformed hopping matrices are then: The question of interest is whether we can define a transformation, defined by a specific set T jk = e −iv j ·σ T jk e iv k ·σ =t jk I 2×2 (14) This process can be repeated indefinitely until the loop is about to be closed. The global transformation is consistent only if the string of transformations is compatible with the last T mi also being proportional to the identity. Since the product is invariant: it holds that:T Therefore, the final hopping matrix is automatically made diagonal by this string of transformations, provided Supplementary Eq. (10) holds. This completes the proof for the existence of a transformation that sets allλ ij → 0. The practical implication is that SOC effects can be completely gauged away already at the level of the underlying hopping Hamiltonian, from which the spin couplings are derived. As a result, all interactions appearing in the transformed spin Hamiltonians {H (n) } for n > 1 must take an isotropic form.
In the case of a staggered (π, π) pattern, (as for the components D a and D c * in κ- For small canting angles, i.e. weak spin-orbit coupling, we can work with the approximations The utility of such a transformation is that it transfers the explicitly anisotropic terms to the Zeeman Hamiltonian, which contains then a uniform and staggered contribution: The effective staggered field H s in the rotated framework is orthogonal to the uniform field H u .
For small canting angles, the two field terms become: where we introduced the matrix: In the main text, these total terms are discussed as the total g-tensors in the rotated framework, which are for small canting angles: Note that due to the structure of R and of the g-tensors Supplementary Eq. (5) This gives rise to odd order contributions in perturbation theory with a dependence on Φ = ∂S A · dl. This quantity is independent of the local spin coordinates, and therefore is invariant under the transformation described above.
The dominant three-spin term, illustrated in Supplementary Figure 2(c), is the so-called scalar spin chirality term, with the exchange term given by (up to order t 3 ): where Φ is proportional to the magnetic flux enclosed by the triangular plaquette ijk .
With the assumption of a homogeneous magnetic field (A = 1 2 r × H) we may use the approximation where n is the out-of-plane unit vector and A triangle the area formed by the triangular plaque- For convenience, it is useful to write the Hamiltonian in analogy with a Zeeman term: in terms of the unitless plaquette parameter: which is defined for each closed triangle plaquette ijk . For the triangle pictured in Supplementary Figure 4, we estimated j Φ ≈ 0.039. As noted in the main text, the operator S i ·(S j ×S k ) is isotropic, but these interactions provide explicit contributions to the magnetic torque through the appearance of (H · n) in the coupling.
Finally, we have also considered higher order 4-spin ring-exchange couplings: The distinct four-site plaquette are labelled according to a horizontal (K h ), vertical (K v ) and diagonal (K d ) interaction (shown in Supplementary Figure 2(d)). In previous works in which the ring-exchange terms have been considered, the approximation has been typically taken that K h = K v = K d and K = J J K 6,7 . However, these relations are not enforced by symmetries. Interestingly, as shown in Supplementary Table 1, we find that such relations do not hold when considering the full electronic structure of the dimers. While the effects of such terms on the ground state are a matter of intense study, the isotropic ring-exchange terms do not explicitly contribute to the torque, and therefore may be the subject of future studies.

SUPPLEMENTARY NOTE 2: GENERALIZED TORQUE EXPRESSIONS
As mentioned in the main text, in general the magnetic torque is the derivative of the energy E = H with respect to a reference angle θ: This expression holds strictly in the T → 0 limit, as entropic contributions are omitted for simplicity. After employing the site-dependent rotations described in Supplementary Note 1, the only terms in the Hamiltonian contributing to the magnetic torque are the uniform and staggered Zeeman terms and chiral 3-spin interactions. Here, we show the derivation of the bulk torque contribution for the uniform Zeeman term explicitly.
In terms of the uniform g-tensorG u and laboratory field H, the uniform Zeeman Hamiltonian is: In general, we assume that iS i = 0 at zero field. In the presence of a finite field, there are several subtleties that arise due to anisotropy inG u . For example, the Zeeman energy is minimized when the spinsS i are parallel to the effective field given by: where G T denotes the transpose of G. As a result, we may write: in terms of a general susceptibility χ u . Since the Hamiltonian governing the response of the spinsS i is otherwise isotropic, and we consider a regime with no spontaneously broken symmetry, the susceptibility χ u is isotropic with respect to the effective field, and therefore depends only on the magnitude |H eff,u |. We therefore assume that the susceptibility scales as a power law in terms of the magnitude of the effective field: Combining these expressions, the uniform Zeeman energy is given by: To compute the torque, we need to take the angular derivative of this expression. The θ-dependence arises from the norm of the effective field. For the derivative of the norm we use the following relation: The torque we then obtain is given by: where h is a unit vector in the direction of the laboratory field H, and H is the magnitude of the laboratory field H = |H|. In the main text, this expression is simplified by separating the H and θ dependencies: Analogous expressions follow for the staggered and chiral contributions to the torque.
It is useful to see that these expressions reproduce the conventional sin 2θ dependence in the case when ζ = 0. To show this, we consider the torque in the a − c * plane, with θ being the angle between H and a within this plane. In this case, the corresponding torque is: Assuming the g-tensor is diagonal in the a − c * coordinates, this gives:

SUPPLEMENTARY NOTE 3: IMPURITY SCALING
In this supplementary note, we present the derivation of the approximate scaling expressions for the impurity contributions to the magnetic torque discussed in the main text.
Generically, the coupling of such "orphan spins" to the external field H is governed by the Zeeman Hamiltonian: whereG I,m is the effective impurity g-tensor.
The response of the randomly coupled orphan spins can be understood with reference to the "strong disorder renormalization group" (SDRG) approach [8][9][10][11][12][13] to studying problems with quenched disorder. In the application to spin systems, the central quantity is the distribution of exchange couplings ρ(J). An initial energy scale Ω is set by the strongest interaction within the network, which couples impurity spinsS I,1 , andS I,2 . The relative degrees of freedom associated with these impurity spinsS I,1 andS I,2 are then integrated out, yielding a new "cluster" with total effective spin S eff = |S I,1 ∓S I,2 |, depending on the sign of J 12 . This process modifies the effective interactions, yielding a new distribution ρ Ω (J) of interactions between clusters that is dependent on the energy scale Ω. As Ω is successively lowered, the effective interactions between remaining spin clusters tend towards a fixed point power law distribution ρ Ω (J) ∼ J d/z−1 , which gives rise to power law behaviour in relevant physical observables. Here, d is the effective dimension, and z is the dynamical critical exponent.
In order to derive approximate expressions for response in this scaling regime, it is useful to recast the summation over impurity spins as a summation over the clusters C.
where N C (Ω) denotes the total number of such clusters and n C denotes the number of spins within the cluster C. As impurity spins become successively coupled at lower energies, N C decreases as Ω is lowered. At each energy scale Ω, it is assumed that the distribution ρ Ω (J) is sufficiently broad, that each impurity cluster is approximated as an independent "spin", with effective moment size given by S C,eff . As a result, each cluster is described by a thermodynamic partition function: The contribution of each cluster to the torque is evaluated by taking the angular derivative of the free energy G C = −k B T ln Z C , such that: This yields: In order to simplify this expression further, we introduce the cluster average moment: and approximate the cluster sum by: Similarly, we identify the energy scale with the largest of either the thermal energy or typical Zeeman energy of a cluster: These expressions lead to the impurity torque expression given in the main text by Eq. (14).
In practice, for the purpose of plotting, we use a soft maximum approximation, max(A, B) ≈ (A p + B p ) (1/p) with p = 2. When plotted over several orders of A/B, the resulting functions are largely insensitive to the choice of p.
In order to evaluate the torque expression, the specific scaling of N C and S avg eff with Ω is required. This depends on the nature of the disordered fixed point 10,11 . At any given energy scale, the number of independent clusters scales as N C (Ω) ∼ Ω 0 ρ Ω (J)dJ ∼ Ω d/z , which is an increasing function of Ω. For purely antiferromagnetic and unfrustrated interactions, the pairs of spins integrated out at any given energy scale would always form S = 0 singlets. As a result, no clusters of large moment would be formed as the energy is lowered, and S avg eff would remain fixed. In contrast, in the presence of ferromagnetic or frustrated interactions 10,13 , the average cluster moment must grow as Ω is successively lowered, scaling as S avg eff (Ω) ∼ Ω −κ for some exponent κ. For purely ferromagnetic interactions, the average cluster moment would be directly proportional to the cluster size S avg eff ∝ N −1 C . For interactions with mixed signs, or frustrated antiferromagnetic couplings, the cluster moments increase more slowly Therefore, for κ-Cu and other frustrated systems, this suggests d/z = 2κ is applicable: for some constants S 0 and N 0 .
In the low temperature or high-field limit k B T µ B S avg eff |G T I · H|, solving Supplementary Eq. (53) and (55) leads to: In this limit, all cluster moments should be saturated, leading to the relation for the mag- . Thus, the relation between the exponents is given by: Including prefactors, the reduced torque susceptibility follows from the definition of the torque as τ = µ B d(H · m I )/dθ, together with the above expressions: In the high temperature limit k B T µ B S avg eff |G T I · H| the energy scale is according to Supplementary Eq. (53): Of particular note, this expression (for suitable constants) reproduces the lowest two orders in the expressions for the susceptibiltiy derived in Supplementary Ref. 14 for random 1D chains.
Regarding the NMR linewidth, as discussed in Supplementary Ref. 15-17, the orphan spin impurities contribute through the staggered moment induced in the surrounding bulk around each impurity. As a result, a nuclear spin at site i in the bulk experiences a different effective local field, which is given by H i = H +S i ·Ã i , in terms of the local hyperfine coupling tensorÃ. We assume that the impurity-induced local magnetization is given by S i = a i S I,m , where m labels the impurity closest to the dimer site i. The constants a i are determined, for example, by the distance between i and m. Thus, finite impurity moments will lead to a distribution of local fields, which then broadens the NMR lines according to the specific distribution of a i and S I,m . An important observation is that the magnitude of this broadening depends explicitly only on local quantities, rather than the cluster averages appearing in the total impurity torque. This leads to a different scaling of the NMR linewidth ν with field and temperature. However, within a given cluster, the impurity moments are assumed to remain perfectly correlated. As a result, the external field is able to orient the local impurity spins only through coupling to the total moments of the clusters. These observations can be summarized by the approximation: where the NMR linewidth scales as the root-mean-square impurity magnetization. Here, N is the total number of impurities. As before, we recast the summation in terms of clusters C: where N C gives the total number of clusters, and n C gives the number of original impurity spins within cluster C. In analogy with Supplementary Eq. (49), the contribution per cluster is approximated by: Note that n C appears as a prefactor here instead of S C,eff due to the fact that S I,m · S I,m > 0. We then introduce that average cluster size as: such that n avg N C = N . Finally, making the approximation: provides to the proposed expression: where S avg eff appears only in the argument of the Brillouin function. Finally, it is important to consider experimentally relevant values for the nonuniversal exponents appearing in the scaling forms. As noted above, the torque response of the orphan spin defects is parameterized by the nonuniversal exponent ζ I , which is related to the low-energy distribution of effective interactions, ρ Ω (J) ∼ J 2/ζ I −3 . Although ζ I is unknown a priori, practical considerations restrict 1 z/d ≤ ∞, which corresponds to a narrow range 2/3 ζ I ≤ 1. In general, ζ I is likely to be sample dependent, and should tend to decrease with increasing frustration of the bulk interactions, and/or uniformity of the impurity distribution within the sample. Both such qualities lead to less singular fixed point distributions ρ Ω (J). For example, the limit ζ → 1 (i.e. z/d → ∞) corresponds to the infinite randomness limit, in which ρ Ω (J) is maximally singular. At any given energy scale, the vast majority of spins remain essentially decoupled, leading to a Curie-like response χ ∼ 1/T , up to logarithmic corrections. Such a fixed point describes, for example, the random singlet phase (RSP) 8,9 for purely antiferromagnetic but random interactions in d = 1. In this case, the interactions are not frustrated, and strongly interacting pairs of spins always form singlets, such that large spin clusters do not form at low energies, suggesting κ = 0. In contrast, the opposite limit of a flat energy distribution, given by ζ → 2/3 (i.e. z/d → 1), corresponds to a so-called spin-glass fixed point (SGFP) 10 . This fixed point has been found in d = 2 via the SDRG appraoch for both geometrically frustrated lattices with random but purely antiferromagnetic interactions, as well as bipartite lattices with mixed ferro-and antiferromagnetic couplings 10 . In both cases clusters with large S eff are generated at lower energies. Lying between these extremes are the so-called "large-spin" fixed points (LSFP), which have 0 < κ < 1/2. For example, the inclusion of both ferro-and antiferromagnetic couplings in d = 1 leads to a LSFP 12-14 with κ = 0.21, ζ = 0.83. Similarly, a LSFP also describes randomly site-diluted models in d = 2 with purely antiferromagnetic interactions, yielding κ ∼ 0.1−0.2, depending on the degree of dilution. This corresponds to ζ ∼ 0.8−0.9.
In principle, the impurities in κ-Cu should correspond to a random d = 2 lattice with both site dilution and random ferro-/antiferromagnetic couplings. To the best of our knowledge, the appropriate exponents have not yet been studied for this case. However, it should be emphasized that a relatively large variation in z/d leads to narrow range of susceptibility exponents 2/3 ζ ≤ 1. On this basis, we conclude that the experimental values of ζ exp = 0.76 − 0.83 observed for κ-Cu fall well within the range expected for impurity effects. The observed variance ∆ζ = 0.07 corresponds to about 20% of the realistic range, which may be an indication of strong sample dependence.