D T Liu

F J Burnell

L D C Jaubert

J T Chalker

Classical spin liquids in stacked triangular-lattice Ising antiferromagnets

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The triangular lattice Ising antiferromagnet is arguably the simplest model of a highly frustrated magnet and was probably the earliest such system to be studied in detail 1 . At low temperatures it is both highly fluctuating and strongly correlated; indeed, it remains disordered down to zero temperature and has a macroscopically degenerate ground state. The combination of fluctuations with correlations is typical more generally of highly frustrated magnets, which in this regime have been termed cooperative paramagnets or classical spin-liquids [START_REF] For | Highly Frustrated Magnetism[END_REF] .

In this paper we consider three-dimensional (3D) Ising antiferromagnets built from triangular layers that are stacked in such a way that nearest-neighbour interlayer interactions are frustrated, and make comparisons with the unfrustrated stacking. We focus on lowtemperature behaviour in systems with weak interlayer coupling, where correlations within each layer are necessarily strong but correlations between layers are controlled by a competition between fluctuations and interactions. Using a combination of perturbative and nonperturbative analytical techniques and Monte Carlo simulations, we show that this competition leads to a classical spin liquid regime, in which strong correlations exist without long range order.

Models for frustrated magnets can be classified at the mean-field level according to the properties of the matrix of exchange interactions. In this approach, the eigenvectors associated with the minimum eigenvalues of the interaction matrix provide candidate ordering patterns. These minimum eigenvalues appear at isolated points in reciprocal space for unfrustrated systems, but may be highly degenerate for frustrated systems. For example, for nearest neighbour interactions on the kagome and pyrochlore lattices, the subspace of minimum eigenvalues forms a flat band that spans the entire Brillouin zone [START_REF] For | Highly Frustrated Magnetism[END_REF][START_REF] Chalker | [END_REF][4] . Other cases display intermediate behaviour: on the diamond lattice with nearest and next-nearest neighbour interactions, the minimum eigenvalues form a two-dimensional surface in the 3D Brillouin zone 5 . The systems we discuss here are distinctive in having minimum eigenvalues that lie on lines in the 3D Brillouin zone 6 . One of our central findings is that these systems have a cooperative paramagnetic regime in which they develop strong correlations that are centred near these reciprocalspace lines.

The three different ways of stacking triangular layers that we compare in this work are indicated in standard notation by aaa, abc, and abab: see Fig. 1. Of these, the first provides a reference model without interlayer frustration, while the abc stacking yields minimum eigenvalues along helices in the Brillouin zone, and the abab stacking gives minimum eigenvalues on a ring around the Brillouin zone corner. The abc stacking with equal in-plane and interlayer interactions is equivalent to a nearest-neighbour model on the face-centered-cubic (fcc) lattice, while the abab stacking forms the hexagonalclose-packed (hcp) lattice.

Moving beyond a mean-field classification, the theoretical understanding of stacked triangular lattice Ising antiferromagnets (TLIAFMs) that we develop here is based on the height model description of low-temperature states for a single layer 7,8 . This long-established model represents ground states of a layer in terms of an emergent height field, with a simple effective Hamiltonian that captures the entropy of fluctuations. A spin-flip excitation fractionalises into an unbound vortex-antivortex excitation pair in the height field, and the vortex separation sets the correlation length at finite temperature in the single-layer model. In the following we derive and study height models for weakly coupled multilayer systems, showing how the interplay of interlayer coupling and vortex excitations allows strong correlations to develop between layers, without long-range order. We also use the results of extensive Monte Carlo simulations to test these conclusions and to examine behaviour when interlayer coupling is not weak.

Our study is motivated in part by observations 9,10 of charge ordering in the materials LuFe 2 O 4 and YbFe 2 O 4 .

FIG. 1. The three different ways of stacking triangular lattices that are considered in this paper: aaa (top left), abc (top right), abab stacking (bottom). In-plane interactions J and interlayer interactions J ⊥ are indicated with full and dashed lines, respectively.

The charge states of Fe 2+ and Fe 3+ ions in these systems can be represented using Ising pseudospins, with antiferromagnetic coupling between pseudospins arising from screened Coulomb interactions 9,11 . The pseudospins occupy the sites of an abc-stacked triangular lattice, though with an alternating layer spacing that is not included in the models we study. Experimental studies 9,10,12 , in particular of YbFe 2 O 4 10 , find helices of scattering intensity in a temperature range above a three-dimensional charge-ordering transition. These helices mirror in their reciprocal space location the positions of minimum eigenvalues of the interaction matrix discussed above. While an accurate description of these materials would require treating additional (magnetic) degrees of freedom 12 , the results we present in this paper demonstrate how strong interlayer correlations can arise over an extended temperature range without long-range order.

Past theoretical work on charge ordering in these materials has included quite detailed mean-field treatments 9,11 and Monte Carlo simulations of a bilayer model 13 , but has not made use of the understanding of single-layer TLIAFMs provided by height models, or used simulations to study correlations in the paramagnetic phase with the detail we present here.

TLIAFMs with other stackings have been examined previously in a variety of contexts. Treatments of the abab case include mean-field theory, a low temperature expansion, and Monte Carlo simulations [14][15][16] . That work has probed the ordering transition, but without examining the limit of weakly coupled layers or correlations in the paramagnetic phase. TLIAFMs with unfrustrated (aaa) stacking have been of long-standing interest 17 . They display a continuous phase transition that, strikingly, is in the 3D XY universality class despite the absence of a microscopic continuous symmetry 18,19 . The two components of the order parameter represent ordering at the two inequivalent Brillouin zone corners, and the XY symmetry is broken in the ordered phase by dangerously irrelevant six-fold anisotropies. This model and transition are also important as an imaginary time representation of the quantum dimer model on the hexagonal lattice 20 .

The remainder of the paper is organised as follows. We introduce the models studied and give an overview of their physical behaviour in Sec. II. We describe Monte Carlo results in Sec. III. We introduce height models in Sec. IV and analyse their behaviour in Secs. V and VI. Results from our different approaches are compared in Sec. VII. Some technical details are described in a series of appendices. An outline of some of the results has been presented previously in Ref. 21.

II. MODELS AND OVERVIEW

The starting point for our investigation is the nearest neighbour Ising antiferromagnet on stacked triangular layers with anisotropic couplings. Each spin is coupled to its six in-plane neighbours with an exchange constant J > 0 and to the closest spins in the layers above and below with an exchange constant J ⊥ (see Fig. 1). The Hamiltonian is

H = J ij ,z
σ i,z σ j,z + J ⊥ {ij},z σ i,z σ j,z+1 + H (1) (1) where H (1) indicates further-neighbour interactions, which may be present in the bare Hamiltonian or may represent terms generated under renormalisation. Here σ i,z = ±1, the notation i, j denotes nearest neighbour pairs of sites from the same layer, and {i, j} nearest neighbour pairs from adjacent layers. The sign of J ⊥ may be taken positive without loss of generality, since it can be reversed by the transformation:

σ i,z → σ i,z = (-1) z σ i,z .
We are concerned with the statistical mechanics of these models as a function of temperature T and the interaction strength ratio J ⊥ /J. At J ⊥ /J = 1, one expects ordering below a temperature T c ∼ J, while for J ⊥ /J = 0 the system of uncoupled layers remains disordered at all temperatures. A schematic phase diagram obtained by interpolating between these limits has the form shown in Fig. 2. For J ⊥ /J 1 the paramagnetic phase extends to temperatures T J. In this regime, spins are highly correlated within each layer. Our objectives are to understand interlayer correlations and the form of the phase boundary for small T /J and J ⊥ /J, in each of the three stackings. For the two frustrated stackings we find that at small J ⊥ /J there is a low-temperature regime in which the correlation lengths, both in-layer and interlayer, are much larger than the lattice spacing. A system in this regime is termed a cooperative paramagnet or classical spin liquid. This regime is smoothly connected to the conventional paramagnetic state at T J but distinguished from it by strong correlations.

For orientation it is useful to have a simple approach that gives an initial indication of likely behaviour. Mean field theory can often be employed in this way but fails here, wrongly predicting an ordering temperature set by J, even for small J ⊥ . An alternative that has been widely applied in geometrically frustrated magnets is the selfconsistent Gaussian approximation (SCGA) 22 . It is wellcontrolled only for n-component spins at large n, but is known in some instances to be quite accurate even for Ising systems 23 . In the SCGA, correlations are given in terms of the interaction matrix J and the inverse temperature β by

σ i σ j = (βJ + λI) -1 ij . (2)
Here, λ is a parameter fixed by the consistency condition |σ i | 2 = 1, which can be satisfied throughout the paramagnetic phase. Using a spectral decomposition of J in terms of its eigenvalues l q and eigenvectors u l q (α), where α labels sites within a unit cell and l labels the bands of J, the SCGA expression for the structure factor is

S(q) = 1 N i,j (βJ + λI) -1 ij e iq•(ri-rj ) = l,α,α u l * q (α) u l q (α ) β l q + λ (3)
From this it is apparent [barring cancellations in the sum α,α u l * q (α) u l q (α )] that maxima in S(q) arise from minima in l q . Applying the SCGA to stacked triangular lattice antiferromagnets, the paramagnetic phase extends to temperatures T J if J ⊥ J, and in this regime the maxima in S(q) are sharply defined. To find the location of these maxima in reciprocal space, we examine the minima of l q . We take axes with ẑ perpendicular to the triangular layers, unit spacing between neighbouring layers for the aaa and abc stackings, and unit spacing between neighbouring a-layers in the abab stacking, which has two sites per primitive unit cell. We choose in-plane lattice vectors

a 1 = (1, 0, 0)
and

a 2 = (1/2, √ 3/2, 0), (4)
The corresponding in-plane reciprocal lattice vectors are A 1 = 2π(1, -1/ √ 3, 0) and A 2 = 2π(0, 2/ √ 3, 0). We use δ to denote the separation vector between neighbouring sites in adjacent layers. Hence δ = (0, 0, 1), (1/2, 1/(2 √ 3), 1) and (1/2, 1/(2 √ 3), 1/2) for the aaa, abc and abab stackings, respectively.

The contribution to l q from in-plane interactions has a minimum at the K-points of the triangular lattice Brillouin zone:

K = ( 4π 3 , 0) and K = ( 2π 3 , 2π √ 3 ). ( 5
)
Upon inclusion of small J ⊥ , these minima evolve in different ways for each of the stackings we consider. For the aaa stacking, they lie at isolated points, undisplaced inplane and at q z = π. For the frustrated stackings, their locations can be specified in terms of the wavevectordependent complex scalar ζ = 1 + e iq•a1 + e iq•a2 . In the abc case they lie on the curve

ζ = - J ⊥ J e iq•δ (6)
and in the abab case they lie on

|ζ| = J ⊥ J , q z = 0 . (7)
These conditions respectively define helices and rings centred on the zone corners, as shown in Fig. 3. Further discussion of the interaction matrix eigenvalues is given in Appendix A.

III. MONTE CARLO SIMULATIONS

We use extensive Monte Carlo simulations to find the ordering temperature for all three models and to study correlations in the paramagnetic phase of models on the abc and abab stacked lattices. The primary observables computed are the energy E, specific heat C, and the structure factor S(q), which is obtained from the Fourier transform of magnetisation

σ(q) = i e iq•ri σ i (8)
as Because of the complex energy landscape arising from geometrical frustration, we employ a parallel tempering algorithm with single-spin-flip Metropolis dynamics 24,25 . Specifically, we simulate N r replicas (taking N r ∼ 100) at geometrically spaced temperatures, with the highest temperature ∼ 5J. A Monte Carlo sweep involves one single-spin-flip attempt per site, followed by one parallel tempering swap attempt between replicas at adjacent temperatures. A system consists of L z rhombic layers, each of size L × L lattice constants, with periodic boundary conditions in all directions. A typical simulation treats ≈ 10 5 sites (L = 72 -200, L z = 12 -48) using 10 5 sweeps. We measure E and C each sweep, and S(q) every N r sweeps. Further details of the data analysis are presented in Appendix B.

