Theory of partial quantum disorder in the stuffed honeycomb Heisenberg antiferromagnet

# Theory of partial quantum disorder in the stuffed honeycomb Heisenberg antiferromagnet

## Abstract

Recent numerical results [Gonzalez et al., arXiv:1804.06720; Shimada et al., J. Phys. Conf. Ser. 969, 012126 (2018)] point to the existence of a partial-disorder ground state for a spin-1/2 antiferromagnet on the stuffed honeycomb lattice, with 2/3 of the local moments ordering in an antiferromagnetic Néel pattern, while the remaining 1/3 of the sites display short-range correlations only, akin to a quantum spin liquid. We derive an effective model for this disordered subsystem, by integrating out fluctuations of the ordered local moments, which yield couplings in a formal expansion, with being the spin amplitude. The result is an effective triangular-lattice XXZ model, with planar ferromagnetic order for large and a stripe-ordered Ising ground state for small , the latter being the result of frustrated Ising interactions. Within the semiclassical analysis, the transition point between the two orders is located at , being very close to the relevant case . Near quantum fluctuations tend to destabilize magnetic order. We conjecture that this applies to , thus explaining the observed partial-disorder state.

## I Introduction

Frustration counteracts the tendency of a spin system to order magnetically at low temperatures, and, if strong enough, can lead to quantum disordered ground states devoid of spontaneous symmetry breaking. Such states are conceptually fascinating, a prime example being gapped quantum spin liquids with intrinsic topological order and fractionalized excitations, evoking continued both theoretical and experimental interest.savary_rop17 (); ropfru () A prototypical frustrated spin system is the triangular-lattice antiferromagnet.anderson73 (); ms01 () Following Anderson’s proposal that the spin- Heisenberg model may host a resonating-valence-bond (RVB) state, recent numerical studies confirm that upon amending the nearest-neighbor model with further frustrating interactions, such as second-neighbor exchange, indeed stabilizes a quantum spin liquid.whiche07 (); zhuwhi15 (); sheng15 (); cherny18 ()

While quantum disorder, non-trivial topology, and fractionalization are typically properties of the whole system, intriguing cases with coexisting topological and non-topological components have been discussed in different settings.flst1 (); flst2 (); semoe15 () In particular, it has been suggested that frustration may lead to partially disordered states, where a disordered subsystem coexist with a magnetically ordered part. While these studies have been mostly concerned with Ising modelsmeka77 (); blan84 (); moto12 () or classical systemsdiep97 () at finite temperature, a recent density-matrix-renormalization-group (DMRG) study by Gonzalez et al.gonz18 () has discovered such phenomenology in a SU(2)-invariant quantum spin system.

Ref. gonz18, considered a spin- Heisenberg antiferromagnet on the triangular lattice with a superstructure (also dubbed stuffed honeycomb lattice), i.e., a honeycomb lattice with exchange coupling , supplemented by spins in the center of each hexagon which are coupled by to the surrounding spins, Fig. 1. For a range of finite small couplings , the central spins were found to be decoupled from the Néel-ordered honeycomb subsystem. In this partial-disorder phase, the central spins were found to be short-range correlated, as opposed to previous models of partial disorder (PD), in which the disordered subsystem had a classical extensive degeneracy. Importantly, PD in the stuffed honeycomb lattice is a quantum effect, as the (semi-)classical ground state at any finite involves ferromagnetic order of the central spins, with the honeycomb spins canting in direction of the ordering axis of the central spins. The existence of PD physics in the stuffed honeycomb lattice is supported by an exact diagonalization (ED) study which shows a phase with vanishing magnetization at weak .shima18 ()

The goal of this paper is to identify the mechanism which causes the PD state in the stuffed honeycomb-lattice Heisenberg model. The fact that PD exists in the regime of small allows us to attack the problem perturbatively: We integrate out the degrees of freedom of the honeycomb subsystem, assuming collinear Néel order of the latter. We work to second order in and apply a systematic expansion such that our study is formally controlled in the semiclassical (large-) limit. We obtain an anisotropic effective spin model for the central spins on an emergent triangular lattice. The interactions among the central spins induced by magnons are given by predominantly ferromagnetic transverse () interactions and frustrating Ising interactions for the -components, with the leading terms arising at different orders in : The leading coupling scales as , while the leading Ising interaction scales as . As a result, the effective model displays ferromagnetic in-plane order at large , while at small the Ising interactions become dominant, favoring an out-of-plane stripe-ordered ground state. With all effective couplings calculated to order , we find the transition between the two competing ground states to be located at , being remarkably close to the case of considered in the numerical studies. Moreover, the transition is masked by a small window of more complicated (most likely incommensurate) order, and magnetization corrections grow near . We therefore argue that fluctuations in this frustrated effective model destabilize magnetic long-range order for , which then leads to a PD phase for the stuffed honeycomb antiferromagnet.

