Emergent moments and random singlet physics in a Majorana spin liquid

We exhibit an exactly solvable example of a SU(2) symmetric Majorana spin liquid phase, in which quenched disorder leads to random-singlet phenomenology. More precisely, we argue that a strong-disorder fixed point controls the low temperature susceptibility $\chi(T)$ of an exactly solvable $S=1/2$ model on the decorated honeycomb lattice with quenched bond disorder and/or vacancies, leading to $\chi(T) = {\mathcal C}/T+ {\mathcal D} T^{\alpha(T) - 1}$ where $\alpha(T) \rightarrow 0$ as $T \rightarrow 0$. The first term is a Curie tail that represents the emergent response of vacancy-induced spin textures spread over many unit cells: it is an intrinsic feature of the site-diluted system, rather than an extraneous effect arising from isolated free spins. The second term, common to both vacancy and bond disorder (with different $\alpha(T)$ in the two cases) is the response of a random singlet phase, familiar from random antiferromagnetic spin chains and the analogous regime in phosphorus-doped silicon (Si:P).

Quantum spin liquids [1][2][3] are of interest as topological states of magnets. Specifically, their low-energy physics contains unusual degrees of freedom, such as emergent gauge fields and fractionalised excitations carrying corresponding gauge charge [4]. Although elusive in the ground state of a clean system, these excitations are perhaps the most direct signature of the presence of a spin liquid. As these need not correspond to known elementary particles, spin liquids can even make quasiparticles with desired quantum numbers uniquely available [5]. In particular, these excitations can then exhibit interesting, new and unusual low-energy behaviour of their own, where the bulk of the spin liquid can act as a matrix mediating emergent interactions between the defects.
Besides thermal excitations, disorder turns out to be a particularly revealing probe, as it can lead to the appearance ('nucleation') of gauge-charged defects already at T = 0, and especially for gapless excitations, can rearrange the low-energy spectral weight even when the total amount of disorder is small [6][7][8][9][10][11][12][13][14]. In reverse, probing the low-energy physics of an experimental compound can provide insights into amount and nature of disorder present in a particular material [15][16][17][18][19][20].
Our work weaves together several of these threads. We study an SU(2)-invariant version [21] of Kitaev's model [22], on a decorated honeycomb lattice, subject to random site (dilution) and bond disorder. This fractionalised magnet exhibits Majorana fermion excitations coupled to an emergent Z 2 gauge field.
Both types of disorder realise the physics of a strongdisorder random singlet phase, with a characteristic divergence of the low-T susceptiblity, as captured by the second term in with positive α(T ) vanishing with T . Such a random singlet phase is familiar from the physics of random antiferromagnetic spin chains [23][24][25][26][27][28]  ing a random pattern of valence bonds between moments, and a broad distribution of triplet excitation energies for the valence bonds. A random singlet description [11][12][13] has been argued to provide a reasonably good phenomenological fit [15] to recent experiments on several interesting S = 1/2 magnets with geometric frustration and quenched disorder [16][17][18][19][20], provided one includes the effects of spin-orbit coupling whenever present. Here, we employ the exact solubility of the model combined with arbitrary-precision numerics to pin down the low-T structure in considerable detail. Dilution, on top of this random-singlet physics, generates the Curie tail, C/T (Eq. 1), in the susceptibility. This is predominantly due to zero modes of the Majorana fermions, rather than isolated spins. Its presence can thus be used to distinguish between the two types of disorder, while the random singlet signal -subdominant in the case of dilution -exhibits a non-universal exponent depending on the disorder strength.
The integrable model [21] studied here provides an SU(2) symmetric Majorana spin liquid [21,32,33] in which the Z 2 fluxes are static and gapped. It has S = 1/2 moments S R on sites R of a decorated honeycomb lattice in which each honeycomb site r is replaced by a trian-gle r consisting of three sites R 1/2/3 ( r) (see Fig. 1). The three spins of a given triangle are coupled to each other with a large antiferromagnetic exchange J, the largest energy scale in the Hamiltonian. As a result, the lowenergy physics is controlled by states in which each triangle is in one of two S tot = 1/2 doublet states. These two doublets are distinguished by the "orbital" quantum number τ z r = ±1. Neighbouring triangles are coupled by a multi-spin interaction of strength K that is sensitive to the "orbital wavefunction" of the total spin state of each triangle.
When K J, the low-energy Hamiltonian can be written in terms of spin-half variables S r = σ r /2 (where σ r are Pauli-matrices) representing the total spin of each triangle: Here, the first sum is over all nearest-neighbour links r r λ on the honeycomb lattice, λ = x, y, z denotes the orientation of the nearest neighbour link connecting the A-sublattice site r to the B-sublattice site r of the honeycomb lattice, τ z r is the orbital quantum number introduced above, and τ x,y,z r are Pauli matrices in the orbital Hilbert space at each r.
Nonmagnetic impurities, corresponding to missing spins in the original model on the decorated honeycomb lattice, are represented by missing sites in this low-energy Hamiltonian since a single vacancy on triangle r leads to a nondegenerate (inert) singlet state for this triangle, S r = 0. In other words, a nonmagnetic impurity on triangle r of the original model leads to the corresponding site being deleted in H Y L (along with the three K-bonds connecting it to neighbours on the honeycomb lattice). In addition, bond disorder in the coupling K of H Y L arises as a consequence of quenched disorder in the strength of the multispin interaction, so long as the intra-triangle exchange J remains the largest coupling in the system. In what follows, we will assume this to be true and analyze the effects of bond and site disorder in H Y L .
We begin by noting that the Hamiltonian H Y L admits an exact solution [21,22] in terms of a Majorana representation [34,35]: σ z r = −ic x r c y r , τ z r = −ib x r b y r , and cyclic permutations thereof, where c λ r and b λ r are Majorana (real) fermion operators. In the physical Hilbert space, which is characterized by the local constraint: Defining Z 2 gauge fields u r r ≡ b λ r b λ r on λ-links, the problem reduces in these variables to three flavours of c fermions hopping on the honeycomb lattice while coupled to a common Z 2 gauge field u that has no quantum dynamics: where B = hẑ. In other words, the model reduces to a collection of free fermion problems, one for each static configuration of Z 2 fluxes threading faces of the lattice. In consequence, the temperature-dependent susceptibility for the spin model, including the effects of exchange disorder and vacancies, can be determined from the density of states of an associated free fermion system.
A dictionary relating the free fermion problem to the magnetic one is as follows. Physical properties of the SU(2) symmetric model are controlled by the behaviour of a triplet of Majorana fermions hopping on the honeycomb lattice. The hopping matrix K has matrix elements ±iK r A r B u r A r B with the u corresponding to the ground state flux sector (here r A and r B are the A and B sublattice sites connected by the corresponding link of the honeycomb lattice). We introduce canonical (complex) fermions, defined as f r = ±(c x r + ic y r )/2, with the plus (minus) sign on A (B) sublattice sites. The f -fermion Hamiltonian is then a tight-binding model with hopping matrix K, single-particle eigenenergies , and density of states ρ( ). The magnetic field h acts as a chemical potential for the f -fermions, since S z r = −ic x r c y r = 1/2 − f † r f r . The key conclusion is hence that the magnetic susceptibility of the spin model is equivalent to the compressibility of the fermion system and determined by ρ( ). This conclusion must be augmented with a discussion of the consequences of the constraint D r and of the flux sector selected at low temperature [22,[36][37][38]: we omit details for brevity. In general, chiral symmetry of K ensures that ρ( ) is an even function of , and for the clean system ρ( ) ∝ | |. Disorder affects ρ( ) in two possible ways. Exchange randomness generates additional low-energy fermion states with a continuous distribution of energies and a density ρ reg ( ) that diverges as → 0. Site dilution produces both a continuous contribution ρ reg ( ) (similarly diver-gent for → 0, albeit with a different functional form) and also a finite density of zero modes, so that in this case ρ( ) = ρ 0 δ( ) + ρ reg ( ). It is convenient to measure fermion energy and temperature in units of the disorderaveraged exchange coupling K av and to use the parameterisations Γ ≡ log 10 (K av / ) and Γ T ≡ log 10 (K av /T ). In addition, we express the integrated density of states in terms of Γ via N (Γ ) = − d ρ reg ( ).
An elementary calculation relates the susceptibility to the compressibility of the canonical fermions. For T K av we obtain At low temperature, Γ T is large and N (Γ T ) reflects the form of ρ reg ( ) near = 0. Consequences of this mapping for magnetic properties are summarised in Fig. 3. Zero modes of the fermion system yield a Curie tail in the susceptibility, with coefficient C = ρ 0 /4. That this is a cooperative response and not the result of individual free spins is reflected in the form of the corresponding fermion eigenfuctions, which extend over multiple sites. Disorder-induced low-energy fermion states with density ρ reg ( ) are responsible for a second contribution to χ(T ), also divergent for T → 0 but subleading.
These statements can be backed up in considerable detail by appealing to a mapping of the spin liquid to a random bipartite hopping problem [39][40][41][42][43], from which results can be transposed into the present context. Justification of the framework appeals to an underlying infinitedisorder fixed point as the origin of the random-singlet physics. The T -dependence of the second term in Eq. 4 turns out to be N (Γ T )/T = T α(T )−1 , with α(T ) vanishing slowly but non-universally as T → 0, so that different α(T ) arise for site and bond disorder.
For site dilution, our numerical results are summarized below. We find α(T ) ∼ y(n v ) log(Γ T )/Γ T at not-too-low temperature, which crosses over to α(T ) ∼ 1/Γ 1/3 T below the crossover temperature T c ∼ K av 10 −Γc(nv) . Interestingly, for the site-diluted case, y(n v ) and T c (n v ) both decrease quite rapidly with decreasing concentration n v of vacancies, implying a correspondingly stronger singularity in the random singlet form of the susceptibility for lower values of vacancy density. We note that, as all the key features of the low-energy physics of the random bond-disordered case are already present here, further adding such disorder is not expected to lead to any qualitative changes.
The results summarised above are supported by extensive numerical calculations, as illustrated in Fig. 3. Crucially, these calculations use the methods of Ref. [45] to extend to much lower energies previous studies [6,7], which were limited to systems with vacancies but no exchange disorder. High resolution in the low-energy density of states is essential to identify the behaviour described above for N (Γ) and to distinguish exact zero modes from very low-energy contributions to ρ av ( ). The form of K relevant to low-temperature properties of the spin system with site dilution has a Z 2 flux bound to every vacancy, and differs in this way from the tight-binding models for disordered graphene studied in Ref. [45][46][47][48]. Interestingly, the main features of the low-energy density of states are common to systems with and without Z 2 fluxes.
To obtain insight into the physical origin of the Curie contribution to Eq. (4), which arises in systems with vacancies, it is necessary to probe in detail the nature of zero modes in the fermion description. To this end it is useful to compare behaviour in systems with different numbers of vacancies, and with or without disorder in the hopping amplitudes. A single vacancy necessarily gives rise to a single zero mode, because it creates an imbalance sublattice site numbers for the bipartite hopping problem. In contrast, zero modes at finite vacancy density are a complicated multivacancy effect, spread over many unit cells. As noted in Ref. [45], the density of these vacancy-induced zero modes has contributions from both "fragile" zero modes (which are sensitive to the values of the hopping matrix elements) and generic "disorderrobust" zero modes, which remain pinned to zero energy independent of the disorder in the hopping amplitudes. Fig. 2 visualises such a set of local zero modes, obtained using the methods of Ref. 49. Based on a comparison of the ρ 0 computed here with the corresponding results of Ref. [45] without flux attachment, we conclude that such disorder-robust zero modes dominate over fragile ones. Thus, we expect the Curie constant C to be largely determined by the vacancy density and relatively insensitive to bond disorder.
These results fit snugly into the random singlet phenomenology advocated [15] as a reasonably good description of the low temperature physics of a variety of interesting S = 1/2 magnets with geometric frustration and quenched disorder [16][17][18][19][20], with due allowance for spinorbit coupling as necessary. From this perspective, our work adds two interesting new angles. Firstly, it transfers the exact solubility of Kitaev models to the SU(2) random singlet case, thereby permitting us to make controlled statements about a vast temperature range for unusually large systems, compared to what is commonly possible for disordered quantum magnets. Second, it adds a Curie tail to the random singlet physics, with the tail arising from spatially extended vacancy-induced spin textures, rather than the response of isolated free spins. This in turn extends a similar phenomenon from the classical realm, known there as orphan spins [50,51], whose behaviour in a quantum context has remained a puzzle.
Given the rather comprehensive theoretical understanding we have described, the question of experimental implications naturally follows. First of all, the insensitivity of our results with regard to the choice of flux sector -the phenomenology is explicitly the same with and without the bound fluxes -implies that perturbations the only role of which is to favour a different groundstate flux configuration will likely not affect the basic phenomenology. In addition, this may provide some robustness against endowing the fluxes with dynamics of their own, as happens when one tunes away from the exactly soluble point.
For any such perturbations, one would as is customary for gapless frustrated magnets, expect a crossover to different behaviour at an energy scale set by the strength of the perturbation. We emphasize that two central aspects of our preceding discussion should be robust: firstly, the response above this scale; and secondly, the integrated spectral weight below this scale, comprised of the random singlet density of states and that of the Curie tail.
What happens below this scale is then an interesting many-body problem of its own, which a priori needs to be addressed on a case-by-case basis. For appropriately chosen, or perhaps even naturally occurring, peturbations, this may even yield as yet unexplored variants of cooperative physics of the low-energy emergent degrees of freedom in a disordered, strongly interacting topological magnet.
Finally, on the technical side, we note that our work also suggests interesting questions for future work. Can one construct a strong-disorder RG procedure directly in spin language for this exactly solvable model? How are emergent moments seeded in this description, and what prevents them from being quenched by singlet valence bonds? And finally, can this RG approach be used to perturb away from the solvable model?
Acknowledgements: We thank F. Evers, J. Knolle, O. I. Motrunich, N. Perkins and M. Vojta for engaging discussions, and R. Bhola and S. Biswas for assistance with Fig. 2. KD gratefully acknowledges a fruitful collaboration [49] with R. Bhola, S. Biswas, and M. Islam on the random geometry of disorder-robust zero modes. SS and KD gratefully acknowledge use of departmental computational resources of DTP-TIFR and ICTS-TIFR for all the computational work described here. The computational study described here formed part of the Ph.D thesis work of SS, funded by a fellowship from the TIFR. KD gratefully acknowledges partial support from the Infosys-Chandrasekharan Center for Random Geometry at the TIFR. JTC thanks the Dept. of Theoretical Physics of the TIFR for hospitality during a visit in which this collaboration was initiated, and gratefully acknowledges travel support from the Oxford-India Theoretical Physics Network for this visit; his work is also supported by EPSRC Grant No EP/S020527/1. This work was in part supported by the Deutsche Forschungsgemeinschaft under grants SFB 1143 (project-id 247310070) and the cluster of excellence ct.qmat (EXC 2147, project-id 390858490).