S (q) = 1 L 2 L z |σ(q)| 2 . (9)

A. Ordering Transition

Phase diagrams as a function of T and J ⊥ are shown in Fig. 4 for both the unfrustrated (aaa) and the frustrated (abc and abab) stackings. For a given strength of interlayer coupling, the ordering temperature (determined from the maximum of the heat capacity) is much lower in the systems with frustrated stackings compared with the unfrustrated one. In addition, over most of the range of J ⊥ /J studied, the transitions in the systems with frustrated stackings are strongly first order: the probability distribution of the energy is strongly bimodal at the transition unless J ⊥ /J 1. The discontinuity in the energy at the transition decreases with decreasing J ⊥ , and for J ⊥ 0.05J the order of the transition is not discernible from the simulations. Differences in transition temperature between the two frustrated stackings are very small for J ⊥ /J ≤ 1. Our results for the abc stacking at J ⊥ = J can be compared with earlier work on the fcc lattice, and are in good agreement with the transition temperature of T c ≈ 1.72J found in Refs. 26 and 27. Examples of the energy distribution at different temperatures are shown in Fig. 5. We monitor the overlap of distributions at adjacent temperatures in the parallel tempering scheme, as substantial overlap is a requirement for effective exchange of replicas. The top panel demonstrates that this is the case in our simulations. At a firstorder transition, the energy distribution is bimodal. The middle panel illustrates this. Finite size shifts in our estimates of the transition temperature are a few percent, as indicated by a comparison of the middle and lower panels.

B. Correlation functions

A characteristic feature of classical spin liquids is the presence of strong correlations and a large correlation length, without long-range order or proximity to a critical point. In this subsection we present correlation functions and correlation lengths for TLIAFMs with frustrated stackings, determined from Monte Carlo simulations.

The abc stacking

The behaviour of the structure factor for a system with abc stacking in the classical spin-liquid regime is illustrated in Fig. 6. Combining information from the series of slices in reciprocal space that are shown in this figure, it is apparent that maxima in S(q) lie on helices in reciprocal space. The axes of these helices pass through corners of the triangular-lattice Brillouin zone.

To analyse this behaviour quantitatively, we extract a reciprocal-space radius Q for the helix and a correlation length ξ ⊥ by fitting data for S(q) near the maxima to a sum of in-plane Lorentzians S (q) = I ξ 2 ⊥ (q ⊥ -q 0 ⊥ (q z )) 2 + 1 (10) from each helix. Provided any dependence of |q 0 (q z )| on q z is weak, we can make the identification Q = |q 0 (q z )|.

(See Appendix B for further discussion.) Results are shown in Fig. 7. The correlation length ξ ⊥ increases rapidly with decreasing T for T J, as demonstrated in Fig. 7(a). It reaches large values within the paramagnetic phase if J ⊥ /J is small. Its dependence on J ⊥ at fixed T is very weak, because its value is determined by the density of vortices in the height field [see Sec. V] and for J ⊥ J this in turn is controlled mainly by the value of T /J. The variation of the helix radius Q with J ⊥ and T is illustrated in Fig. 7(b). Its value is given quite accurately by the SCGA, Eq. ( 6), for T J, and shows a small increase with decreasing temperature.

In the ordered phase, Bragg peaks are expected in the structure factor, in place of a continuous distribution of weight on helices. We probe the evolution between the two behaviours by computing

S avg (q z ) = 1 L 2 qx,qy S(q) . (11)
Results in Fig. 8 show the rapid development of Bragg peaks as temperature is lowered through the transition.

Although we believe that the transition is first order for the value of J ⊥ /J studied here, discontinuities are not apparent in the temperature dependence of S avg (q z ), presumably because of finite-size rounding. Indeed, since evaluation of correlation functions is more computationally demanding than calculation of energy distributions, the results presented in Fig. 8 are for smaller system size than those in Fig. 5; we find (data not shown) that the energy distribution at the transition is not bimodal for the smaller size.

The abab stacking

Because the abab-stacked lattice has two sites in a primitive unit cell, the relation between fluctuations and correlations is less direct than for the abc stacking, in which the unit cell has a single site. More specifically, the form of S(q) is affected by interference between contributions from the two sites. Within the SCGA, this is apparent from Eq. (3), where contributions involving a given eigenvalue l q of the interaction matrix are weighted by a sum α,α u l * q (α) u l q (α ) that includes both sitediagonal (α = α ) and interference (α = α ) terms. In order to eliminate these interference effects and expose fluctuations in the abab stacking in a simple way, we compute the structure factor using contributions only from

q z = 0 q x 0 2 4 6
-2 0 2 4 6 q y q z = 1.3 q z = 2.6 q z = 3.8 q z = 5.1 q z = 6.4 FIG. 6. Cross-sections of structure factor at constant qz in a system with abc stacking. For each qz, sharp maxima in S(q) occur near the Brillouin zone boundary, which is shown as a green dashed line. As qz increases, the maxima precess around the zone corners without significant change in intensity, indicating that they form helices in the three-dimensional reciprocal space. Parameter values are J ⊥ = 0.2J, T = 0.8J, L = 72, Lz = 12; for this value of J ⊥ , Tc = (0.68 ± 0.01)J.

one of the two sites in each unit cell, by restricting the sum in Eq. ( 8) to this set of sites. We expect from Eq. ( 7) that this single-sublattice structure factor will have its maxima lying on closed loops in the q z = 0 plane. An overview of our data, illustrating this behaviour, is given in Fig. 9.

A simple way to extract a correlation length ξ ⊥ is by fitting data for q z = 0 and q x , q y close to a selected Brillouin zone corner to the functional form

S (q) = I ξ 2 ⊥ (Q -|q ⊥ -K|) 2 + 1 , (12)
where K denotes the location of the Brillouin-zone corner and Q specifies the reciprocal-space radius of the ring of intensity. This fitting function provides a good description of the data for small values of J ⊥ /J, where the maximum in the structure factor lies on a circle, but it does not capture the triangular distortions for larger J ⊥ /J that are apparent in the left-most panel of Fig. 9.

As shown in Fig. 10, and as for the abc stacking, the resulting values of ξ ⊥ increase rapidly with decreasing temperature but vary little with J ⊥ .

Self-consistent Gaussian Approximation

As discussed in section II, the SCGA provides a useful description of frustrated magnets in the strongly correlated regime. In particular, it offers a simple theoretical prediction for S(q), which we now show to be a good representation of our simulation data. We use the functional form of Eq. (3) in two ways, which are distinct in principle but yield very similar results. One of these treats the variable λ as a fitting parameter with respect to simulations; the other fixes its value using the SCGA condition

|σ i | 2 = 1.
The SCGA form for S(q) is especially helpful at larger values of J ⊥ /J, when detailed lattice effects are impor-tant. The results of these lattice effects for the abc stacking include a dependence of the helix radius [q 0 ⊥ (q z ) in Eq. ( 10)] on q z . For the abab stacking they generate correlations that are not represented using the circular maximum in S(q) implied by the fitting function given in Eq. (12). The SCGA gives a good description of this physics. Most notably, for the abab stacking the SCGA fits are effective in capturing the triangular distortion of the rings, as demonstrated in Fig. 11.

Once the value of λ is obtained from the fit, the correlation length can be extracted from the model. The results for ξ ⊥ are shown in Fig. 12. They agree to ∼ 10% with those obtained by fitting the functional forms given in Eq. (10) and Eq. ( 12) for the abc and abab cases respectively (see Figs. 7a and10). Alternatively, the value of λ can be determined without reference to simulations, using the SCGA condition, yielding a theoretical prediction for ξ ⊥ . From Fig. 12, it is apparent that both approaches to determining λ yield very similar results. 6). Results were obtained in a system of size L = 72, Lz = 12. Data for each value of J ⊥ extend to the lowest temperature employed in parallel tempering that was above Tc.

IV. HEIGHT MODEL

h = 1 3 (hA + hB + hC ) hA hB hC σA σB σC 0 0 1 5 + -- 1 0 1 2 + -+ 2 3 1 2 --+ 3 3 4 2 -+ + 4 3 4
We now turn to an analytical treatment of stacked triangular lattice Ising antiferromagnets. Although the SCGA, as demonstrated, provides a good approximate description, it is formally correct only for n-component spins in the large-n limit. It is therefore not a natural starting point for a systematic approach. By contrast, the height model provides a representation of a single-layer TLIAFM that is known to capture exactly the physics at low temperatures and long distances. Here we use the height model to construct a description of the multilayer system that allows for a controlled treatment of weak interlayer interactions.

Following Blöte et al. 7 and Zeng and Henley 8 , we map ground states of a single layer Ising model onto states of a height model in such a way that spin configurations with long-range three-sublattice order correspond to flat height configurations. Because of frustration, domain walls can be introduced without energy cost between regions with different types of three-sublattice order. These domain walls correspond to steps in the height field. In a coarse-grained description, steps are represented by a gradient in the height field, and a large value for this 11)] as a function of qz at four selected temperatures near the transition, in a system with J ⊥ = 0.1J. Data are for T = 0.45J, 0.5J, 0.52J and 0.6J, in order of decreasing peak intensity, and the transition temperature is Tc ≈ 0.54J. Bottom: Savg(qz) as a function of T for the three values of qz that are marked with vertical dashed lines in the top panel. Results for both panels were obtained in a system of size L = 36, Lz = 48.

gradient carries an entropy penalty. The mapping is conveniently described in two stages. First we define heights at the sites of the triangular lattice, as in Ref. 7. Second, following Ref. 8, we average these site heights to define heights at the centres of triangles, obtaining a height model that is easily coarsegrained.

To map from a spin configuration to heights at lattice sites, we first assign height zero to a reference site. The heights on all other sites of the lattice are then fixed by the requirement that the height difference between the neighbouring sites i and j is +2 if σ i = σ j , and -1 if σ i = -σ j going anticlockwise around an up-triangle (or clockwise around a down triangle): see Fig. 13a. Heights at triangle centres are defined as the averages of site heights at vertices. The advantage of this locallyaveraged height field is that ground states with threesublattice order are exactly flat in these variables: see Fig. 13c. In the following, we use the term 'height field' exclusively for the locally-averaged quantity.

This mapping is summarised for a single triangle in FIG. 9. Cross-sections of structure factor at qz = 0 for systems with varying J ⊥ in the abab stacking. Intensity is maximum on a closed loop, which is approximately circular for small J ⊥ /J but develops triangular distortions with increasing J ⊥ /J. Data (from left to right) are for T = 1.14J, 0.71J, 0.64J, 0.57J, obtained in systems of size L = 72, 90, 90, 204 and Lz = 12, 12, 30, 6. Note the changing intensity scale and increasing maximum intensity as J ⊥ and T decrease. The ordering temperatures are Tc/J = 0.99 ± 0.008, 0.680 ± 0.014, 0.602 ± 0.007 and 0.502 ± 0.01. Table I. Here the sites of the triangular lattice are divided into three sublattices, labelled A, B, C and indicated by the three colours of dots at the vertices in Fig. 13b. With this convention, the six ground states of each triangle are specified by the orientation of the spin on sublattice A and the location of the frustrated bond. The groundstate spin configuration of a triangle fixes the value of the height h at its centre modulo 6.

J ⊥ = 0.4J q x 2 3 4 5 6
The mapping is unique up to labelling conventions. Permuting the choice of A, B and C sublattices (which results from lattice translations or rotations by 2π/3 about the centre of a triangle) corresponds to a global shift h → h + 2. (By contrast, rotations about an axis passing through a site leave the labelling and hence the height field invariant.) Shifting h → h + 3 corresponds to a global spin flip operation. The remaining possibilities (shifting h by 1 or 5) correspond to a combination of the global spin-flip and re-assignment of the three sublattices.