The remainder of the paper is organized as follows: In Section II we introduce the Hamiltonian and review previous numerical results. Section III outlines the derivation of the effective spin model, which is discussed in Section IV. A discussion (Section VI) closes the paper.

## Ii Model

We consider local moments on the stuffed honeycomb lattice, consisting of a honeycomb lattice with central spins placed in each hexagon, as depicted in Fig. 1(b). The situation corresponds to a supermodulation of a triangular lattice.

### ii.1 Heisenberg Hamiltonian

The local moments on the honeycomb lattice (with , sublattices) are coupled antiferromagnetically with a coupling . Each central spin (sublattice ) couples with interaction with its nearest neighbors on the sublattices, such that the Hamiltonian of the full model reads with

 HJ =J∑⟨ij⟩→Si,A⋅→Sj,B HJ′ =J′∑⟨ij⟩s=A,B→Si,s⋅→Sj,C. (1)

At , the model reproduces the Heisenberg antiferromagnet on the isotropic triangular lattice, while for the system reduces to antiferromagnetically coupled spins on the honeycomb lattice, with the central spins being decoupled. For the subsequent analysis, it will be convenient to consider a unit cell containing sites, with the primitive lattice vectors given by .

### ii.2 Previous numerical results

The classical ground state of the Eq. (II.1) on the stuffed honeycomb lattice for is depicted in Fig. 1(b), with the -sublattice spins canting at an angle to the direction transverse to the sublattice. This leads to the familiar triangular-lattice order at and collinear order of the honeycomb spins in the limit . The classical system thus displays ferrimagnetic behaviour at , with the ground state having a non-zero magnetization in the spin direction of the central spins.

The authors of Ref. gonz18, focus on the case of small and find that the classical order persists in linear spin-wave theory (LSWT) even for infinitesimal , with the sublattice magnetization taking more than of its classical value. Their DMRG results for the spin- system, however, show clear signatures of a PD phase at , with vanishing and collinear Néel order on the honeycomb subsystem. Notably, in PD the central spins are found to be ferromagnetically correlated, hinting towards the presence of an effective interaction between these spins. For a canted state is found, with canting angles closely matching the LSWT result. At the central-spin magnetization changes abruptly (as do nearest-neighbor spin correlators). Altogether, this points toward a first-order transition between a PD state and a semiclassical canted state as function of , Fig. 1(c).

Furthermore, recent ED results by Shimada et al.shima18 () for systems up to sites indicate that the spontaneous magnetization of the semiclassical ferrimagnetic phase for vanishes non-continuously below some non-zero . This would be consistent with the existence of PD, however the value of the critical and precise nature of the correlations below appear to be difficult to obtain in ED due to limited accessible system sizes.

## Iii Derivation of an effective model

The PD phase is realized for small , suggesting that the coupling between the honeycomb spins and central spins can be treated perturbatively. In the limit it appears justified to neglect feedback effects of the central spins on the honeycomb spins, such that the natural effective model for the system involves central spins on an emergent triangular lattice, with effective interactions resulting from the coupling to the honeycomb system.

The goal of this section is thus to derive this effective spin model for the central spins by perturbatively integrating out quantum fluctuations around the Néel-ordered ground state of the honeycomb spins, which we treat in a expansion. For a detailed derivation we refer the reader to Appendix A.

We emphasize that the magnons in the honeycomb system are gapless, and integrating out theses gapless degrees of freedom may in general induce singular interactions. We thus apply a staggered field in direction of the Néel order on the honeycomb sublattices, giving rise to a gap in the magnon dispersion, and take the massless limit at the end.

Without loss of generality, we parameterize the semiclassical Néel order of the honeycomb system by and such that the Hamiltonian with a staggered field reads . We consider magnon excitations on top of this ordered ground state by representing the spin algebra in terms of Holstein-Primakoff bosons , admitting a systematic expansion in . Expanding , where contains boson operators. The leading order represents the classical ground-state energy. The bilinear piece corresponds to LSWT and can be written as

 H(2)J,h=JS∑q[f(q)aqb−q+h.c.+3(1+h)(a†qaq+b†qbq)], (2)

where with defined as above. can be diagonalized by means of a Bogoliubov transformation, yielding two magnon modes (having neglected a global energy shift and employing inversion symmetry )

 H(2)J,h=∑qϵ(q)(α†qαq+β†qβq), (3)