The inverse mapping, from a height configuration to a spin configuration, can be expressed in terms of a function f (h) and a constant s α . The function f (h) ≡ f (h+6) takes the values f (h) = +1 for h = -1, 0, 1 and f (h) = -1 for h = 2, 3, 4. The constant s α takes values s A = 0, s B = 2 and s C = -2 on sublattices α = A, B or C. The spin orientation is then given by

σ α = f (h + s α ) ≡ f α (h). ( 13
)
For integer h we can represent this function as f (h) = 4 3 cos πh 3 -1 3 cos πh. Note that since each spin is part of six triangles, to fully specify the mapping we must choose which triangle's height dictates which spin. Reassuringly, one can verify that this choice is unimportant: when the height configurations are integers, and can change by at most 1 between any pair of adjacent triangles, every convention yields the same spin configuration.

Excitations of the spin model consist of triangles in which all spins are up, or all are down. They are represented by vortices in the height field, which is multivalued in their presence: it increases by 6 on going anticlockwise around an upward-facing excited triangle, and decreases by 6 around a down-facing triangle. An excited state produced from a ground state by reversing a single spin necessarily contains a vortex-antivortex pair, which may be separated by additional spin flips without further energy cost.

H = K 2 d 2 r |∇h(r)| 2 + d 2 r Ṽ (h) . (14)
We can determine the value of K (and verify that (14) captures the correct physics) by comparing the correlation functions of this model with Ṽ (h) = 0 to those of the exact solution for the 2D TLIAFM. Stephenson 28 has shown that at long distances

σ α (r)σ β (r ) ∼ ω s |r -r | + c. c. , (15)
where s = (s α -s β )/2 and ω = e i2π/3 . The dominant terms in the expression for the intra-sublattice spin-spin correlation function in terms of the height fields are

σ α (r)σ β (r ) ∼ e i π 3 [h(r)-h(r )] ω s + c. c. ∼ exp - 2π 36K ln |r -r | (ω s + ω -s ) ∼ |r -r | -2π 36K (ω s + ω -s ) (16)
Hence at zero temperature, to reproduce the longwavelength properties of the exact solution, we take

K = π/9.
What about the potential term, which we ignored in the above calculation? Microscopically the heights are integers; we can account for this by including the potential Ṽ (h) = -v cos(2πh). At short distances v is large and positive. At longer length scales the effective value of v is determined by the scaling dimension of the operator cos 2πh, which can be deduced from the 2-point function

cos(2πhr) cos(2πhr ) ∼ |r -r | -2π K implying d 2 r cos(2πhr ) ∼ L 2-π K . (17)
This yields the scaling dimension 2 -π K = -7 at T = 0; hence the effective value of the coefficient v decreases rapidly as we probe the system at longer length-scales, and its effect on the long-wavelength correlations is negligible.

Finally, we can ask about behaviour at finite temperature. To describe the system at finite temperature we must include the possibility of vortices in the height field. Dropping Ṽ (h) in Eq. ( 14) but including vortices, we recover the physics of the 2D xy model at an effective temperature that is set by the value of K. The scaling dimension of the vortex can be computed by estimating its free energy: for v = 0 the entropic cost of the gradients in the height field required to insert a single vortex into a triangular layer of side length L is δH = 9K π ln L/a, where a is the lattice constant. The number of ways to place the vortex in the system is L 2 /a 2 . Together, these contributions to the free energy of a single vortex are

δF = 9K π -2 ln L a . ( 18
)
For K = π/9 this grows more negative with increasing L. We are therefore in the high-temperature phase of the xy model, where vortices are unbound. The vortex density, determined by the fugacity associated with the vortex excitation energy 4J, sets the correlation length. This reflects the fact that the triangular layer, which is critical at T = 0, is a paramagnet at any finite temperature. Hence the height model ( 14) correctly reproduces the phase diagram and correlations of an isolated triangular layer. The potential Ṽ (h) is an irrelevant operator and can be dropped from the long-wavelength analysis; however the vortices arising at finite temperature are relevant, making the system paramagnetic for any T > 0.

B. Coupled layers in the height model description

We now turn to the situation of interest, in which spins in triangular layers are coupled to their nearest neighbours in the planes directly above and below. We will derive expressions for these couplings in the height language, and discuss their effect on the physics of the system.

Frustrated interlayer coupling favours domain walls in the three-sublattice order that is represented by flat configurations of the height field. To minimise the interlayer exchange energy, these domain walls should stack in such a way that a domain wall consisting of up spins sits in the adjacent layer to a domain wall consisting of down spins, as shown in Fig. 14. To find the functional form of the interlayer coupling in height language, we use Eq. ( 13) to express it in terms of the height fields. We then find the scaling dimensions of the various contributions to determine which of these play an important role in the long-wavelength physics. We will show that, as in the SCGA treatment, for frustrated stackings the relevant terms in the nearestneighbour model lead to one-parameter sets of degenerate ground states in the height models, whose symmetry can be broken by including further-neighbour couplings.

!! ! !! ! !! ! !! !

Unfrustrated stacking

It is instructive to begin by studying the unfrustrated stacking. For the aaa stacking, the interlayer coupling is

J ⊥ (σ A,z σ A,z+1 + σ B,z σ B,z+1 + σ C,z σ C,z+1 ) = 8J ⊥ 3 cos π 3 (h z+1 -h z ) + J ⊥ 3 cos πh z cos πh z+1 + . . . (19)
where . . . represents terms of quadratic and higher order in the derivatives, which we drop as they are irrelevant in the scaling sense. The most relevant term is cos π 3 (h z+1h z ), which has a scaling dimension of 3/2 for K = π/9. The term cos πh z cos πh z+1 has scaling dimension -5/2 and can be neglected. Hence the effective Hamiltonian of the height model for the aaa stacking is

H (aaa) = K 2 z d 2 r |∇h z (r)| 2 +κ 3 cos π 3 (h z+1 -h z ) , (20)
with κ 3 = 16βJ ⊥ /3K. The ground states

h z (r) = γ (21)
of this effective model have a U (1) symmetry under changes of the constant γ. This symmetry is broken down to a six-fold discrete symmetry by the interaction Ṽ (h), which is irrelevant in the scaling sense at the fixed point describing uncoupled layers, and dangerously irrelevant at the three-dimensional ordering transition 19 .

Frustrated stackings

For both the abc and the abab stackings, we consider two neighbouring layers as shown in Fig. 15. There is a coupling between each site on the black lattice and the three sites around it from an up-triangle on the red lattice, or equivalently between each site on the red lattice and the three sites around it from a down triangle on the black lattice. We denote heights on the black lattice by h z+1 (r), and ones on the red lattice by h z (r). The coupling is

H ⊥ = J ⊥ r∈A σ A (r)[σ a (r + e 1 ) + σ b (r + e 2 ) + σ c (r + e 3 )] + symmetry-related terms = J ⊥ r∈A f A (h n+1 (r))[f a (h n (r + e 1 )) + f b (h n (r + e 2 )) + f c (h n (r + e 3
))] + symmetry-related terms, (22) where 'symmetry-related terms' have B or C in place of A, and a corresponding permutation of the vectors e i . These are defined in terms of the lattice vectors [Eq. (4)] by e 1 = 2 3 a 2 -1 3 a 1 , e 2 = 2 3 a 1 -1 3 a 2 and e 3 = -1 3 a 1 -1 3 a 2 , and are illustrated in Fig. 15 . Expanding h(r) in a Taylor series, we obtain

H ⊥ = - 4πJ ⊥ 9 √ 3 r cos π 3 (h z+1 (r) -h z (r))∂ x h z (r) -sin π 3 (h z+1 (r) -h z (r))∂ y h z (r) + . . . (23)
where . . . indicates RG-irrelevant terms. Thus keeping only the relevant inter-layer couplings leads to the effective Hamiltonian for the abc stacking

H (abc) = K 2 z d 2 r ∂ x h z -κ ⊥ cos π 3 (h z+1 -h z ) 2 + ∂ y h z + κ ⊥ sin π 3 (h z+1 -h z ) 2 - κ ⊥ K 2 (24)
with κ ⊥ ∝ βJ ⊥ .

For the abab stacking, the derivation is identical except that the vertical unit cell contains two layers, with the layers above and below offset in opposite directions. We use integer z to label unit cells in the vertical direction and µ = 1, 2 to label layers within each unit cell. The effective Hamiltonian is

H (abab) = K 2 z d 2 r µ |∇h z,µ | 2 -κ ⊥ ∂ x (h z,1 + h z,2 ) cos π 3 (h z,2 -h z,1 ) -∂ y (h z,1 + h z,2 ) sin π 3 (h z,2 -h z,1 ) +∂ x (h z+1,1 + h z,2 ) cos π 3 (h z,2 -h z+1,1 ) -∂ y (h z+1,1 + h z,2 ) sin π 3 (h z,2 -h z+1,1 ) . ( 25

C. Symmetries and further-neighbour couplings

For both frustrated stackings, emergent continuous symmetries not present in the lattice models are displayed by the effective Hamiltonian of Eqns. ( 24) and (25) if terms irrelevant at the J ⊥ = 0 fixed point are omitted. Both models have a U (1) × U (1) symmetry. One U (1) symmetry is associated with global shifts in the height field. It results from the discrete symmetry of the microscopic model related to global shifts in h, which -as for the single-layer height model -is enhanced to become a continuous symmetry because the pinning potential Ṽ (h) is RG-irrelevant and has been omitted. As in the unfrustrated case [see Eq. ( 21)] we parameterise it with γ. The second U (1) symmetry is associated with real-space rotations and is reduced to the discrete rotational symmetry of the lattice by irrelevant terms. We parameterise it with θ.

In detail, these symmetries take the following form. Let R θ denote a rotation in the xy plane through the angle θ and write r = R θ (r). Then H (abc) is invariant under the transformation

h z (r) → h z (r) = h z (r ) + 3zθ π + γ . ( 26
)
Similarly

H (abab) is invariant under h z,µ (r) → h z,µ (r) with h z,1 (r) = h z,1 (r ) - 3θ 2π + γ and h z,2 (r) = h z,2 (r ) + 3θ 2π + γ . ( 27
)
Ground state configurations of the height model for the abc stacking have the form

h z (r) = κ ⊥ (x cos θ -y sin θ) + 3zθ π + γ . ( 28
)
For the abab stacking the ground states are

h z,1 (r) = κ ⊥ (x cos θ -y sin θ) - 3θ 2π + γ h z,2 (r) = κ ⊥ (x cos θ -y sin θ) + 3θ 2π + γ , (29)
together with a second symmetry-related set.

The symmetry under continuous changes of θ is not a feature of the microscopic model: it is broken by the leading irrelevant terms in Eq. (23). For the abc stacking these have the form

H b = κ b z d 2 r (∂ x h z (r)) 2 -(∂ y h z (r)) 2 cos δh z (r) + 2∂ x h z (r)∂ y h z (r) sin δh z (r) , (30)
where we introduce the notation δ p h z (r) = π 3 [h z+p (r)h z (r)] and δh z (r) ≡ δ 1 h z (r). (The form for the abab stacking follows the obvious equivalent pattern.)

Significantly, it may also be broken by relevant furtherneighbour couplings, if these are present microscopically, or are generated under renormalisation. For the abc stacking, some relevant and marginal couplings that are not included in Eq. ( 24) are

H m = K m 2 z d 2 r ∇h z (r) • ∇h z+m (r) H 2 = κ 2 z d 2 r ∂ x h z (r) cos π 3 (h z+2 -h z ) +∂ y h z (r) sin π 3 (h z+2 -h z ) H 3 = κ 3 z d 2 r cos π 3 (h z+3 -h z ) . ( 31
)
H 3 is the most relevant of these three: it breaks the degeneracy of Eq. ( 28), selecting ground states for which 3θ = 0 (π) for κ 3 < 0 (κ 3 > 0). H 2 has the same scaling dimension as the bare interlayer coupling. It also breaks the symmetry, again favouring states for which 3θ = 0 (π) for κ 2 < 0 (κ 2 > 0). H m is marginal, and does not break the degeneracy between the ground states identified above, all of which have the same in-plane gradients in each layer. Therefore as well as potentially being broken spontaneously at low temperature, the emergent U (1) spiral symmetry of the abc model can be broken explicitly at a scale set by the coefficients κ 2 and κ 3 . We discuss this scenario in Sec. V.

For the abab stacking, the perturbations of interest are interlayer gradient couplings similar to H m , and also

H 3 = κ 3 z,µ d 2 r cos π 3 (h z+1,µ -h z,µ ), (32)
the unfrustrated coupling between spins two layers apart.

In contrast to the abc case, H 3 is not expected to be important in determining the ordering temperature: the minimum-energy solutions of the abab model have a definite value of h z+1,µ -h z,µ , and so this term does not lift the ground-state degeneracy. Instead, symmetry is broken by the irrelevant coupling H b , Eq. ( 30).

V. BEHAVIOUR OF THE HEIGHT MODEL

To understand the phase diagrams of these coupledlayer height models, we take two successive steps. First we make a perturbative renormalisation group (RG) analysis of the behaviour of weakly coupled layers, as described in Sec. V A. Depending on the values of T and J ⊥ , the model under scaling may remain weakly coupled: this happens in the weakly-correlated paramagnetic regime. Alternatively, it may flow to strong interlayer coupling. In that case a separate analysis is necessary of the influence of vortex pairs, which is presented in Sec.V B. We find that the minimal models with exact U (1) × U (1) symmetry have anomalously soft excitations. For this reason vortex pairs destroy long-range order, establishing instead a paramagnetic regime with strong interlayer correlations. Symmetry-breaking or 'locking' interactions act in competition to vortex pairs, and stabilise the ordered phase when they dominate.

A. Perturbative RG

Our perturbative analysis follows the standard renormalisation-group techniques of Refs. 31 and 38. For small J ⊥ and low T , this allows us to use arguments similar to those of Sec. IV A regarding the phase diagram of these models. If unbound vortices proliferate, the interlayer coupling flows to zero at long distances, while if the coefficient of one of the cosine terms grows large, a strong-coupling analysis is necessary.

The leading-order behaviour of the RG equations is simply determined by the scaling dimensions of the relevant interlayer couplings and vortices. (The intra and interplane gradient terms flow only at higher order.) For the interlayer couplings, these can be calculated either from the two-point functions as described in Sect. IV A, or (as is more appropriate for operators involving derivatives of the height field) using a standard momentumshell RG (see Appendix C 1). Using to denote the shortdistance cut-off and following the notation of Eqns. ( 24), ( 25) and (31), this gives

∂κ ⊥ ∂ ln = (1 -β 1 ) κ ⊥ , ∂κ 3 ∂ ln = (2 -β 1 ) κ 3 and ∂y ∂ ln = (2 -α 1 ) y . (33)
Here, κ ⊥ is the frustrated interlayer coupling that acts between neighbouring layers in the abc and abab stackings, and κ 3 is the unfrustrated inter-layer coupling, which couples nearest neighbour layers in the aaa stacking, second neighbours in the abab stacking and third neighbours in the abc stacking. Finally, y is the vortex fugacity, which dictates the unbound vortex density. For weakly coupled layers we have

β 1 = π 18K and α 1 = 9K π . (34)
For the unfrustrated stacking the bare value of κ 3 is κ 3,0 ∼ βJ ⊥ . For the frustrated stackings the bare value of the interlayer coupling κ ⊥ is κ ⊥,0 ∼ βJ ⊥ . In both cases, the bare value of the vortex fugacity is y 0 ∼ e -4βJ . The initial value of is the lattice spacing, which we set to unity.

Let us now consider what we learn from these scaling dimensions about behaviour in the three different models, keeping only nearest-neighbour interactions and the intralayer gradient interaction K. Using the value K = π 9 appropriate for decoupled triangular layers, we have α 1 = 1, β 1 = 1/2, and single-layer vortices are more relevant than their multi-layer counterparts. Solving the RG equations (33) gives

y = y 0 , κ ⊥ = κ ⊥,0
1/2 and κ 3 = κ 3,0 3/2 .

The calculation reaches its limit of validity at the scale where the largest coupling is of order unity, and the physical state of the system is signalled by which coupling first crosses this threshold. If y ∼ 1 with κ ⊥ and κ 3 1, the system is a weakly correlated paramagnet. If either κ ⊥ ∼ 1 or κ 3 ∼ 1 with y 1, layers are strongly coupled. We turn next to this regime.

B. Strongly coupled layers

To understand behaviour of the height models at large interlayer coupling, we examine the effective Hamiltonian for each type of stacking at quadratic order in an expansion about the ground states given in Eqns. ( 21), ( 28) and (29).

For orientation, consider first the aaa stacking. Let ϕ z (r) denote the deviation of h z from a ground-state configuration and introduce its Fourier transform via ϕ z (r) = 1 (2π) 3 d 3 q ϕ(q)e i(q ⊥ r+qzz) .

The energy cost at quadratic order of this deviation from a ground state is

δH = K 2(2π) 3 d 3 q E(q)|ϕ(q)| 2 (36)
with

E(q) = q 2 x + q 2 y + κ(1 -cos q z ), (37)
where κ⊥ = (π 2 /9)|κ 3 |. Thus, for this unfrustrated stacking, excitations have a dispersion E(q) that is conventional in the sense that it is quadratic in wavevector for all orientations of q.

An equivalent calculation for the abc stacking (for fluctuations around the ground state with θ = 0) yields the quite different dispersion relation

E(q) = q 2 x + (q y -κ⊥ sin q z ) 2 + κ2 ⊥ (1 -cos q z ) 2 , (38)
where κ⊥ = (π/3)κ ⊥ . This is anomalously soft, being quartic in wavevector along the line q y = κ⊥ q z . The soft modes do not give rise to divergent harmonic fluctuations, since

[h n+1 (r) -h n (r)] 2 = 1 K d 3 q (1 -cos q z ) 2 E(q) (39)
is finite provided κ ⊥ = 0. For the abab stacking, since there are two layers within a unit cell, it is necessary to introduce two fields ϕ z,µ (r), with µ = 1, 2. The resulting quadratic Hamiltonian has two eigenvalues, which for θ = 0 are

E ± (q) = q 2 x +q 2 y +2κ 2 ⊥ ±2κ ⊥ | cos(q z /2)| q 2 y + κ2 ⊥ . (40)
In this case as well, the dispersion relation is quartic for one direction, since E -= q 2 x +(κ 2 ⊥ q 2 z +q 4 y /κ 2 ⊥ )/4 for small |q|, but harmonic fluctuations are bounded for κ ⊥ = 0.

C. Destruction of order by defects

Our discussion of harmonic height-field fluctuations around ground states of the multilayer model accounts for spin fluctuations within the ground-state manifold of each triangular layer, but a separate treatment is required to understand the effect of excitations out of this ground-state manifold. That is the subject of this subsection.

The excitations are represented by vortices and antivortices. These are unbound in a single layer, as discussed in Sec. IV A, but acquire a linear confining potential within ordered states of the multilayer systems. More specifically, suppose that the height field in a layer containing a widely separated vortex-antivortex pair has a step of height 6 and width w: its energy cost per unit length is ∼ Kw(w -2 +κ 2 ⊥ ) and is minimised by the choice w ∼ κ -1 ⊥ . Pairs are therefore bound with typical separation w when interlayer correlations are strong. Remarkably, although in other settings bound vortex pairs are typically irrelevant at large scales, we find that they exert a controlling influence in multilayer height models with frustrated stackings.

Height fields in the presence of vortices are in general multivalued, but can be taken to be single-valued in a domain that excludes a core around each vortex-antivortex pair. The presence of these pairs influences the height field far from the cores. A convenient alternative to an explicit treatment of multivalued height fields is to impose a potential that couples linearly to the height field and has the same effect on the far field as a votex-antivortex pair. In order to demonstrate the required form of this potential, consider a single layer containing a pair centred at the origin with separation vector b. This pair is described by the height field configuration

h(r) = 3 π arctan 2x + b • x 2y + b • ŷ -arctan 2x -b • x 2y -b • ŷ .
For |r| |b| we have

h(x, y) ≈ 3 π ẑ • (b × r) r 2 (41)
or equivalently

h(q) ≈ 6i ẑ • (q × b) q 2 . ( 42
)
The same far-field height configuration can be induced by adding a potential term v(q) to the effective Hamiltonian for the height field. Specifically, for an isolated layer, the effective Hamiltonian (K/[2π] 2 ) d 2 q 1 2 E(q)|ϕ(q)| 2 -ϕ(-q)v(q) has the minimum energy configuration

ϕ(q) = v(q) E(q) (43)
with E(q) = q 2 for a single layer. Thus choosing a potential

v(q) = 6iẑ • (q × b) (44)
we recover the desired far-field configuration.

To examine the effect of many pairs j with locations r j , z j and separations b j we impose on the multilayer system the potential

v tot (q) = 6i j ẑ • (q × b j ) e -i(q ⊥ rj +qzzj ) . ( 45
)
The ground state in the presence of these pairs is again given by ( 43), but now with the multilayer form for E(q). We compute the mean square amplitude of the fluctuations these pairs generate, averaged over bound pair positions with a Poisson distribution at a density ρ, obtaining

[ϕ z (r)] 2 = ρ (2π) 3 d 3 q |v(q)| 2 E 2 (q) (46)
where . . . indicates an average over pair separations b. This integral is convergent at small q for the unfrustrated stacking but divergent for the frustrated systems. Moreover, corrections to a Poisson distribution arising from correlations between pairs appear only at higher order in ρ. Vortex-antivortex pairs in the absence of locking interactions therefore destroy long-range order in the frustrated systems.

We can estimate the correlation length in this disordered state by determining the small-wavevector cut-off for which [ϕ z (r)] 2 ∼ 1. We write |b| 2 ∼ 2 , where is the cut-off scale at which the system reaches the strongcoupling regime with κ ⊥ ∼ 1. This scale is ∼ (βJ ⊥ ) -2 . Then for the abc stacking the correlations lengths in the in-plane and z-directions are

ξ ⊥ ∼ κ -1 ⊥ ( 2 ρ) -2 and ξ z ∼ ( 2 ρ) -1 . ( 47
)
For the abab stacking the corresponding expressions are

ξ ⊥ ∼ κ -1 ⊥ ( 2 ρ) -1/2 and ξ z ∼ ( 2 ρ) -1/2 . ( 48
)
The phase transition to a long-range ordered state involves a competition between this disordering effect of bound vortex pairs, and the opposite tendency produced by locking interactions. A simple estimate for the location of the phase boundary is obtained demanding that the locking interaction at the scale , integrated over the correlation volume, is of order unity.

The most RG-relevant locking interaction for the abc stacking is κ 3 [see Eq. ( 31)]. As this is a coupling between layers three apart, it is not present in the bare description of a system with only nearest-neighbour interactions. It is however generated under the first steps of RG, so that the initial value can be taken to be κ 3,0 ∼ (βJ ⊥ ) 7 (see Sec. VI). At the scale the locking interaction is hence κ 3 ∼ (βJ ⊥ ) 4 . Note that an important role is played by the fact that κ 3 is generated only at high order: if instead one had κ 3,0 ∼ (βJ ⊥ ) 3 as might naively have been expected for a third-neighbour coupling, then the value of κ 3 at scale would be O(1) and independent of J ⊥ . This would leave no scope for a regime with strong interlayer correlations but no long-range order.

For the abab stacking, we have not found locking interactions that are RG-relevant. The leading (least irrelevant) locking term in this case is κ b , given in Eq. (30). At the scale it is of order βJ ⊥ -1/2 ∼ (βJ ⊥ ) 2 .

D. Phase diagram

Combining results from our discussion of RG for weakly coupled layers with our results on the effect of defects in strongly coupled layers, we can determine regimes of behaviour and phase boundaries for systems with each type of stacking, in the limit J ⊥ J. The phase boundaries determined theoretically in this section are compared with Monte Carlo results in Sec. VII.

For the unfrustrated stacking, bound vortex pairs have no important effects. The phase boundary is the point at which y ∼ κ 3 ∼ 1. From the results of Sec. V A, this implies ∼ e 4βJ and βJ ⊥ e 6βJ ∼ 1. Solving approximately in the limit J ⊥ J, the phase boundary is at J ⊥ ≈ Je -6βJ . Interlayer correlations are weak for J ⊥ Je -6βJ while the system has long-range order for J ⊥ Je -6βJ . Within the minimal model of Eq. ( 20), the set of ordered states has a U(1) symmetry, as displayed in Eq. ( 21). This is broken by the (RG-irrelevant) interaction Ṽ (h), introduced for a single layer in Eq. ( 14). It selects integer values of the height field, corresponding to six possible types of three-sublattice spin order.

In contrast, for both types of frustrated stacking, the condition y ∼ κ ⊥ ∼ 1 implies J ⊥ ≈ Je -2βJ . Interlayer correlations in this case are weak for J ⊥ Je -2βJ . The paramagnetic regime with only weak interlayer correlations therefore extends to parametrically lower temperatures and larger values of J ⊥ in these systems than in the unfrustrated stacking. Moreover, because of the effect of bound vortex pairs in systems with frustrated stacking, long range order appears at still lower temperatures or larger values of J ⊥ than strong interlayer correlations.

In the case of the abc stacking, if long range order is stabilised by generation of the RG-relevant third-neighbour coupling κ 3 the condition κ 3 ξ 2 ⊥ ξ z ∼ 1 implies order for J ⊥ Je -5βJ/3 . Alternatively, order may be stabilised by residual contributions from the RG-irrelevant coupling κ b . Specifically, RG flow stops on the scale at which κ ⊥ ∼ 1. At this scale, interactions (whether RG-relevant or RG-irrelevant) that break the U (1) × U (1) symmetry of Eq. ( 26) down to a discrete one will act coherently over a correlation volume. This ordering tendency competes with the disordering effect of bound vortex-antivortex pairs. Since κ 3 is generated rather slowly under RG, RG-irrelevant interactions turn out to be the dominant cause of locking if microscopic interactions are just nearest neighbour. 29 The condition κ b ξ 2 ⊥ ξ z ∼ 1 implies order for J ⊥ Je -20βJ/11 .

For the abab stacking, locking is driven only by irrelevant interactions. Taking into account the dependence of ξ ⊥ and ξ on ρ for the abab stacking, the condition (βJ ⊥ ) 2 ξ 2 ⊥ ξ z ∼ 1 yields a boundary for long range order at J ⊥ ≈ Je -5βJ/3 . In summary, with J ⊥ J, the classical spin liquid regime, in which correlations are strong both within and between layers, extends for both types of frustrated stacking over the interval e -2βJ J ⊥ /J e -cβJ (49) with c = 20/11 for the abc stacking and c = 5/3 for the abab stacking.

E. Spin correlations from the height model

In the classical spin liquid regime, in which interlayer correlations are strong but there is no long-range order, the system is approximately ordered within each correlation volume ξ 2 ⊥ ξ z but different correlation volumes are essentially independent. We can compute correlations approximately in this regime as an average over all ground states. The starting point for this calculation is the expression (13) for spin variables in terms of height fields, and the expressions ( 28) and ( 29) for ground states of the minimal height models in the systems with frustrated stackings.

We require Fourier components of the spin density at wavevectors that are close in-plane to either of the corners K and K of the triangular-lattice Brillouin zone. To obtain the leading contribution at long distance it is sufficient to use the approximation σ j,z ∼ cos π 3 (h z (r j ) + s α ), omitting higher harmonics in h z (r j ).

Recalling that s α = 0, ±2 on the three sublattices, we have for the aaa stacking e iK•rj,z = e iπsα/3 and e iK •rj,z = e -iπsα/3 . The same result holds for the abab stacking on one of the two layers in the unit cell, but for the abc stacking it is necessary to take account of the relative displacement e 1 of neighbouring sites on the same sublattice in successive layers. We have (modulo 2π)

(K + n 1 A 1 + n 2 A 2 ) • r j,z = π 3 s α -z(K + n 1 A 1 + n 2 A 2 ) • e 1 = π 3 (2pz + s α )
with p = n 1 + n 2 , and

(K + n 1 A 1 + n 2 A 2 ) • r j,z = π 3 (2p z -s α ) with p = 2 + n 1 + n 2 .
Retaining only smoothly varying contributions, we can then write for q ⊥ small but q z arbitrary j σ j,z e i(K+q)•rj,z ∼ d 2 r e -i π 3 hz(r) e i(q ⊥ •r+[qz+ 2π

3 p]z)

and j σ jz e i(K +q)•rjz ∼ d 2 r e i π 3 hz(r) e i(q ⊥ •r+[qz 2π 3 p ]z)
where we can include the aaa and the a-layers of the abab stacking by setting p = p = 0 in these cases.

We use these expressions to evaluate

S(K + q) = j,z σ 0,0 σ j,z e i(K+q)•rj,z (50)
and the equivalent with K in place of K, computing the average . . . over ground states [Eqns. (28) and ( 29)].

For the abc stacking this gives

S(K + q) ∝ δ(q x - π 3 κ ⊥ cos[q z + 2π 3 p]) × δ(q y + π 3 κ ⊥ sin[q z + 2π 3 p]) (51)
and

S(K + q) ∝ δ(q x + π 3 κ ⊥ cos[q z + 2π 3 p ]) × δ(q y + π 3 κ ⊥ sin[q z + 2π 3 p ]) . (52)
For the abab stacking, following our discussion in Sec. III B 2, we focus on the contribution to the structure factor from sites on only one of the two sublattices by restricting z,µ to the layer µ = 1. This gives

S(K + q) = S(K + q) ∝ δ(q z )δ (2) (q 2 ⊥ -[ π 3 ] 2 κ 2 ⊥ ) . ( 53
)
It is reasonable to expect that the main consequence of finite correlation lengths ξ ⊥ and ξ z will be broadening of the delta functions in these expressions for S(q). Making that allowance, we see that the height model calculation produces results similar to the ones from the SCGA and from Monte Carlo simulations.

VI. RENORMALISATION GROUP FLOWS BEYOND LEADING ORDER

Our calculation of RG flow is perturbative in interlayer coupling and vortex fugacity. We can improve the estimates of the previous section by including terms to higher order. Qualitatively, this has two potentially important consequences. First, the in-plane stiffness K becomes scale-dependent and interlayer gradient couplings are generated under the RG flow. This in turn modifies the dimensions of the various operators discussed above. Second, for the abc stacking, the relevant furtherneighbour couplings that break the U(1) symmetry under spatial rotations are generated from the irrelevant contribution to the nearest-neighbour interlayer coupling, Eq. (30).

A. Simply stacked triangular layers

To set the stage, it is instructive to consider the case of aaa-stacked triangular layers. The model [Eq. (20)] is simply a 3D XY model, in which the coupling between neighbouring layers is much weaker than the intralayer coupling. For small interlayer couplings there is a regime where the RG flows are well-described by those of a system of coupled 2D XY models 35,36 . Though this treatment is not adequate to describe the transition between the low-temperature ordered phase and the hightemperature paramagnet, which is in the 3D XY universality class, it represents behaviour well so long as the renormalised interlayer coupling is not strong.

For uncoupled layers, two different ways exist to derive RG equations. The original work by Kosterlitz and Thouless 30,31 on the 2D XY model used a real-space calculation, integrating out vortex-antivortex pairs separated by less than a minimum length scale , and this method has been extended to include models analogous to (20) with vortices [START_REF] Minnhagen | [END_REF] . Somewhat later, the momentumshell RG approach was applied to these systems 34 and we use this second approach, which is more transparent in the case of the frustrated abc and abab stackings. We review the method and give technical details of our calculations in Appendix C; here we discuss the physical implications of the results.

Including the most relevant interlayer couplings, the marginal gradient couplings introduced in Eq. ( 31), and a new second-layer coupling term cos π 3 (h z+2 (r) -h z (r)) with coefficient g 2 , the RG equations additional to (33) to quadratic order in κ 3 and y are

∂K ∂ ln = c 1 κ 2 3 -y 2 K 2 ∂K 1 ∂ ln = -c 1 κ 2 3 ∂g 2 ∂ ln = g 2 (2 - π 18K ) -c 2 κ 3 3 . (54)
Here we have allowed for the effect of fluctuating bound vortex pairs on the stiffness. A deficiency of the momentum-space approach is that this correction cannot be evaluated easily, and so we take the result computed in the real-space RG using Coulomb gas methods [START_REF] Nienhuis | Phase Transitions and Critical Phenomena[END_REF] . The constants c 1 and c 2 are given in Eq. (C37).

The RG flow described by Eqs. ( 33) and (54) includes several important effects. First, at this order the stiffness K flows towards smaller values if vortices dominate. As the interlayer coupling κ 3 is irrelevant if K is sufficiently small, this ensures that the paramagnetic phase is stable to weak interlayer coupling. Second, new interlayer couplings are generated from κ 3 : the marginal gradient coupling K 1 and the relevant second-neighbour coupling g 2 . The latter contributes to stabilising long-range order if vortices are not dominant.

Interlayer gradient couplings change the scaling dimensions of other interlayer couplings and of the fugacity for multilayer complexes of vortices. The scaling dimensions of Eq. ( 34) become more generally

β 1 = π 18 π -π dk z 2π 1 -cos k z K 0 + p K p cos pk z α 1 = 9 π i,j σ i σ j K |i-j| , (55)
where σ i is the vortex strength in layer i.

A striking consequence of interlayer gradient couplings that follows from these results for scaling dimensions is the possibility of a sliding phase, 37 in which for appropriate values of {K p } neither vortices nor interlayer cosine couplings are relevant. The window of stability of this phase is however quite narrow, and it does not seem likely that it would be reached by RG flow starting from stacked TLIAFMs with only nearest-neighbour interactions, whether frustrated or not.

B. abc stacking

We now consider the abc stacking. As for the aaa stacking, under RG at second order the stiffness K flows and further-neighbour interactions are generated. The most important of these are shown in Eq. ( 31) with coupling constants denoted by κ 2 and κ 3 . As they break the spatial U(1) symmetry of H (abc) [see Eq (24)], their generation involves the RG-irrelevant nearest-neighbour interaction κ b appearing in Eq. (30). The coupled RG equations

∂K ∂ ln = c 3 κ 2 ⊥ -y 2 K 2 ∂K 1 ∂ ln = -c 4 κ 2 ⊥ ∂κ b ∂ ln = - π 18K κ b ∂κ 2 ∂ ln = κ 2 1 - π 18K + c 5 κ ⊥ κ b ∂κ 3 ∂ ln = κ 3 2 - π 18K + c 6 κ ⊥ κ 2 (56)
and values of the constants c 3 , c 4 , c 4 and c 6 are given in Eq. (C37); both c 5 and c 6 are proportional to K 1 for small K 1 . For each coupling, we have included the flow due to its scaling dimension, as well as (for those not initially present in the nearest-neighbour model) the leading-order term that generates it. Flow of the vortex fugacity y is given in Eq. (33). For the in-plane stiffness K, we have included the leading-order non-vanishing contributions to its RG flow, demonstrating that this is slow.

The most important physical effect captured by this second-order calculation is the generation of the locking interaction κ 3 from κ b (which appears microscopically in a nearest-neighbour model) via the coupling κ 2 . Since κ 3 is more strongly RG-relevant than κ 2 (which has the same scaling dimension as κ ⊥ ), it is the key interaction. It is generated only in the presence of non-zero K 1 , itself produced from the nearest-neighbour interaction κ ⊥ . Combining these steps, we find for a system with initial values κ ⊥ = κ ⊥,0 , κ b ∼ κ ⊥,0 and K 1 = κ 2 = κ 3 = 0, that κ 3 ∼ (κ ⊥,0 ) 7 is generated after an RG scale change of order one. As discussed in Sec. V D, this locking interaction stabilises long-range order if it dominates over the disordering effects of vortex-antivortex pairs.

We have not examined RG for the abab stacking in detail beyond leading order, since we have not identi-fied RG-relevant interactions that break the continuous ground-state symmetry of the minimal model. Symmetry is instead broken by RG-irrelevant nearest-neighbour interactions that are present microscopically, as discussed in Sec. V D.

VII. DISCUSSION

The results from the three approaches we have presented -the self-consistent Gaussian approximation, Monte Carlo simulations, and analysis of height modelsestablish a consistent picture. They show that triangular lattice Ising antiferromagnets with frustrated stackings exhibit classical spin liquid behaviour over an extended temperature range if interlayer coupling is weak. In this regime, there are strong correlations within and between layers, but without long-range order.

The most significant weakness of the SCGA is that it fails to capture the ordering transition, giving instead a finite correlation length at all non-zero temperatures. The SCGA also predicts a temperatureindependent value for the helix radius Q, while within the height model Q is a function of βJ ⊥ . Small increases in Q with decreasing T at fixed J ⊥ are apparent in Fig. 7(b), although the anticipated continuum behaviour is not fully-developed.

Some more detailed comparisons between Monte Carlo simulations and height model calculations are possible. The prediction of Sec. V D that the ordering transition is at larger values of J ⊥ and smaller temperatures in systems with frustrated stacking compared to the unfrustrated case (J ⊥ ≈ Je -20βJ/11 or J ⊥ ≈ Je -5βJ/3 compared with J ⊥ ≈ Je -6βJ ) is clearly consistent with simulation results shown in Fig. 4. For a quantitative test, we fit the phase boundaries determined in simulations to the form J ⊥ = AJ exp(-cβJ). We obtain c = 1.90 ± 0.08 for the abc stacking, c = 1.63±0.11 for the abab stacking, and c = 5.44 ± 0.2 for the unfrustrated case, in striking agreement with analytical results. Values of the other fitting parameter are A = 2.87 ± 0.2, A = 2.16 ± 0.27 and A = 6.43 ± 0.5, respectively.

Appendix B: Analysis of Monte Carlo results

In this appendix we discuss in further detail our Monte Carlo results for S(q) and the fitting procedures used to analyse them.

As a simple check, we start by considering uncoupled layers, which are expected to display power-law correlations at low temperature with S(K + q) ∝ q -3/2 . The behaviour illustrated in Fig. 16 matches this quite accurately. Interlayer interactions produce significant changes in S(q), and no clear remnant of the 3/2 power law is identifiable even for the smallest values of J ⊥ /J that we have investigated. Instead, we find for non-zero J ⊥ that S(q) is well-represented using Lorenztian functions of wavevector. FIG. 16. Illustration of power-law behavior without interlayer coupling: line has slope -3 2 ; system parameters are L = 72 and T = 0.31J.

Correlations for the abc stacking

The data displayed in Fig. 6 show helices of high intensity with axes passing through the K-points of the triangular-lattice Brillouin zone. In broad terms, we extract the correlation length ξ ⊥ and the helix radius Q by analysing simulation results for S(q) separately at each q z , and fitting data near the maximum to a sum of Lorentzian contributions, one from for each helix that intersects the plane.

In detail, we consider values of S(q) at fixed q z with (q x , q y ) spanning one Brillouin zone. To focus on the maxima, we retain the N largest values of S(q) from a total of L 2 points within each q z -plane. If N is too large, some points are included that are too far in reciprocal space from the helix to be well-represented by the fitting function; if N is too small, statistical accuracy is sacrificed. Results are insensitive to the choice of N in the range 20 ≤ N ≤ 200, and we use N = 50. Referring to Fig. 17, the form of S(q) near the K-points labelled a and b should be dominated by helices with their axes passing through these K-points, but may also be influenced by helices with axes passing through the four K-points c - f if the helix radius is large. Our fitting function

F 4nn (q ⊥ ) = i I ξ 2 ⊥ (q ⊥ -q ⊥,i ) 2 + 1 (B1)
therefore includes six terms, labelled by i. Since the different values of q ⊥,i are related by symmetry, it contains four real scalar fitting parameters. The quality of fit we obtain in this way is illustrated in Fig. 18.

In principle, one expects S(q) to be characterised by two distinct correlation lengths, ξ ⊥ and ξ z , as discussed in Sec. V. In practice, we have been unable to extract a second correlation length from our Monte Carlo data for the abc stacking, for reasons we now discuss. Consider first the ideal form of correlations, reached in the limit of divergent correlation lengths: S ideal (q) = δ q x -q 0

x (q z ) δ q y -q 0 y (q z ) .

(B2)

The consequences of finite correlation lengths can be represented by convolving S ideal (q) with a form factor that is characterised by its width in two directions transverse to the line q 0 x (q z ), q 0 y (q z ). The fitting function F 4nn (q ⊥ ) corresponds to a choice for this form factor that has circular contours in the q x -q y plane. More general possibilities have elliptical contours; we have made fits of this type, Fixed q ⊥ , data Fixed q ⊥ , fit FIG. 19. S(q) vs qz for fixed qx, qy in the abc stacking, comparing data and fitting function. J ⊥ = 0.1J, L = 36, Lz = 48, T = 0.56J. q z FIG. 20. S(q) vs qz, for fixed qx, qy passing through the maximum, in the abab stacking: data (red); fit to SCGA (green); sum of Lorentzians (black). J ⊥ = 0.20J, T = 0.73J.

but find they do not show significant in-plane anisotropy. As a demonstration that the form F 4nn (q ⊥ ) is an adequate representation of our data, we show in Fig. 19 a comparison of it with Monte Carlo data, as a function of q z at fixed q x , q y , on a line passing through the helix.

The close match indicates that the broadening within the q x -q y plane that is contained in F 4nn (q ⊥ ) also accounts for the broadening of the helix along q z .

Correlations for the abab stacking

For the abab stacking, our fitting of S(q) as a function of q x and q y follows similar steps to the ones used for the abc stacking, but analysis of the dependence on q z has new features. For this stacking the peak width of S(q) as a function of q z yields directly the interlayer correlation length ξ z . An example of a fit is shown in Fig. 20 and the resulting values of ξ z are displayed as a function of J ⊥ and T in Fig. 21 The unit of length is the spacing between successive a-layers.

Appendix C: RG calculations

Here we present technical aspects of our RG calculations, following a standard momentum-shell approach 38 . The general method is as follows. Our objective is to evaluate correlation functions or the partition function

Z = D[h]e -(H0+H1) (C1)
with an initial momentum cutoff Λ = 1/ , where is the lattice constant. Here H 0 is a quadratic effective Hamiltonian, which may include both in-plane and interplane gradient terms:

H 0 = 1 2 z d 2 r K (∇h z (r)) 2 + p>0 K p ∇h z (r) • ∇h z+p (r) . (C2)
We divide the height field into short-wavelength and long-wavelength modes by writing

h > z (r) = Λ/s<|q|<Λ d 2 q h z (q)e iq•r and h < z (r) = |q|≤Λ/s d 2 q h z (q)e iq•r . (C3)
A new effective Hamiltonian H eff with a reduced cutoff Λ/s is obtained by integrating out the short-wavelength modes, and then re-scaling all in-plane lengths by s. Note that we retain the layer index z as a discrete variable, and coarse-grain only the in-plane co-ordinates. Expanding in powers of

H 1 e -H eff = z D[h > z ]e -(H0+H1) (C4) ≈ z D[h > z ]e -H0 1 -H 1 + 1 2 H 2 1 + ... .
To quadratic order, the effective Hamiltonian with the reduced cutoff Λ/s is

H eff = H 0 + H 1 0 - 1 2 H 2 1 0 + ... (C5)
where all terms are functions of only the long-wavelength fields h < z (r), the average . . . 0 is over short-wavelength fields with weight e -H0 and H 0 is obtained from H 0 by omitting the short-wavelength fields. As a final step, lengths in H 0 , and in the expectation values on the right are re-scaled according to r → sr.

First-order calculation

We derive the first-order RG equations as follows. Consider the interlayer coupling

H 1 = κ ⊥ d 2 r [∂ x h z (r) cos δh z (r) -∂ y h z (r) sin δh z (r)] = κ ⊥ Im d 2 r ∂ ζ z e iδhz(r) (C6)
where we have introduced ζ z = x z + iy z , and z denotes the layer with which the coordinates x, y are associated.

We have

H 1 0 = κ ⊥ Im d 2 r ∂ ζz e iδh<(r,z) e iδh>(r,z) 0 . (C7) Defining F n ≡ h > (r, z)h > (r, z + n) = γ n 4π 2 K Λ>|q|>Λ/s d 2 q 1 q 2 = γ n 2πK log s (C8) with γ n = 1 2π 2π 0 dk z K cos nk z K + p>0 K p cos pk z and β n = π 18K (γ 0 -γ n ) , (C9)
we obtain

e iδnh(r,z) = exp - π 2 18 δ n h(r, z) 2 = exp - π 2 9 (F 0 -F n ) = s -βn . (C10)
The re-scaling r → sr gives

H eff = κ ⊥ d 2 r Im ∂ ζz e iδh<(r,z) s 1-β1 (
C11)

In the continuum limit s → 1 we have

∂κ ⊥ ∂ ln = (1 -β 1 )κ ⊥ . (C12)
Scaling dimensions of the other operators can be deduced in a similar way.

Second-order calculation

At second order, we must evaluate the quadratic terms in Eq. (C5). It is useful to introduce some notation. Let H n (z) denote a contribution to interlayer coupling involving the height differences δ p h z (r) ≡ π 3 [h z+p (r)h z (r)] and define

∆ m,n,z-z = H m (z)H n (z ) 0 -H m (z) 0 H n (z ) 0 .
We are primarily interested in two types of such term: those that contribute to the most relevant interlayer couplings, and those that contribute corrections to the marginal gradient couplings.

a. Corrections to gradient couplings

We first compute corrections to the gradient couplings that are generated by ∆ n,n,0 for various n. An example is

∆ 3,3,0 = H 3 (z)H 3 (z) 0 -H 3 (z) 0 H 3 (z) 0 = (κ 3 ) 2 2 d 2 r d 2 r cos δ 3 h < z (r) + δ 3 h < z (r ) e iδ3h > z (r) e iδ3h > z (r ) 0 -s -2β3 + cos (δ 3 h < (r, z) -δ 3 h < (r , z)) e iδ3h>(r,z) e -iδ3h>(r ,z) -s -2β3 (C13)
We write

e iδ3h > z (r) e iδ3h > z (r ) 0 -s -2β3 = s -2β3 e -4πβ3G(r-r ) -1 , (C14)
where

G(r) = Λ/s<|q|<Λ d 2 q 4π 2 e iq•r q 2 . (C15)
Assuming that e 4πβ3G(r-r ) -1 is small unless |r -r | 1, we expand the long-wavelength height fields in R = rr to obtain

∆ 3,3,0 = (κ 3 ) 2 2 d 2 r d 2 r cos 2δ 3 h < z (r) + (r -r ) • ∇δ 3 h < z (r) sin 2δ 3 h < z (r) + ... s -2β3 e -4πβ3G(r-r ) -1 + 1 2 d 2 r d 2 r 1 - 1 2 (r -r ) • ∇δ 3 h < z (r) 2 s -2β3 e 4πβ3G(r-r ) -1 (C16)
The terms in the first line are new, less relevant couplings between spins 3 layers apart, and can be ignored. The first term in the second line is a constant, and the second term in the second line is the contribution to the gradient energy that we are interested in. After performing the angular integration, only terms of the form (∇δh < ) 2 remain and we obtain

∆ 3,3,0 = - (κ 3 s -β3 ) 2 8 B d 2 r |∇(δ 3 h < z (r))| 2 + const + ..., (C17)
where the ellipsis represents the less relevant interlayer couplings, and

B = d 2 R R 2 e 4πβ3G(R) -1 . (C18)
The other important corrections to the gradient energy come from ∆ 1,1,0 and ∆ 2,2,0 . These contain two types of terms, with the forms

∆ 1,1,0 (+) = -(κ ⊥ ) 2 d 2 r d 2 r ∂ ζ z ∂ ζ z e i(δhz(r)+δhz(r )) 0 -∂ ζ z e iδhz(r) 0 ∂ ζ z e iδhz(r ) 0 + c.c. and ∆ 1,1,0 (-) = (κ ⊥ ) 2 d 2 r d 2 r ∂ ζ z ∂ ζ z e i(δhz(r)-δhz(r )) 0 -∂ ζ z e iδhz(r) 0 ∂ ζ z e -iδhz(r ) 0 + c.c. . (C19)
Terms of the first type generate new (but irrelevant) inter-layer couplings that do not lift the helical degeneracy; they are not important for our analysis. We are interested in terms of the second type, which reduce to

∆ 1,1,0 (-) = (κ ⊥ ) 2 2 d 2 r d 2 r ∂ ζ z ∂ ζ z cos δh < z (r) -δh < z (r ) s -2β1 e 4πβ1G(r-r ) -1 (C20)
Differentiating both slow and fast fields, and expanding for small R, we obtain the four terms

∆ 1,1,0 (-) = (κ ⊥ ) 2 2 s -2β1 d 2 r C 0 |∇h < z (r)| 2 -C 1 ∇h < z (r) • ∇δh < z (r) + C 2 |∇δh < z (r)| 2 + C 3 + irrel. (C21)
where

C 0 = d 2 R e 4πβ1G(R) -1 , C 1 = 4πβ 1 d 2 R (R • ∇ R G(R)) e 4πβ1G(R) = -C 0 C 2 = - 1 4 d 2 R R 2 ∇ 2 R e 4πβ1G(R) -8π 2 β 2 1 |∇ R G(R)| 2 e 4πβ1G(R) = - 1 4 C 0 + π 2 9K (γ 1 + γ 0 ) d 2 R R 2 ∇ 2 R G(R)e 4πβ1G(R) C 3 = d 2 R π 2 9K γ 0 ∇ 2 R G(R) + 2πβ 2 1 |∇ R G(R)| 2 e 4πβ1G(R) (C22)
and we have exploited symmetries in the integration over R. Summing over layers, the contribution to the gradient energy is

δH = (κ ⊥ s -β1 ) 2 2 z d 2 r 2C 2 |∇h < z (r)| 2 -(2C 2 -C 0 ) ∇h < z (r) • ∇δh < z+1 (r) . (C23)
A similar contribution arises from ∆ 2,2,0 . Although the leading irrelevant terms in the interlayer coupling also renormalise the gradient energy, we will neglect their effect here as it influences only the initial part of the RG flow.

b. Generation of new interlayer couplings

We now turn to the most important part of our RG calculation, which is to determine at what order in κ ⊥ the relevant inter-layer coupling H 3 is generated in a microscopic theory with only nearest-layer couplings. (Recall from Eq. ( 31) that H 3 couples layers three apart in the abc stacking.) Importantly, we show that though one might expect H 3 to be produced at order κ 3 ⊥ , in fact this is not the case: generating this interaction requires K 1 = 0, and it consequently appears at order κ 7

⊥ . If H 3 is absent, then to generate it we must keep the leading irrelevant term that breaks the U(1) symmetry, which is the interlayer coupling H b (z) of Eq. (30). Then H 2 [Eq. (31)] is generated by the bilinear

∆ b,1,1 = H 1 (z)H b (z + 1) -H 1 (z) H b (z + 1) (C24)
and H 3 is generated by

∆ 2,1,1 = H 1 (z)H 2 (z + 1) -H 1 (z) H 2 (z + 1) . (C25)
Other cross-terms, such as H 1 (z)H 1 (z + 1) and H 3 (z)H 1 (z) also generate new inter-layer couplings. However, for our purposes these can be ignored: they are either less relevant than the terms listed above, or equally relevant but appear at a higher order in κ ⊥ . We have

∆ 2,1,1 = -κ ⊥ κ 2 d 2 r d 2 r ∂ ζ z ∂ ζ z+1 e i(δhz(r)+δ2hz+1(r )) -∂ ζ z e iδhz(r) ∂ ζ z+1 e iδ2hz+1(r ) + h.c. + ... ≈ -κ ⊥ s -β1 κ 2 s -β2 C 3 (2) d 2 r cos (δ 3 h z (r)) + less relevant terms (C26)
where +... represents a contribution that generates terms of the form exp i π 3 (h z+3 (r ) + h z (r) -h z+1 (r ) -h z+1 (r)) , which produce inter-layer couplings less relevant than the terms of interest, which have been neglected in the second line. Additionally, in the second line we have kept only terms in which all derivatives are applied to the fast height fields, as these generate the most relevant inter-layer coupling, and neglected all but the leading order term in a derivative expansion of the argument of the cosine term for small R = rr . The coefficient is where again, +... generates interlayer couplings of the form exp i π 3 (h z+2 + h z -2h z+1 ) , which we omit as they are less relevant. In this case, because H 2 involves one derivative of the slow height fields, we must calculate two terms. First, applying all three of the derivatives to h > gives the terms Next, we Taylor expand for small R = rr . After integrating over R, the leading-order term vanishes, and the most relevant term that we are left with is

C 3 (n) = π 2 9K d 2 R -γ 1 ∇ 2 R G(R) + π 2 9K (γ 1 -γ n+1 ) (γ 1 -γ 0 ) |∇ R G(R)|
κ ⊥ s -β1 κ b s -β1 C 4 d 2 r ∂ x δ 2 h < z (r) cos δ 2 h < z (r) + ∂ y δ 2 h < z (r) sin δ 2 h < z (r) (C28)
where with C 3 (n) as defined in Eq. (C27). The remaining contributions, which come from applying two or three derivatives to h < , necessarily produce irrelevant couplings, and are safely omitted from our RG calculation. In total, we therefore obtain

C 4 = π 2 9K 3 (γ 1 -γ n+1 )(γ 1 -γ 0 ) 2 d 2 R R x ∂ x G(R)|∇G(R)|
∆ b,1,1 = κ ⊥ s -β1 κ b s -β1 d 2 r C 3 (1) - π 3 C 4 ∂ x h < z (r) cos δ 2 h < z (r) + ∂ y h < z (r) sin δ 2 h < z (r) + π 3 C 4 ∂ x h < z+2 (r) cos δ 2 h < z (r) + ∂ y h < z+2 (r) sin δ 2 h < z (r) + irrelevant terms (C31)
Though the first terms in the second line are just as relevant as the terms in the first line, they do not contribute to generating H 3 until inter-layer kinetic terms K n are generated for n > 1 (see Eq. (C2)). Hence we have neglected them in our discussion, as their impact on the other couplings in the RG is very weak. We will also see presently that the coefficient C 4 is negligibly small compared to C 3 (1). We emphasise that both C 3 (1) and C 4 are of order at least κ 2 ⊥ , since for K 1 K and K n = 0, n > 1, we have

γ 0 = 1 1 -K1 K 2 ≈ 1 , γ 1 = K K 1   1 - 1 1 -K1 K 2   ≈ K 1 2K , γ n = 0 , n > 1 , (C32)
and K 1 is generated only at order κ 2 ⊥ . Therefore, in summary, κ 3 is generated not at order κ 3 ⊥ , as one might naively have expected, but at order κ 7 ⊥ . A similar effect was noted for frustrated couplings in Ref. 39.

Evaluation of coefficients

To proceed further, we must evaluate the coefficients B, C 0 , C 1 , C 2 , and C 3 (n). In order to compute the relevant integrals, we expand the exponentials for small G(r). (We will justify this expansion presently). To ensure that all integrals are absolutely convergent, we take our system to have a finite size L x = L y = L, and use periodic boundary conditions. In this case the first-order terms vanish after integration, since

d 2 r G(r) = 1 L 2 d 2 r Λ/s<q<Λ e iq•r q 2 = 1 L 2 Λ/s<q<Λ 1 q 2 d 2 r e iq•r = Λ/s<q<Λ
1 if q = 0 0 else .

A similar derivation applies for derivatives of a single power of G, which also vanish. The leading-order contributions are therefore quadratic in G. Keeping only these terms, the integrals of interest are The two integrals not involving explicit powers of r are easily evaluated as

I 0 = d 2 r (G(
I 0 = d 2 r (G(r)) 2 = 1 L 4 d 2 r q,k e i(q+k)•r q 2 k 2 = 1 L 2 q 1 q 4 ≈ 1 4π 2 Λ Λ/s d 2 q q 4 = 1 2πΛ 2 s 2 -1 2 ≈ ds 2πΛ 2 I 3 = d 2 r G(r)∇ 2 r G(r) = 1 L 4 d 2 r
q,k e i(q+k)•r q 2 ≈ 1 (2π) 2 Λ Λ/s d 2 q q 2 = log(s) 2π ≈ ds 2π .

We note, somewhat surprisingly, that it is the terms quadratic in G(r) -rather than the linear terms -that are proportional to ds. Our Taylor expansion is nevertheless justified: for higher powers of G, the δ-function constraint takes the form δ( n i=1 k i ). In practice, this means that non-zero contributions to momentum integrals require both that |k i | is within the momentum shell for each i and also that | n-1 i=1 k i | lies in this shell. This leads to a strong phase-space suppression of the relevant angular integrals for n > 2, justifying the quadratic approximation used here.

Evaluating I 1 and I 2 , we encounter a second difficulty: the resulting integrals retain an explicit dependence not only on the cutoff Λ but also on the system size L. This stems from the factors of r 2 in the integrands, which arise from Taylor expansions of the type f (h < (R) -h < (R + r)) ≈ f (r • ∇h < (R)) + ... (C33)

followed by an expansion of f for small r. The expansion is justified if G(r) falls off sufficiently quickly in r that only small values of r contribute; however, in the cases at hand this is not so.

To circumvent this difficulty, we instead expand the function f to quadratic order in the difference h < (R)-h < (R+r) of the height fields, without making a Taylor expansion of the height fields in powers of r. The approach amounts to the substitution d 2 q q 2 h < q h < -q d 2 r r 2 F (r) ≈ d 2 q h < q h < -q d 2 r 2(1 -cos q • r)F (r) .

Using this, we obtain

q 2 I 1 = 2 d 2 r (1 -cos q • r) (G(r)) 2 = 1 L 4 k1,k2 d 2 r 2 e i(k1+k2)•r k 2 1 k 2 2 - e i(k1+k2+q)•r k 2 1 k 2 2 - e i(k1+k2-q)•r k 2 1 k 2 2 ≈ 1 4π 2 d 2 k k 2 2 k 2 - 1 |k + q| 2 - 1 |k -q| 2 , ≈ -q 2 ds πΛ 4 (C34)
where in the last line we have kept terms only to quadratic order in q, as higher-order terms are RG-irrelevant. Similarly, we may evaluate

q 2 I 2 = 2 d 2 r (1 -cos q • r)G(r)∇ 2 r G(r) = 1 L 4 k1,k2 d 2 r 2 e i(k1+k2)•r k 2 1 - e i(k1+k2+q)•r k 2 1 - e i(k1+k2-q)•r k 2 1 ≈ 1 4π 2 d 2 k 2 k 2 - 1 |k + q| 2 - 1 |k -q| 2 ≈ -q 2 ds πΛ 2 . (C35)
Note that in both of these evaluations, we have ignored an important constraint, which is that for the terms involving q, we must have Λ/s ≤ |k + q| ≤ Λ, in addition to Λ/s ≤ k ≤ Λ. However this constraint, if included, will modify the result by a factor of order unity, provided that q is not large compared to Λ -Λ/s. The final results are

B = -16πβ 2 3 ds Λ 4 , C 0 = 1 2 (4πβ 1 ) 2 I 0 = 4πβ 2 1 ds Λ 2 , C 4 = 0 (C36) C 2 = - C 0 4 + π 2 9K (4πβ 1 )(γ 1 + γ 0 )I 2 = -π β 2 1 + 8 π 18K 2 γ 2 0 -γ 2 1 ds Λ 2 C 3 (n) = - π 2 9K 2 [γ 1 (γ 0 -γ 1 + γ n+1 -γ n ) + (γ 1 -γ n+1 )(γ 1 -γ 0 )] I 3 = ds 2π π 2 9K 2 [γ 1 γ n -γ 0 γ n+1 ] .
Here the factors of Λ in each coefficient reflect the total engineering dimension of the couplings involved; these factors can be eliminated by defining appropriate dimensionless couplings. From these expressions, we can extract values for the constants appearing in Eqns.

1 FIG. 2 .

 12 FIG. 2. Schematic phase diagram for stacked triangular lattice Ising antiferromagnets. The full line represents the phase boundary, and the dashed line indicates a smooth crossover.

FIG. 3 .

 3 FIG. 3. Location of surfaces on which eigenvalues of the interaction matrix are constant and close to the minimum, for (top) the abc stacking, and (bottom) the abab stacking, at J ⊥ /J = 0.2.

FIG. 4 .

 4 FIG. 4. Phase boundaries for the unfrustrated (aaa) and frustrated (abc and abab) stackings. Points: data from Monte Carlo simulations. Lines: fits to theory of Sec. V D; see discussion in Sec. VII.

  FIG. 5. Distributions P (E) of energy E for the abc stacking with J ⊥ = 0.1J. Top: temperatures in the range 0.4J ≤ T ≤ 0.65J. Middle: temperatures T = 0.50J, 0.51J and 0.52J close to the transition. System size L = 72, Lz = 12. The distribution closest to the transition is solid red outlined in black, centered around E = -1.011 in the top panel, and is the middle temperature in the middle panel. Its bimodal form indicates a first-order transition. Bottom: finite size effects, illustrated for L = 96 and T = 0.52J.

FIG. 7 .

 7 FIG. 7. (a) Correlation length ξ ⊥ and (b) helix radius Q, as a function of temperature for various values of J ⊥ in the abc stacking. ξ ⊥ is measured in units of lattice spacing, Q in units of inverse lattice spacing. Dashed lines are SCGA predictions for Q from Eq. (6). Results were obtained in a system of size L = 72, Lz = 12. Data for each value of J ⊥ extend to the lowest temperature employed in parallel tempering that was above Tc.

3 FIG. 8 .

 38 FIG.8. Development of Bragg peaks in the ordered phase for the abc stacking. Top: Savg(qz) [Eq.(11)] as a function of qz at four selected temperatures near the transition, in a system with J ⊥ = 0.1J. Data are for T = 0.45J, 0.5J, 0.52J and 0.6J, in order of decreasing peak intensity, and the transition temperature is Tc ≈ 0.54J. Bottom: Savg(qz) as a function of T for the three values of qz that are marked with vertical dashed lines in the top panel. Results for both panels were obtained in a system of size L = 36, Lz = 48.

FIG. 10 .

 10 FIG.10. Correlation length, ξ ⊥ , as a function of temperature for various values of J ⊥ in the abab stacking, obtained by fitting to the functional form given in Eq.(12).

3 FIG. 11 .

 311 FIG.11. Comparison of SCGA and simulation results for S(q) in the abab stacking. J ⊥ = 0.4J, T = 1.36J

FIG. 12 .

 12 FIG. 12.Correlation length ξ ⊥ as a function of temperature for various values of J ⊥ as obtained from the SCGA. Top: abc stacking. Bottom: abab stacking. The data labeled 'SCGA' have been derived by imposing the condition |σ| 2 = 1 while the variable λ is used as a fitting parameter in the other curves. For clarity, results from the first of these approaches are shown only at one value of J ⊥ ; agreement is similar at other values of J ⊥ .

FIG. 13 .

 13 FIG. 13. Mapping from Ising spins to heights on the triangular lattice. (a) The height field decreases by 1 (increases by 2) along an unfrustrated (frustrated) bond as an upward-facing triangle is traversed in the counter-clockwise direction. This ensures that the net change in height field around each triangle is zero provided the triangle is in one of its ground states. (b) and (c): Sample patterns of frustrated bonds and height fields. Green (red) edges on the triangular lattice represent frustrated bonds between pairs of up (down) spins; blue edges correspond to unfrustrated bonds. The number at the centre of each triangle indicates the value of the corresponding height variable; the different shades highlight regions with different heights. (b) shows a maximally tilted configuration (height variables at triangle centres decrease as rapidly as possible from left to right), corresponding to the true ground state for the abc and abab stackings; (c) shows a flat, three-sublattice ordered configuration with a single domain wall (height variables differ only along the domain wall). Our convention for the three sublattices of Table I is indicated by the coloured circles: A= solid blue; B= yellow with dashed border; C= open white.

FIG. 14 .

 14 FIG. 14. Energetically preferred domain wall stacking. Arrows at sites of a triangular lattice represent the spin configuration in one layer. The height in this layer increases by 1 moving from the blue region to the white region. The dashed parallel green and red lines indicate the energetically favourable domain walls in a neighbouring layer, with spin orientations as illustrated. The height difference between adjacent layers determines the orientation of the domain walls.

FIG. 15 .

 15 FIG.15. Two stacked layers, with sublattice labels and definitions of the vectors e1, e2, and e3.

FIG. 17 .FIG. 18 .

 1718 FIG. 17. Brillouin zone for the triangular lattice, with Kpoints labelled a-f .

  .

FIG. 21 .

 21 FIG.21. ξz vs T for different values of J ⊥ in the abab stacking. The unit of length is the spacing between successive a-layers.

2 e π 2 9K

 2 (γ0-γ1+γn+1-γn)G(R) . (C27)The frustrated second-layer coupling is generated by∆ b,1,1 = -4κ ⊥ κ b d 2 r d 2 r ∂ ζ z e iδhz(r) (∂ ζ z+1 e iδhz+1(r )/2 ) 2 -∂ ζ z e iδhz(r) (∂ ζ z+1 e iδhz+1(r )/2 ) 2 + h.c. + ...

κ ⊥ κ b d 2 9K 3 (γ 1 - 2 × d 2 r d 2 r

 23122 rd 2 r e i(δh < z (r)+δh < z+1 (r )) lim r →r ∂ ζ z ∂ ζ z+1 ∂ ζ z+1 e iδh > z (r) e iδh > z+1 (r )/2 e iδh > z+1 (r )/2 + h.c. = κ ⊥ s -β1 κ b s -β1 π 2 γ n+1 )(γ 1 -γ 0 ) e i(δh < z (r)+δh < z+1 (r )) ∂ ζ z G(rr )(∂ ζ z G(rr )) 2 + h.c. e -π 2 9K (γ0+γ2-2γ1)G(R) .

r)) 2 , I 1 = d 2 r r 2 (

 212 G(r)) 2 , I 2 = d 2 r r 2 G(r)∇ 2 G(r) and I 3 = d 2 r G(r)∇ 2 G(r)d 2 r|∇G(r)| 2 .

2 c 3

 23 ds = C 2 c 4 ds = C0 2 -C 2 c 5 ds = -1 2 C 3 (1) c 6 ds = 1 2 C 3 (2) .

  Second, applying one derivative to h < in Eq. (C28) generates the contribution toH 2 4κ ⊥ κ b d 2 r d 2 r e iδh < z+1 (r )/2 ∂ ζ z e iδh < z+2 (r )/2 ∂ ζ z e iδh > z (r) (e iδh > z+1 (r )/2 ∂ ζ z e iδh > z+1 (r )/2 ) + h.c. = κ ⊥ κ b d 2 r d 2 r ∂ ζ z e iδh < z+1 (r ) ∂ ζ z e iδh > z (r) ∂ ζ z e iδh > z+1 (r ) ) + h.c. ≈ κ ⊥ s -β1 κ b s -β1 C 3 (1) d 2 r ∂ x h < z (r) cos(δ 2 h < z (r)) + ∂ y h < z (r) sin(δ 2 h <

	z (r))	(C30)