where with the dimensionless dispersion . The quartic part of the Hamiltonian corresponds to magnon-magnon interactions,

 H(4)J,h=−J4∑i,δ[a†iaiaibi+δ+aib†i+δbi+δbi+δ +2a†iaib†i+δbi+δ+h.c.]; (4)

we note that cubic boson terms are absent in the present collinear case. The quartic terms give rise to self-energy corrections to the magnon Green’s function. In general, the self-energy corrections can be split into a static Hartree-Fock contribution arising from “balloon” diagrams as shown in Fig. 2(b) and a frequency-dependent contribution represented by “sunset” diagrams, with the latter scaling as .zhitom09 (); fn:dimless () An explicit expression for is given in Appendix A.

The interaction vertices for the coupling of the magnons to the central spins (omitting the label for brevity) can be determined by expanding in , yielding up to third order

 HJ′=J′∑q[(√S2Γaqaq+14√2S∑k,pΓ3aq;k,pa†k+p−qakap+√S2Γbqb†−q+14√2S∑k,pΓ3bq;k,pb†−kb†−pb−(k+p−q))S−−q+h.c. +∑k(Γaaq,ka†k+qak+Γbbq,kb†k+qbk)Szq], (5)

where we have made the scaling for the respective vertices explicit (for definitions of the various vertices we refer the reader to Appendix A). Note that Eq. (III) does not contain a term linear in , which corresponds to the fact that mean fields on the central spins vanish – apparently a prerequisite to obtain a PD phase.

We are now in position to integrate out the bosonic modes, yielding an effective model for the central spins which are arranged on a triangular lattice. At order and next-to-leading order in the effective action takes the form

 Seff=(J′)2J∫dτ∫dτ′ ∑qS+q(τ)jxy(q,τ−τ′)S−−q(τ′) + Szq(τ)jz(q,τ−τ′)Sz−q(τ′) (6)

with dimensionless couplings . The relevant processes contributing to this order are shown in Fig. 2: The coupling as a result of transverse magnon modes is given by the magnon propagator and subleading Hartree-Fock corrections for both the propagator and vertex, while the longitudinal results from magnon bubble diagrams.

Note that the effective model breaks spin rotation symmetry down to a symmetry of rotations in the plane as a result of the Néel ordering of the honeycomb spins in -direction.

The time dependence of the couplings in Eq. (III) is a consequence of retardation effects. In order to obtain a “static” spin Hamiltonian for the central spins, we employ the instantaneous approximation , where is given by the limit. This approximation is justified as long as there is a separation of energy scales between the fast (high-energy) magnon modes at an energy scale and the low-energy central spins at a scale .

We thus find the static transversal coupling to be given by

 jxy(q)=−12ωq∑δ=α,β[|Γδ|2(1−ΣHFqSωq)+R[ΓδΓ3δ∗q]4S], (7)

and the longitudinal coupling

 jz(q)=−12S∑k⎡⎢⎣|Γαβq,k|2+|Γβαq,k|2ω(k)+ω(q+k)⎤⎥⎦, (8)

with the explicit form of the vertices detailed in Appendix A. The couplings given above implicitly depend (through the dispersion and the Bogoliubov factors appearing in the vertices) on the staggered field strength . In order to extrapolate and to identify the most important spin interactions, it is convenient to work in real space. Fourier-transforming thus leads to the effective Hamiltonian

 Heff=(J′)2JN∑n=0∑⟨ij⟩n[jxyij(SxiSxj+SyiSyj)+jzijSziSzj], (9)

where denotes a bond between -th nearest neighbors on the triangular lattice of central spins. Working to order , the transverse coupling contains both a leading and subleading contribution, while the longitudinal coupling scales as .

We truncate the generated (long-ranged) interactions after the -th nearest neighbor, finding that the properties of the effective model, discussed in Sec. IV, do not significantly depend on the truncation after the most dominant interactions are taken into account. The limit is taken by first evaluating the expressions in Eqs. (7) and (8) for fixed , and then fitting the results to a leading-order scaling form in obtained by a continuum approximation. For details on the extrapolation, we refer the reader to Appendix B. The obtained values of the couplings are shown in Table 1.

We note that the above scheme naturally yields on-site spin couplings (), corresponding to (retarded) self-interactions of the spins mediated by magnons. In the instantaneous approximation, the on-site couplings combine to a single-ion anisotropy for any , while simply yielding a global energy shift at . In order to allow for a consistent large- study, we henceforth include the single-ion anisotropy in the analysis of the effective model (i.e. in classical iterative minimization and linear spin-wave theory).fn:nosingleion ()

## Iv Analysis of the effective model

The anisotropic effective model obtained in the previous section contains an interaction which has a dominant ferromagnetic coupling on nearest neighbors to order . As visible from Table 1, this coupling gets significantly reduced at order by self-energy and vertex corrections. The most important -Ising coupling is an antiferromagnetic second-neighbor interaction.

### iv.1 Qualitative discussion of dominant processes

We quickly discuss the physical picture behind the most important effective interactions. The nearest-neighbor (transverse) interaction is pictorially shown in Fig. 4(a). Two nearest-neighbor central spins can exchange their spin orientation with an intermediate honeycomb-spin flip (which, loosely speaking, can be understood as a one-magnon excitation). As the energy cost of the intermediate state is , second-order perturbation theory leads to matrix elements of the form

 ⟨↑↓|Heff|↓↑⟩∼−(J′)2J, (10)

corresponding precisely to the scaling of the effective action in Eq. (III). Note that this exchange process constitutes the leading-order contribution to . It can only be realized for antiparallel nearest-neighbor spins in the initial and final configurations, thus leading to an effective ferromagnetic interaction, as also observed in the effective model. We note that a dominant ferromagnetic nearest-neighbor coupling in the effective model is consistent with the observed ferromagnetic spin correlation for PD in the DMRG data.gonz18 ()

The dominant contribution of the (longitudinal) -Ising interaction is antiferromagnetic on next-nearest neighbors. The computation of this coupling (cf. Sec. III) involves the evaluation of a magnon-bubble diagram, as displayed in Fig. 2, corresponding to the coherent two-magnon excitation. This process can be illustrated by considering flipping two neighboring spins on the Néel-ordered honeycomb lattice (we emphasize that the spin-flip is not equivalent to a magnon excitation, which is rather a coherent superposition of the spin flips due to the Bogoliubov transformation), as shown in Fig. 4(b). This double spin flip induces antiparallel fields along the spin quantization axis on next-nearest hexagons, therefore favoring an antiparallel alignment of the central spins in the respective two hexagons. This process hence induces an antiferromagnetic Ising interaction of next-nearest neighbors on the effective triangular lattice. It appears likely that this interaction is further enhanced when taking into account magnon-magnon interactions, because of nearest-neighbor magnon attraction. However, this only contributes to at order and is not included in the present calculation.

### iv.2 Classical ground states

Since the strengths of and interactions display a different scaling as function of , we expect the nature of the ground state to change as function of . We hence inspect the family of effective models parameterized by : For each fixed we have a specific triangular-lattice XXZ model, which we again analyze in a formal expansion.

In this subsection, we start by considering the classical ground states. Noting that at only the interaction in (9) is of relevance (and subleading corrections are small), while at sufficiently small the terms are dominant, it is legitimate to consider the two interactions separately. First, we determine the order favoured by the interaction by finding the wavevector corresponding to minima of according the Luttinger-Tisza method. One finds a unique minimum at for any , corresponding to ferromagnetic order with spins lying in the -plane (subsequently dubbed “XY-FM”).

The Luttinger-Tisza method cannot straightforwardly applied to the case of pure Ising interactions. Instead, we approach the problem of finding the ground state of The Ising part of the Hamiltonian by iteratively minimizing the energy on finite lattices of various sizes. We find that the system orders in a two-up-two-down stripe phase, commonly dubbed ,fiselke80 () which spontaneously breaks the triangular lattice rotation symmetry, corresponding to the ordering wave vector (and -symmetry related). The phase diagram of the Ising model with competing (i.e. ferro-antiferro on nearest and next-nearest neighbors) interactions on the triangular lattice has been explored extensively by Kaburagi and Kanamori, confirming as the ground state in the parameter range relevant to our model.kaka74 (); kaka78 () It has been noted that in general competing short-range ferromagnetic and longer-range dominantly antiferromagnetic interactions often favour striped phases,lieb11 (); giu16 () and can lead to rich phase diagrams with exotic critical properties.fiselke80 (); mumaku95 (); jinsensan12 ()

Returning to the combined XXZ model for a given in Eq. (9), we study the (-dependent) competition of transverse and longitudinal interaction terms at the classical level. Numerical results on various system sizes show that above , the classical ground state is given by XY-FM, while for smaller the lowest-energy configuration is given by , see also Fig. 5. This appears to be consistent with the qualitative discussion above in which we noted that at large the ferromagnetic interaction is dominant, while the coupling becomes important compared to for sufficiently small .