2 e -π 2 9K (γ0+γ2-2γ1)G(R) (C29)

ACKNOWLEDGEMENTS.

We thank F. H. L Essler, O. A. Starykh and especially P. G. Radaelli for discussions. FJB is supported by NSF-DMR 1352271 and Sloan FG-2015-65927. JTC is supported in part by EPSRC Grants Nos. EP/I032487/1 and EP/N01930X/1. LDCJ is supported by the Okinawa Institute of Science and Technology Graduate University.

Appendix A: Reciprocal-space form of interaction

In this appendix we discuss the reciprocal-space form of the interaction. This is input for SCGA calculations and is illustrated in Fig. 3. Definitions of the lattice vectors, reciprocal lattice vectors, and K, K points are given in Eqns. ( 4) and (5).

The contribution for all three stackings from in-plane couplings is

For the aaa stacking, the interplane interactions contribute J ⊥ (q) = cos q z and the combined minima of J(q) ≡ J 2D (q) + J ⊥ (q) are isolated points in reciprocal space, at ( 4π 3 , 0, π) and ( 2π 3 , 2π √ 3 , π). For the abc stacking, setting ζ = 1 + e iq•a1 + e iq•a2 , we can write the interplane coupling as J ⊥ (q) = J ⊥ (ζe -iq•δ + c.c.)/2. The in-plane coupling can also be expressed in terms of ζ, as J 2D (q) = J(|ζ| 2 -3)/2. The combined interaction can hence be put into the form

From this it is clear that the minima of J(q) lie on the lines ζ = -(J ⊥ /J)e iq•δ . If J ⊥ J, these lines are helixes with axes passing through K-points [Eq. ( 5)] of the triangular-lattice Brillouin zone: for J(k) with

where p = n 1 + n 2 , as in Sec. V E. For k = K + n 1 A 1 + n 2 A 2 + q the line is

where p = 2 + n 1 + n 2 . At larger values of J ⊥ /J, the helix is deformed, acquiring triangular projection in the xy plane, but the degeneracy of the line of minima is not lifted. As the abab stacking has two sites per unit cell, the combined interaction in this case is represented by a matrix

Minima lie on the line q z = 0, |ζ| = J ⊥ /J. For J ⊥ J they form circles around the K-points of the triangularlattice Brillouin zone, as shown in Fig. 3.