We note that there is a small window around in which both XY-FM and become unstable in favor of an incommensurate phase (dubbed “IC”). To map out the extent of this window, it is more convenient to use LSWT, as described in the following section.

### iv.3 Magnetization corrections in linear spin-wave theory

The fact that the transition between ferromagnetic and stripe order occurs at , very close to the case of considered in the numerical studies, indicates that the physics of PD is closely linked to the appearance of this transition. It is also clear that the competition of XY-FM and at small (and thus PD) is an inherently quantum effect, driven by spin fluctuations of the honeycomb lattice. To further explore the competition between the and interactions in and the resulting quantum fluctuations, we consider magnon excitations on top of either ferromagnetic or the striped ground state of the -dependent effective model (9) in LSWT.

Expanding about the ferromagnetic state, we find that the spin-wave dispersion becomes imaginary for . The dispersion closes at the (incommensurate) wave vector (and -symmetry related). Complementary, expanding around the classical reference state, the spin-wave dispersion becomes imaginary for at a very small . The finite window in which LSWT for both XY-FM and breaks down hints at an additional phase, with a possible incommensurate ordering as indicated by the gap closing in LSWT in the XY-FM. Iterative minimization in the IC window, e.g. at indeed yields configurations in which the spins align ferromagnetically in the -plane and have a (incommensurately) modulated -component. Inspecting the static structure factor, we find the ordering wave vector to be given by (to our numerical accuracy). The ordering pattern appears to be related to , as also seen from the fact that is close to . We have checked on finite-size lattices that the classical energy of the incommensurate configuration is indeed lower than the competing commensurate reference states. This incommensurate phase thus masks a critical value at which the energies XY-FM and are degenerate, as we demonstrate employing a toy model in Appendix C. The phase diagram resulting from this discussion is in Fig. 6.

In addition, we have computed the magnetization corrections in LSWT for the two respective ground states, these are also shown in Fig. 6. We note a strong increase of the magnetization correction as the spin-wave theory breaks down when tuning towards IC from both sides. Therefore we consider it likely that long-range order disappears completely near once quantum fluctuations are fully taken into account.

Given that the derivation of the effective model is based on a expansion (and involves further approximations, cf. Sec. III), we conjecture that the true is close to the physical value , and the central spins are in a quantum disordered state for . This then yields the PD state observed in numerical studies.gonz18 (); shima18 () Clearly, higher orders in as well as higher orders in , the latter generating multi-spin exchange interactions which themselves tend to destroy long-range order,motru04 () may be important for a fully quantitative understanding.

It is worth emphasizing that the theme of supplementing a frustrated Ising model by (ferromagnetic) transverse interactions, yielding strong quantum fluctuations, is also key for stabilizing a spin liquid in the easy-axis kagome-lattice spin model studied by Balents, Fisher, and Girvin.bfg02 ()

## V Correlated partial disorder beyond the effective model

In this section we comment on features of the putative correlated PD phase beyond our effective model.

### v.1 Energetic competition between PD and (semi-)classical canting

As visible from the phase diagram in Fig. 1(c), the PD phase competes energetically with semi-classical canted state. We thus discuss the energy contributions of the two competing phases as a function of .

In the classical ground state of the model (as discussed in Sec. II), the honeycomb spins cant at an angle with the central spins pointing in a direction transverse to the Néel order. For this yields a canting angle . Considering that the interaction between honeycomb and central spins is of the form and using that for small canting angles, we thus find that the mean-field energy of central spins in the canted phase scales as . Further we have . Taken together, we see that the energy gain of the canted state relative to the decoupled (collinear) state at scales as .

The energetics of the PD phase is determined by the effective model for the central spins, while mean-field energies between the two subsystems vanish. The effective model was derived in second-order perturbation theory, see Eq. (III). It is thus clear that the energy gain of the PD phase, relative to the decoupled state at , also scales as .

With the two phases having the same energetic scaling, no further qualitative arguments can be made as to which phase has a lower energy. We conclude that – if the PD phase has lower energy than the canted phase for small , as indicated by the numericsgonz18 (); shima18 () – the first-order transition at will occur at of order unity, where contributions in higher order in become important.

### v.2 Topological properties

In PD, the central spins are in a correlated quantum disordered phase. Although the precise nature of the phase cannot be obtained on the level of our analysis, it appears likely that this disordered state of the central spin subsystem by itself possesses topological order. Note that this implies that the notion of topological order is expected to be applicable to the entire system, in the sense of the existence of superselection sectors.

We mention analogies to two systems in which a subsystem is in topologically state. First, a topological spin-glass phase has been proposed in diluted spin ice, in which defect-induced “ghost spins” eventually freeze while the remaining system stays in the Coulomb phase.semoe15 () However, in the present PD phase the situation is reversed: The fluctuations of the ordered “bulk” state stabilize a disordered phase of a subsystem. Second, the model for correlated PD discussed here also bears similarities to fractionalized Fermi liquids (FL), in which conventional electronic charge carriers coexist with local moments which itself form a fractionalized spin liquid.flst1 (); flst2 () One conceptual difference is that the two subsystems of FL can be adiabatically decoupled in two-band models, while such a decoupling is not possible in the present PD phase.

## Vi Summary and Outlook

We have provided an effective theory for the correlated partial disorder phase detected in Ref. gonz18, in the stuffed honeycomb antiferromagnet, by deriving an effective model for the central disordered spins in the presence of collinear background order. Our perturbative treatment is controlled in the limit and includes subleading corrections in . Within our semiclassical analysis, the effective triangular-lattice XXY model undergoes a transition from an in-plane ferromagnet to an out-of-plane stripe phase, the latter driven by competing Ising interactions, at a critical , with an intermediate incommensurate phase. With magnetization corrections growing near the transition, we argue that fluctuations due to the competition of ferromagnetic and stripe orders lead to a quantum disordered ground state for the central spins. Given that is close to , we conjecture that this mechanism drives the correlated partial disorder phase observed numerically.gonz18 (); shima18 ()

Our analysis calls for further numerical investigations. First, it would be of particular interest to investigate the obtained effective triangular-lattice XXZ model beyond the semiclassical techniques employed here. A DMRG study could provide further insights into the nature of the quantum disordered phase anticipated for . Second, one can modify the behavior of the full stuffed honeycomb-lattice model, by introducing additional interactions between the central spins. Based on our results, we predict that explicit antiferromagnetic interactions of order either between first-neighbor or second-neighbor sites on the central-spin lattice lead to stripe order on the sublattice.

On the experimental front, the -distortion of triangular lattice leading to a stuffed honeycomb lattice is of significance for modelling basal planes in magnetic perovskite compounds. Notably, shows signs of PD at low temperatures, however due to a large easy-plane anisotropy the local moments on the plane are rather thought to be described by an effective spin-1 XY model.cope97 (); atmmh83 () The stuffed honeycomb lattice Heisenberg antiferromagnet has also been studied as a model for the triangular lattice spin- antiferromagnet LiZnMoO, in which 2/3 of the spins are in a disordered state below a critical temperature.mcq12 () It has been suggested that the honeycomb local moments form a spin liquid which is stabilized by the central spins as magnetic impurities (as opposed to the case discussed in this work, in which the central spins enter a quantum disordered phase). However, in these studies, a frustrating next-nearest neighbor coupling is crucial for destabilizing the magnetic order of the honeycomb lattice.flint13 (); flint18 ()

###### Acknowledgements.
We thank L. Janssen for useful discussions. This research has been supported in part by the DFG via SFB 1143 (project A01).

## Appendix A Derivation of the effective model

We integrate out the fluctuations of the honeycomb spins in the form of magnons by employing a path integral approach. Denoting bosonic magnon modes as and the central spins by , the partition function in imaginary time for the full system is given by

 Z=∫D[a,b,S]e−SJ−SJ′. (11)

Expanding the exponential in , we obtain

 (12)

with the effective action given by

 Seff=⟨SJ′⟩J+12(−⟨S2J′⟩J+12⟨SJ′⟩2J)+O(J′3) (13)

where denotes the expectation value with respect to the honeycomb magnon action, and the corresponding partition function. The actions noted above are to be taken in imaginary time.

### a.1 1/s expansion

We find the magnons of the honeycomb system in a expansions by employing the Holstein-Primakoff (HP) representation with bosonic modes ,

 S+A =√2S−naa,S−A=a†√2S−na (14) SzA =S−na, (15)

and similarily for , after having picked a local basis by rotating the spins on the -sublattice by around the -axis. After expanding in powers in , the Hamiltonian reads , where corresponds to the classical ground state energy and , which is bilinear in the bosons, is given by (2). After performing a gauge transformation with , the Hamiltonian can be diagonalized by a Bogoliubov transformation

 aq =uqαq−vqβ†−q (16) b†−q =−vqαq+uqβ†−q, (17)

where and with , yielding Eq. (3). The free magnon Green’s functions are then given byfn:dimless ()

 Gαq(iω)=Gβq(iω)=J(iω−ϵq)−1. (18)

#### Self-energy correction to magnon propagator

The leading-order correction to the magnon propagator due to quartic interactions of the HP bosons as given in Eq. (III) can be obtained by normal ordering the magnon operators , after a Bogoliubov transform, or equivalently by a static mean-field decoupling scheme. To this end, we note that the non-vanishing boson bilinear expectation values in momentum space are given by

 ⟨a†qak⟩ =⟨b†−qb−k⟩=v2kδq,k (19a) ⟨aqb−k⟩ =⟨a†qb†−k⟩=−ukvkδq,k. (19b)

Decoupling into the above channels and Bogoliubov-transforming, one obtains the quadratic Hamiltonian

 :H(0)J: =J∑qΣHF(q)(α†qαq+β†qβq) (20)

with the Hartree-Fock self-energy given by

 ΣHF(q)=−16 9(1+h)−|f(q)|2ωq ×∑k[9(1+h)−|f(k)|2ωk−3], (21)

which coincides with Oguchi’s result for bipartite antiferromagnets.oguchi60 ()

Having obtained the magnon self-energy, we note that inverting Dyson’s equation for the magnon Green’s function, yields , with the static self energy given by . It should be emphasized that by summing up all 1-particle irreducible contributions to the propagator correction through the use of Dyson’s equation, one automatically takes all powers of into consideration, which leads to an inconsistency in our systematics. To remedy this inconsistency we use the fact that the self-energy is subleading ( while ) and expand the interacting Greens function in powers of , obtaining

 Gq(iω)=Gq(iω)+Gq(iω)ΣHFqGq(iω)+O(1/S3). (22)

#### Magnon-central spin interaction

Expanding in yields the Hamiltonian given in Eq. (III), with the vertices defined as

 Γaq =fA(q)e−iϕq,Γbq=fB(q), (23a) Γaaq,k =−e−iϕkf∗A(q),Γbbq,k=f∗B(q), (23b) Γ3aq;k,p =−fA(q)eiϕk+iϕp−iϕk+p−q, (23c) Γ3aq;k,p =−fB(q). (23d)

The structure factors and for coupling the and sublattices to the central spins are given by and . We now rewrite in terms of the Bogoliubov modes , , yielding

 HJ′J′ =∑q{[(√S2Γαq+14√2SΓ3αq)αq +(√S2Γβq+14√2SΓ3βq)β†−q]S−−q+h.c. +∑k[Γααq,kα†k+qαk+Γααq,kβ−(k+q)β†−k Γαβq,kα†k+qβ†−k+Γβαq,kβ−(k+q)αk]Szq−Sz0Γbb0,q}, (24)

where we have dropped in the interest of brevity all arising three-boson terms after normal ordering in anticipation of the Hartree diagram shown in Fig. 2(c): The boson loops will then just give , occupation numbers, which vanish at . However, the normal ordering does give rise to corrections to the vertices and as shown above. Note that the last summand arises from re-ordering the term for notational convenience and is cancelled after normal ordering . For further considerations this term will be neglected as it does not contribute to any connected diagrams. The vertices in (A.1.2) are given by

 Γαq =Γaquq−Γbqvq,Γβq=−Γaqvq+Γbquq (25a) Γ3αq =(2Γ3aquq−2Γ3bqvq)∑kv2k (25b) Γ3βq =(−2Γ3aqvq+2Γ3bquq)∑ku2k (25c) Γααq,k =Γaaq,kuk+quk+Γbbq,−(k+q)vk+qvk (25d) Γββq,k =Γaaq,kvk+qvk+Γbbq,−(k+q)uk+quk (25e) Γαβq,k =−Γaaq,kuk+qvk−Γbbq,−(k+q)vk+quk (25f) Γβαq,k =−Γaaq,kvk+quk−Γbbq,−(k+q)uk+qvk. (25g)

For the subsequent analysis it is useful to note the (and analogous).

### a.2 Longitudinal coupling

We can now compute the connected diagrams that contribute to Eq. (13) at quadratic order. The two magnon diagrams relevant for interactions are shown in Figs. 2(a) and (b). As explained above, the magnon loop to be computed in the Hartree correction (b) reduces to a renormalized spin-boson vertex at , so that the leading and subleading contributions to the transversal coupling are simply given by a single magnon contraction. With the Green’s functions for the (free) magnons and similarily for one obtains (setting for convenience)

 jxy(τ,q)=12 {Gα(τ,q)[S|Γαq|2+14R(ΓαqΓ3α∗q)] +(α→β,τ→−τ)11}. (26)

Fourier transforming [with as given in Eq. (22)] and taking the static limit then yields Eq. (7).

### a.3 Transversal coupling

Turning to the longitudinal coupling, we find that at , only particle-particle bubbles for the bosons (as shown in Fig. 2) contribute. For the non-vanishing terms, after changing to the frequency domain (with bosonic Matsubara frequencies ), we find

 jz(iν,q)= −12Jβ∑ω∑k{|Γαβq,k|2Gαq+k(−iω−iν)Gβ−k(iω) +|Γβαq,k|2Gαk(iω)Gβ−k−q(−iω+iν)}. (27)

The Matsubara summations can be evaluated with standard methods to yield, at and using inversion symmetry ,

 1Jβ∑ωGk+q(−iω)G−k(iω−iν)=Jiν+ϵk+ϵk+q, (28)

and analogous for the second term in Eq. (A.3). Taking the static limit yields Eq. (8) in the main text.

## Appendix B Gapless limit by h→0 extrapolation

The expressions for the effective couplings given in Eqs. (7) and (8) depend on the staggered field . We evaluate the couplings numerically for fixed on a lattice with unit cells up to and Fourier-transform to real space according to

 jij=1N∑qj(q)ei→q⋅(→ri−→rj). (29)

Since at any finite the dispersion is gapped and correlations decay exponentially, we perform a finite-size fit of the form , with and as free parameters.

We now discuss the extrapolation which involves non-analyticities. As the gap of closes for , the momentum summations involved in computing split in a lattice contribution for momenta , and a continuum contribution for , where and are the corresponding IR and UV cutoffs. The leading-order behaviour for small will thus be given by the leading-order in of the continuum contribution. Note that the gap acts as an effective IR cutoff, such that we can either evaluate the integral at and then take , or equivalently evaluate the momentum space integral at and then take the limit .

### b.1 Transversal coupling

Proceeding, we separate in Eq. (7) into the leading order and subleading terms. For the leading order, using explicit expressions for the Bogoliubov coefficients yields

 |Γαq|2+|Γβq|2=3(1+h)(|Γaq|2+|Γbq|2)−2|f(q)|R[ΓaqΓb∗q]ω(q). (30)

Expanding the denominator for small up to quadratic order, the continuum contribution to the coupling is of the form (setting )

 j(r)∼∫d2q(2π)2eiq⋅r(27hω(q)2+9(1−h)q2ω(q)2). (31)

The squared dispersion in the long-wavelength limit reads

 ω(q)2≃9h(h+2)+3/2q2. (32)

Counting powers of momenta, we find that the second term gives a regular contribution, while the first integral formally is -divergent. This divergence however is cancelled by the factor originating from expanding the vertex: Explicitly, the first integrand of the form can be evaluated analytically to yield a modified Bessel function (having absorbed numerical prefactors),

 j(r)∼h4πK0(hr2), (33)

such that no divergent terms appear as . Expanding the Bessel function, we thus obtain a scaling function

 jxy,(0)ij(h)∼jxy,(0)ij+Ahlogh+Bh+Ch3logh+Dh3 (34)

for . The prefactors depend on the distance to the respective neighbor, and can be determined by a non-linear fitting routine.

The analysis for the subleading terms in proceeds similarily. Crucially, we note that the vertices and in Eqs. (25b) and (25c) each contribute a constant factor (i.e. independent of ) that involves a momentum summation over the Bogoliubov coefficient. Using that and similarily we find that the leading-order contribution to the scaling due to the summation is of the form (note )

 ∫k<|Λ|d2k1√k2+h=−√h+√h+Λ2∼√h+O(h). (35)

Since these factors are -independent, they also multiply also a -independent terms appearing in the full evaluation of the subleading corrections, multiplying the scaling form of by a factor . Furthermore, we note that the self-energy at

 9(1+h)−|f(q)|2ωqh→0−−→ωq. (36)

The scaling behaviour of the self-energy correction to thus scales similarily to (34), such that in total the following scaling ansatz for the subleading corrections is assumed

 jxy,(1)ij(h)∼jxyij+A hlogh+B√h+Ch+D√h3. (37)

An example for the extrapolation for the subleading term with fit according to the scaling form given above is shown in Fig. 7

### b.2 Longitudinal coupling

We now discuss the scaling behavior of . Since the expression (8) involves odd powers of , which lead to non-analytic behavior (i.e. the continuum limit and do not commute), we consider the case of and work at a finite cutoff as discussed above. For simplicity we consider the first vertex in (8) and expand

 |Γαβq,k|2ωkωq+k∼(2k2+2k⋅q+q2−2|k||k