Spin-liquid behaviour and the interplay between Pokrovsky-Talapov and Ising criticality in the distorted, triangular-lattice, dipolar Ising antiferromagnet

We study the triangular-lattice Ising model with dipolar interactions, inspired by its realisation in artificial arrays of nanomagnets. We show that a classical spin-liquid forms at intermediate temperatures, and that its behaviour can be tuned by temperature and/or a small lattice distortion between a string Luttinger liquid and a domain-wall-network state. At low temperature there is a transition into a magnetically ordered phase, which can be first-order or continous with a crossover in the critical behaviour between Pokrovsky-Talapov and 2D-Ising universality. When the Pokrovsky-Talapov criticality dominates, the transition is essentially of the Kasteleyn type.


Introduction
Spin liquids, which can be found in both quantum and classical systems, are often defined by the absence of symmetry breaking at low temperature. This raises the question of what is happening at low temperature, and the variety of possible behaviour is comparable to that of the vast array of symmetry-broken phases [1]. One of the oldest and best-understood examples of a classical spin liquid is the Ising model on the triangular lattice with nearest-neighbour interactions. It has been known for many years that this fails to order even at zero temperature [2,3], instead forming a critical state with long-range, algebraic spin correlations [4,5]. The key feature of the spin liquid is that there is robust local ordering associated with the requirement that every triangle has only two equivalent spins, but there are exponentially many configurations that respect this local order, resulting in long-range disorder [2].
This behaviour is not confined to T = 0, but holds to a good approximation throughout the region 0 ≤ T J 1 , where J 1 is the nearest-neighbour coupling constant. For T > 0 the correlations between spins are exponentially rather than algebraically decaying, but the correlation length remains large within the low-temperature region. The whole region 0 ≤ T ≤ J 1 can therefore be considered as a spin-liquid, and weakly-correlated, paramagnetic behaviour only occurs for T > J 1 .
The reason that the nearest-neighbour, triangular-lattice, Ising antiferromagnet (TLIAF) is so well understood is that it can be mapped onto a 1D model of free spinless fermions [2,4,6]. This almost magical transformation converts a strongly-interacting spin problem into a non-interacting fermion problem, thus making possible the calculation of virtually all quantities of interest directly in the thermodynamic limit. The key to this "magic" is that the constraints imposed by the strong interactions between spins map directly to the Pauli exclusion principle, and can therefore be dealt with trivially in the fermionic picture.
The situation becomes more difficult, and more interesting, when additional interactions couple spins beyond nearest neighbour. Further-neighbour interactions tend to stabilise an ordered phase at temperatures below a characteristic, further-neighbour energy scale, J fn [7,8,9,10,11], but this leaves open the possibility of spin-liquid behaviour in the temperature window J fn T J 1 .
The difficulty in analysing further-neighbour models is due to the fact that they cannot be mapped onto free-fermion models, and therefore lack simple analytical solutions. Nevertheless, the fermion picture can still provide useful insights, and in particular one can classify different regions of the spin liquid as being weakly or strongly coupled in a fermionic sense. Weak coupling can be expected for J fn T < J 1 where the free-fermion model is only weakly perturbed, and therefore one expects to find the 2D classical equivalent of a Luttinger liquid. On the other hand, for T ∼ J fn the fermionic model is strongly coupled, making simple predictions more difficult.
There is an essentially infinite number of ways to include further-neighbour interactions, and rather than trying to study all possible combinations of couplings, we concentrate on dipolar coupling between spins, where the interactions fall off with the cube of the separation. Nevertheless, we suggest that the results we obtain are very likely to be qualitatively correct for any system in which the coupling constants are monotonically decreasing with distance (if this condition is not respected there are other possibilities [10,11]). The Hamiltonian we consider is therefore given by, where σ i = ±1 denotes an Ising spin at site i, the sum over (i, j) includes all possible pairs of spins with i = j and r i is the position of the ith spin. Our choice to concentrate on dipolar interactions is not made at random, but is motivated by experiment. In particular, artificial spin systems consisting of arrays of out-of-plane nanomagnets arranged on a triangular lattice are starting to be fabricated, and these realise the dipolar TLIAF to a good approximation [12,13]. While artificial systems have been studied for a number of years [14,15,16,17,18,19,20,21,22,23], recent advances have made it possible to make the nano-dots small enough that they remain thermally active at experimentally viable temperatures [24,25], motivating our study of the equilibrium properties of H dip [Eq. 1].
While the primary experimental motivation comes from artificial spin systems, it is worth pointing out that there are many other realisations of TLIAFs with further-neighbour couplings. Examples include crystals of trapped ions [26], the disordered lattice structure of the spin-orbital liquid candidate material Ba 3 Sb 2 CuO 9 [27] and frustrated Coulomb liquids [28].
The only tuneable parameter in H dip [Eq. 1] is the temperature, and this already leads to subtle behaviour. Nevertheless, we also find it interesting to consider a second tuneable parameter, namely a small lattice disortion associated with squeezing the lattice, and this is parametrised by δ as shown in Fig. 1 (we only consider δ > 0).
Such a distortion is both experimentally accessible, since it is possible to build it into the fabrication procedure of artificial spin systems, and convenient, since it can be used to tune the collective transition temperature relative to the single-nano-dot blocking temperature [13].
At least as important is that there is a good theoretical motivation to study the effect of distorting the lattice. For isotropic systems it has been shown that by carefully choosing the relative strengths of the further-neighbour interactions, it is possible to stabilise an intermediate nematic phase that breaks the 3-fold rotational symmetry of the triangular lattice, but not the Ising symmetry [10,11]. This nematic phase can be characterised by a set of fluctuating strings that wind the system and form a disordered grill-like superstructure, and the density of the strings can be controlled by temperature. At low temperatures the transition into an ordered phase takes place via a Kasteleyn transition, and shows Pokrovsky-Talapov critical behaviour, while at higher temperatures the nematic transitions into a paramagnet via a less-interesting first-order transition (see phase diagram in Ref. [11]). The problem with realising such physics in isotropic systems is that one requires a non-monotonically decreasing interaction strength, with, for example, J 5 > J 4 , and this is difficult to find in nature.
By adding anisotropy of the type parameterised by δ, the symmetry distinction between the nematic and the paramagnet is lost, and therefore there is never a high-temperature phase transition between a paramagnet and a nematic, whatever the form of the further-neighbour interactions. However, the anisotropy drives the appearance of the most interesting features of the nematic, even for monotonically decreasing further-neighbour interactions, in particular the stabilisation at low temperature of a state with a tuneable density of fluctuating strings that shows Pokrovsky-Talapov critical behaviour approaching a Kasteleyn transition into a fluctuationless, low-temperature ordered phase. Since in such a situation the isotropic model is expected to have a direct first-order transition from the disordered to the ordered state, the addition of anisotropy also opens up the possibility that there is an unusual tricritical point, as a first-order transition turns into a Kasteleyn transition.
As a foretaste of the results to come, we show in Fig. 1 a simplified phase diagram of H dip [Eq. 1] as a function of T and δ. One can see that there is a spin-liquid region sandwiched between an ordered stripe phase and a weakly correlated paramagnet. The focus of this article will be on the nature of the spin liquid, as well as on the transition from the spin liquid into the ordered phase, and it can be seen that, as expected, this changes from a first-order to a second-order, essentially Kasteleyn, transition via a tricritical point. The boundary between the spin-liquid and the paramagnet is a crossover and not a phase transition, and a naive guess puts this crossover at T ≈J 1 = (J 1A + J 1B + J 1C )/3, where, J 1A , J 1B and J 1C refer to nearest-neighbour interactions along A, B and C bonds (see Fig. 1 for bond labelling). We provide better ways of determining the boundary between the spin-liquid and paramagnetic regions below, but find that they essentially agree with the simple estimate given byJ 1 .
Our results for the dipolar TLIAF are presented in the main text of the article, since this is the most experimentally relevant form of the interactions. The extensive appendicies discuss related, but simpler models, in which the couplings are short range. This allows important features of the TLIAF to be isolated and studied in more detail than is possible for the dipolar model, since there is both more freedom to separate competing energy scales, and the simpler models are more amenable to analytic calculations and larger scale Monte Carlo simulations.

Methods
We employ two complementary methods to study the dipolar TLIAF, Monte Carlo simulation and mapping onto a model of strings/fermions.

Monte Carlo simulations
Monte Carlo is the standard way to simulate 2D Ising systems, but in the case of the dipolar TLIAF proves difficult to equilibriate. To overcome this problem we use a combination of update methods, including parallel tempering, single-spin-flip updates and worm updates.
Equilibriation difficulties are most acute close to the transition temperature, and are related to a vanishingly small density of defect triangles (triangles with three equivalent spins). In consequence local-update algorithms (e.g. single-spin-flip) have freezing problems.
Our solution is to employ a worm algorithm [29,30,31] in which loops are constructed on the dual honeycomb lattice, taking into account the local interactions [32,33,11,34] (see Appendix B for a discussion of dimer configurations on the honeycomb lattice). The sets of Ising spins within these loops are then flipped with high probability, allowing the system to tunnel between very different configurations. For systems with local interactions (e.g. up to 5th neighbour) the loops of the worm algorithm can be constructed such that detailed balance is automatically obeyed, and therefore the algorithm is rejection free [11]. For dipolar interactions the construction of rejection-free updates is prohibitively time consuming, and instead the algorithm uses effective values of the local coupling constants, J worm 1 , J worm 2 and J worm 3 to guide the loop creation. If these are well chosen, then accepting a flip of all the spins within the loop has a reasonable probability. In practice we found that these parameters have to be carefully tuned so as to target the configurations expected just above the phase transition.
For the isotropic, dipolar TLIAF just above the phase transition, an acceptance probability of about 0.035 was found for the worm updates, which dropped to about 0.004 on crossing the transition, before continuing to decrease. By running a dense set of temperatures, parallel tempering steps were accepted with a probability of at least 0.8 across the transition, providing a considerable aid to equilibriation.
Increasing the transition temperature, for example by adding a lattice distortion, simplifies the simulations by increasing the density of defect triangles in the neighbourhood of the transition. Once this density is high enough, it is possible to use a simple single spin-flip algorithm in combination with parallel tempering.
Simulations are run on hexagonal clusters that preserve all the symmetries of the triangular lattice and have periodic boundary conditions. The linear size of the clusters is L and the total number of triangular-lattice sites is N = 3L 2 . For the dipolar TLIAF we typically use clusters sizes of L = 24, L = 36 and L = 48. While in 2D the dipolar energy is convergent, it was found to be useful to use Ewald summation of the interactions to take into account their slow decay [35].
When performing Monte Carlo simulations a number of physical quantities are sampled. This includes the energy, E, and heat capacity, C, as well as the stripe order parameter and its associated susceptibility, where α ∈ {A, B, C}, i labels triangular lattice sites and τ α i is the spin at the ith site for perfect stripe order parallel to the α bond direction. It is also useful to track the triangular average of the winding number, W = (W 1 , W 2 ) (see below for the definition of the winding number, and also Appendix B), given by, which is designed such that W tri = 1 for W = (L, L), W = (−L, 0) and W = (0, −L), while W tri = 0 for W = (0, 0). Alternatively one can track the Monte Carlo average of the density of strings, n string (see Section B.2 for the definition of strings). Correlations can be understood by measurement of the spin structure factor, defined in the usual way in real and reciprocal space as, where r = r i − r j and q is in the Brillouin zone of the triangular lattice.

String/fermion mapping
In order to gain intuitive insights into the TLIAF that complement the Monte Carlo simulations, it is useful to consider some of the mappings that can be made. One option is to make a mapping onto a height model, which describes the configurations of the TLIAF in terms of a coarse-grained height field, and this is particularly powerful way to capture the long-wavelength features of the nearest-neighbour TLIAF [36].
For the questions we are interested in here, we find it more intuitive to make a mapping onto strings [37,38], which can then be interpreted as the worldlines of fermions. The mapping onto strings requires the choice of one of the 3 principal lattice directions as being special, and while this is natural for the anisotropic TLIAF the decision is arbitrary for the isotropic model. The strings live on the dual honeycomb lattice, whose bonds bisect those of the triangular lattice (see Fig. 2 and also Appendix B for the link to dimer mappings). Along the special direction each honeycomb bond bisecting an antiferromagnetically-aligned, triangular-lattice bond is assigned to be a segment of a string, while in the other two directions honeycomb bonds bisecting ferromagnetically-aligned, triangular-lattice bonds are assigned as string segments. The string-free configuration is thus seen to correspond to an Ising stripe configuration, with stripes of aligned Ising spins parallel to the special direction.
Using the above definition, each honeycomb-lattice site can be touched by either 0 or 2 string segments, and this ensures the continuity of the strings. For a system with periodic boundary conditions the strings form closed loops, and no two strings can touch, let alone cross, one another (see Fig. 2). If a reference line is chosen that winds the system, the number of strings crossing it has to be even, and therefore string parity is conserved. If no defect triangles are present, then there is no way for a string to turn back on itself, and strings both wind the system and have a fixed length, and in this sense the strings are taut. Defect triangles act as sources or sinks of pairs of strings, allowing them to turn back on themselves and thus either form local loops or long floppy strings that wind the system.
In the absence of defect triangles the string degrees of freedom provide a way of labelling the Ising configurations according to a pair of winding numbers. Two reference lines can be chosen that wrap around a periodic cluster of linear size L, and the number of strings crossing each reference line is related to the associated winding number according to W i = [L − no. of strings crossing ref. line i] (see Appendix B for the link to dimer representations of the Ising variables). In order to transition between Ising configurations with different winding numbers, it is necessary to either create a pair of defect triangles and transport one of them around the system before they recombine, or to make a non-local change of the Ising configuration.
The properties of the strings allow them to be interpreted as the worldlines of spinless fermions, as is common for 2D statistical mechanics problems [39,40,41,42,43,44]. The strings "travel" in the direction perpendicular to the special lattice direction, and this is interpreted as the imaginary time direction (see Fig. 2). At a mininum, the quantum Hamiltonian has to include a chemical potential measuring the energy cost of creating a string/fermion, a hopping term that allows the strings to move in the direction perpendicular to that of imaginary time and pair creation and annihilation terms that take into account the effect of defect triangles. For the case of the nearest-neighbour triangular lattice these three terms are sufficient, and there is an exact mapping onto, where the parameters of the 1D quantum model can be expressed in terms of those of the 2D classical model (see Appendix C, Appendix D and Appendix G for details). The simplicity of the fermionic model is due to the fact that the string-string interactions in the nearestneighbour Ising model are purely entropic -they arise from the non-crossing constraint -and this maps onto the fermionic Pauli exclusion principle. Once further-neighbour Ising interactions are included, it is necessary to take into account energetic string-string interactions, and these map onto fermionic interactions. However, a phenomenological free-fermion model can still be applicable when the string/fermion density is low, and can provide quantitative insights into the critical behaviour of the Ising model. Futhermore, qualitative insights into the Ising system can be gained from considering the form of the fermion-fermion interactions.

Monte Carlo simulation of the dipolar TLIAF
Here we determine the main physical features of the dipolar TLIAF using Monte Carlo simulation. A more detailed discussion of their physical origin is postponed until Section 5.
The ground state of the dipolar TLIAF is 6-fold degenerate and consists of alternating stripes of equivalent Ising spins running parallel to A, B or C bonds (see Fig. 5). The 3-fold degeneracy associated with the choice of stripe direction is multiplied by a 2-fold Ising degeneracy associated with a global spin flip, giving the overall 6-fold degeneracy.
At low temperature there is a stripe-ordered phase, which is dominated by the ground state configurations. Local fluctuations are highly suppressed because they involve the creation of pairs of defect triangles, and the associated energy cost is large. In principle it is also possible to create strings that wind the system, but these are forbidden in the thermodynamic limit as they cost a finite free energy per unit length [10,11]. The maximum of χ W scales approximately with 1/N , as is standard for a 1st order phase transition, leading to our estimate that in the thermodynamic limit the transition temperature is T1/J1 = 0.1845 ± 0.0010. (b) Energy-histogram analysis of an L = 36 cluster for temperatures close to the χ W -maximum of T1/J1 ≈ 0.183. Energies are measured in units of J1 and energy bins have width 0.0005J1. A sharp low-energy peak, associated with an almost fluctuationless low-temperature phase, is separated from a Gaussian peak associated with the high-temperature phase by an energy gap of approximately 0.03J1 per site. The separation of the two peaks is evidence of a first-order transition.
On further increasing the temperature, there is a transition from the stripe phase into a disordered phase. A previous study determined the transition temperature to be T /J 1 ≈ 0.18, but was not able to determine the nature of the transition or achieve equilibriation across the transition [45]. Our simulations, which do achieve equilibriation, show that the transition is first order, and this is clear from histogram analysis of the energy close to the transition temperature, as shown in Fig. 3. The transition temperature can be determined from the peak in the heat capacity, C, the order parameter susceptibility χ stripe [Eq. 2] or the winding number susceptibility χ W [Eq. 3]. The positions of the peaks in these different quantities coincide for a given system size, L, and we show results for χ W in Fig. 3. The position of the peak shows a weak L dependence, and using the standard scaling of a first-order transition temperature with 1/N , we determine a transition temperature of T 1 /J 1 = 0.1845 ± 0.0010.
The first-order nature of the transition is also clear from simulations of the heat capacity, which are shown in Fig. 4. Integrating C/T from infinity shows that the disordered state just above the phase transition has an entropy per site of S/N ≈ 0.22. While this is less than the Wannier entropy S wan /N = 0.323 . . . associated with the ground state of the nearestneighbour TLIAF [2], it is still considerable. The low-temperature stripe phase is essentially fluctuationless and has S/N = 0, showing that there is a significant entropy release in the first-order transition.
The nature of the correlations can be accessed via the spin stucture factor, and a representative set of examples are shown in Fig. 5. In the stripe-ordered phase there are Bragg peaks Error bars are smaller than the point size unless explicitly shown. (a) The heat capacity (red) shows a broad maximum centred on T = 0.8J1, which corresponds to the freezing out of defect triangles, and a sharp peak at T = 0.185J1 due to a first-order phase transition. This can be compared to the case of the nearest-neighbour TLIAF (blue), which shows a similar broad maximum centred on T = 1.2J1, but no low-temperature phase transition. (b) The entropy per site (red) is calculated by integrating the heat capacity from infinity. The entropy passes though Swan/N , at T = 0.5J1, but does not show an extended plateau, unlike the nearestneighbour TLIAF (blue). At the phase transition there is an entropy jump of ∆S/N ≈ 0.22. (c) The defect triangle density, n def . (Inset) At low temperatures n def follows an activated behaviour, and the blue line is the best fit to Eq. 29 with A = 0.24 and E def = 1.60 (see Appendix A). associated with the three different stripe directions, and these occur at q stripe = (0, 2π/ √ 3) and symmetry related wavevectors. At temperatures just above the transition, the structure factor develops weight on the whole of the Brillouin zone boundary, but remains peaked at q stripe , despite the significant entropy change. Small further increases in temperature result in the growth of sharp structure-factor peaks at q tri = (2π/3, 2π/ √ 3), and the disappearance of peaks at q stripe . Further increasing the temperature results in the structure-factor weight becoming more diffuse, and the peaks at q tri become less sharp.
The crossover from a highly-correlated paramagnet with sharp structure-factor peaks to a weakly-correlated paramagnet with a diffuse structure factor is governed by the presence or absence of defect triangles. The temperature evolution of the density of defect triangles, n def , is shown in Fig. 4, where n def is defined as the total number of triangular plaquettes with three equivalent spins divided by the total number of plaquettes, 2N . Just above the transition temperature the density is very low, and at T /J 1 = 0.2, one finds n def ≈ 10 −4 . On the other hand, in the uncorrelated, infinite-temperature limit the defect-triangle density saturates at n def = 0.25, since triangles can take 8 possible equally probable configurations, 2 of which have three spins aligned. For simplicity we take the crossover from strong to weak correlation to be at n def = 0.025 (i.e. 10% of the saturation value) and this occurs at T = 0.75J 1 , which matches well to the broad peak in the heat capacity (see Fig. 4).
It is possible to estimate the typical energy of an isolated defect triangle by making fits to n def in the temperature range T J 1 . This shows activated behaviour, and a crude estimate of the functional form is derived in Appendix A. The best fit, shown in Fig. 4, corresponds to a defect-triangle energy of E def = 1.60J 1 . This shows that one effect of the further-neighbour interactions is to slightly decrease the typical energy of a defect triangle relative to the nearestneighbour TLIAF, where E def = 2J 1 . 3) and symmetry-related wavevectors, associated with the stripe ordering. (b) Above the transition weight develops around the BZ boundary, but is peaked at q stripe and q tri = (2π/3, 2π/ √ 3), and this is associated with a domain-wall network configuration. (c) Further increasing the temperature results in the peaks at q tri dominating over those at q stripe , and this is associated with a switch from attractive to repulsive string-string interactions. (d) At still higher temperatures the peaks at q tri remain sharp, and the system can be described as a string-Luttinger liquid. (e) Once the temperature becomes comparable with J1 the peaks at q tri lose their sharpness, and this is associated with the proliferation of defect triangles and the loss of significant correlation.

Monte Carlo simulation of the anisotropic, dipolar TLIAF
Next we turn to the distorted triangular lattice, and show that even quite small distortions can lead to significant changes in the physical behaviour compared to the isotropic lattice.
The distortion, which is parameterised by δ, leaves the length of A bonds invariant, while reducing the length of B and C bonds and therefore breaks the 6-fold ground-state degeneracy of the isotropic lattice down to a 2-fold degeneracy. Stripes form parallel to A bonds, and the 2-fold, ground-state degeneracy is simply due to an Ising degree of freedom, associated with a global flip of all the spins.
As in the isotropic case, the transition from the stripe phase to the disordered phase can be located using the peaks in the heat capacity, C, the order parameter susceptibility χ stripe [Eq. 2] or the winding number susceptibility χ W [Eq. 3], and the resulting phase diagram is shown in Fig. 6. Histogram analysis of both n string and E show that the transition changes from first order at low δ (see Fig. 3) to second order at high δ (see Fig. 6), and the change occurs at δ tri ≈ 0.02. However, this type of analysis is not a very precise gauge of δ tri , due to both finite-size effects and the fact that the first-order nature of the transition becomes weaker approaching δ tri . In Section 5 we use finite-size scaling analysis to determine how the critical exponents depend on δ, and thus demonstrate that the change from first to second order occurs via a tricritical point located at δ = δ tri = 0.022 and T = T tri = 0.343.
The spin-liquid region, in which strong local correlation co-exists with long-range disorder, is found to extend until approximately δ ≈ 0.15, with the associated temperature window decreasing with increasing δ. This is shown in Fig. 6, where we continue to use a defecttriangle density, n def = 0.025 (10% of the saturation value) to signify the upper extent of the spin liquid.
The structure factor shows signs of a second order transition for δ > δ tri and can also be used to characterise the disordered region. For δ tri ≤ δ 0.1 satellite peaks appear at the transition either side of q stripe = (0, 2π/ √ 3(1 − δ)), and gradually shift towards q tri = (2π/3, 2π/ √ 3(1−δ)) as the temperature is increased. We will argue below that these follow the string density, and the structure factor is peaked at q = q string = (πn string (T, δ), 2π/ √ 3(1−δ)). In the spin-liquid region these peaks are sharp, and this is associated with a long spin-spin correlation length. In the weakly-correlated region, the structure factor becomes diffuse, signalling the breakdown of strong correlations and a short spin-spin correlation length. For δ 0.15 the spin-liquid region is totally suppressed and, above the transition, the structure factor has peaks at q stripe . In this region the critical behaviour shows the characteristics of a usual second-order Ising transition into a symmetry-broken, stripe phase, with a structure factor peak developing in the disordered region at the ordering vector.

Discussion and analysis
In order to gain physical insight into the Monte Carlo simulation results presented in Sections 3 and 4, it is useful to analyse them in terms of the string model introduced in Section B.2. We first discuss the nature of the phase transitions and then move on to the nature of the correlations within the classical spin-liquid region.

The nature of the phase transitions
Depending on the value of the anisotropy parameter, δ, the phase transition has different character, with a first-order transition at δ < δ tri turning into a second order transition at δ > δ tri via a tricritical point at δ = δ tri (see Fig. 6).
The nature of the second-order transition for δ > δ tri is complicated by the fact that it shows a combination of Pokrovsky-Talapov and Ising criticality, with the details depending on δ and T . Here we show that the criticality is 2D Ising over some potentially narrow temperature window close to the transition, and then crosses over to Pokrovsky-Talapov outside this region (see Appendix D for a discussion of similar behaviour in a simpler setting). at q string that coexist close to the phase transition with an additional peak at q stripe . At higher temperatures the peaks at q stripe are suppressed, leaving the peaks at q string more visible. (Bottom row) At large values of δ (here δ = 0.2) the spin liquid region is absent, and S(q) is peaked at q stripe , displaying the usual behaviour associated with a second-order transition.
The width of the Ising temperature window is exponentially suppressed as δ decreases and for small δ (e.g. δ = 0.05) the transition is to all intents and purposes in the Pokrovsky-Talapov universality class (i.e. is of Kasteleyn type).
First we consider the extreme case where defect triangles are completely absent and the transition is strictly in the Pokrovsky-Talapov universality class (see also Appendix C). While this is never rigorously true in the dipolar TLIAF, it is a good approximation at low T . The partition function of the 2D classical model can be mapped onto that of a 1D quantum model with the Hamiltonian, where T c is the critical temperature of the dipolar TLIAF, and a 0 and b are phenomenological parameters. In the simpler case of the nearest-neighbour TLIAF, the exact 1D quantum Hamiltonian is given in Eq. 5, and matching this to Eq. 6 requires µ + 2t → a 0 (T − T c ), t → b and ∆ → 0. For the dipolar TLIAF such a microscopic matching of parameters is not possible, but the phenomenological dispersion given in Eq. 6 is valid as long as the probability of creating a defect triangle is very low. Furthermore, truncation of the dispersion beyond the q 2 term remains a good approximation as long as the string density is low. For the fermion model, the T that appears in Eq. 6 does not have the meaning of temperature, and is simply a tuneable parameter that controls the transition from an insulating phase with no fermions (T < T c ) to a metallic phase with gapless excitations (T > T c ) at a fermi wavevector q f = −a/b. Once this is mapped back to the dipolar TLIAF, T regains the meaning of temperature, the insulating fermionic phase corresponds to the string vacuum (the stripe phase) and the metallic fermionic phase to the spin-liquid phase where there is a finite density of fluctuating strings that wind around the system.
The density of strings in the 2D classical model is equal to the density of fermions in the 1D quantum model, and can be calculated as, In the fluctuating phase the string density can be expressed as n string ∝ (T − T c ) β , with the critical exponent β = 1/2. This is characteristic of Pokrovsky-Talapov-type critical behaviour [39,40] associated with a Kasteleyn transition [46].
To make a connection with previous work [40], it is useful to determine the free energy of the 2D classical model, which is just given by the energy of the 1D quantum model, resulting in, where n string = q f /π. In Ref. [40] it was shown that the cubic term describes the string-string repulsion.
For the dipolar TLIAF, the above analysis should apply to the second-order phase transition at anisotropy values δ δ tri , where the transition temperature is low enough that there are very few defect triangles in the system. In order to test this, we perform simulations of H dip [Eq. 1] for δ = 0.05 at a range of system sizes. While finite-size effects make it hard to directly measure the exponent β in simulations, it is possible to write down a scaling hypothesis for n string and use this to determine β. We consider, where g PT is an unknown scaling function and ζ is the correlation length in the direction parallel to the domain walls, defined as the lengthscale at which asymptotic values of the structure factor becomes valid [47]. In the critical region it is expected that ζ ∝ (T − T K ) −ν with ν = 1, and this can be compared to the correlation length perpendicular to the domain walls, which is given by ζ ⊥ ∝ (T − T K ) −ν ⊥ with ν ⊥ = 1/2 [40,48,47]. Since we typically consider hexagonal shaped clusters, finite-size effects will be dominated by ζ , since this diverges faster than ζ ⊥ . In order to quantitatively test the goodness of the data collapse according to the scaling hypothesis [Eq. 9], we use the measure proposed in [49]. The best collapse was found for β = 0.44 ± 0.07 and ν = 0.92 ± 0.2 (see Fig. 8), which is consistent with the expected Pokrovsky-Talapov exponents of β = 1/2 and ν = 1. Thus we conclude that at low values Ising tricritical Pokrovsky-Talapov Ising In the tricritical region the best data collapse is found for m stripe at δ = 0.022 using the scaling hypothesis given in Eq. 14 with Ising tricritical exponents β = 1/4 and ν = 1. (b) Data collapse at δ = 0.05 using the scaling hypothesis given in Eq. 9 for Pokrovsky-Talapov critical behaviour and for exponents that minimise the "goodness of collapse" measure proposed in [49]. This gives β = 0.44 ± 0.07 and ν = 0.92 ± 0.2, which are consistent with the expected Pokrovsky-Talapov exponents β = 1/2 and ν = 1. (c) Data collapse at δ = 0.2 using the scaling hypothesis given in Eq. 14 and the Ising critical exponents β = 1/8 and ν = 1.
of δ the second-order transition shows critical behaviour associated with a Kasteleyn-type transition, which is driven by the appearance of non-local strings that wind the system. In reality H dip [Eq. 1] supports a small density of defect triangles at any finite temperature, and even a tiny density of defect triangles drives the transition to be in the Ising universality class (see also Appendix D). However the temperature window over which Ising criticality applies is exponentially suppressed at small δ. In order to better understand the nature of the suppression and the crossover between Ising and Pokrovsky-Talapov criticality, we consider the phenomenological 1D quantum Hamiltonian, where z def = e − E def T and E def is a measure of the energy cost of a defect triangle. Diagonalisation via a Bogoliubov transformation, results in, where a q and a † q are fermionic operators. In terms of fermions the parameter T controls a transition from a gapped, insulating phase at T < T c to a gapped, p-wave-superconducting phase at T > T c via a gapless point at T = T c . In terms of the 2D classical model this maps onto the transition from the stripe phase at T < T c to a phase with fluctuating strings at T > T c .
The Ising/Pokrovsky-Talapov nature of the criticality is encoded in the location of the minimum of ω q [Eq. 11]. Ising criticality is associated with a minimum at q = 0, and this occurs for T c < T < T Is , where, T Is = T c + 8z 2 def /(a 0 b). The 2D Ising nature of the criticality in this temperature window is clear from considering the correlation length, which goes as ξ Is ∼ |T − T c | −1 [43]. For T > T Is the dispersion minimum moves away from q = 0 to q min = (−a/b − 8z 2 def /b 2 ) 1/2 and the system enters the crossover region between Ising and Pokrovsky-Talapov universality. Pure Pokrovsky-Talapov critical behaviour is recovered in the limit T − T c T Is − T c , where q min = (−a/b) 1/2 is recovered, and therefore n string ∝ (T − T c ) 1/2 . In the case T E def , the temperature width of the Ising window is exponentially suppressed due to the z 2 def factor, and the system shows Pokrovsky-Talapov characteristics over all accessible temperatures.
In the critical region the string density and its derivative are given by, In the region T c < T < T Is it is possible to extract analytic expressions for these quantities in terms of elliptic integrals. However, these expressions are not so enlightening, and we instead show a numerical evaluation in Fig. 9. For z def = 0, the string density is finite both above and below the transition and takes the value n string = 2z def /(πb) at T = T c . As such n string is not technically a good order parameter, but it still provides a useful indicator of the transition temperature since dn string /dT shows a logarithmic divergence according to dn string /dT ∝ log |T − T c | (see Fig. 9).
At intermediate values of δ (e.g. δ ≈ 0.1), the dipolar TLIAF should show a crossover between Ising and Pokrovsky-Talapov criticality. As is standard in such situations, an exponent φ can be used to parametrise the crossover [50,51]. This appears in the scaling form of physical quantities, and we consider the defect-triangle density, which is expected to scale as, where α is associated with the Pokrovsky-Talapov critical behaviour and therefore expected to take the value α = 1/2 [48]. As shown in Fig. 10, scaling collapse of Monte Carlo simulation data in the range 0.06 ≤ δ ≤ 0.1 works well for α = 0.5 and φ = 1.12. This compares well to a similar scaling analysis of the nearest-neighbour model, where n def can be calculated directly in the thermodynamic limit, and we find α = 0.5 and φ = 1 [see Appendix D].
At large values of δ there is a high density of defect triangles at the transition, and one expects Ising criticality to apply over a wide temperature window. That this is indeed the case can be shown by analysing simulation data at δ = 0.2. The standard scaling hypothesis for the Ising order parameter, m stripe [Eq. 2], is, and it is expected that data collapse occurs for β = 1/8 and ξ Is ∼ |T − T c | −1 . It can be seen in Fig. 8 this results in good collapse of the simulation data. At δ = δ tri = 0.022 there is a tricritical point, and the critical behaviour is different from that of the second-order transition. Since the transition temperature at the tricritical point is low, one would naively expect that the associated low density of defect triangles would result in Pokrovsky-Talapov tricritical behaviour (see Appendix E for a discussion of Pokrovsky-Talapov tricriticality). Pokrovsky-Talapov tricriticality can be described by the 1D quantum dispersion, where a > 0, c > 0 and b(δ − δ tri ) is an odd function of δ − δ tri that changes sign when δ = δ tri . It follows that exactly at the tricritical point (b=0), resulting in a critical exponent of β = 1/4. In terms of the free energy of the 2D classical model, the tricritical point occurs when the cubic term disappears, resulting in, Since the cubic term controls the string-string interaction, changing its sign is equivalent to going from a repulsive interaction associated with a second-order transition to an attractive interaction associated with a first-order transition.
For the tricritical point to be effectively in the Pokrovsky-Talapov-tricritical universality class, it is necessary that the Ising temperature window is negligible. The problem with this is that the expression we previously calculated, Including the fourth order term in ω q [Eq. 15] results in T Is = T c + 6z The exponential suppression of the Ising temperature region with E def /T is thus less pronounced than at the critical point, and, depending on the value of c, this could in principle lead to a wide temperature window of Ising tricriticality.
Monte Carlo simulations allow us to test whether Pokrovsky-Talalapov or Ising tricriticality dominates, and come down in favour of a significant Ising-tricritical window. This can be seen in Fig. 8, where m stripe [Eq. 2] shows convincing data collapse for 2D Ising tricritical exponents. At higher temperatures the system presumably crosses over to Pokrovsky-Talapov tricriticality, but the Ising temperature window is wide enough that this is difficult to ascertain.
For δ < δ tri the transition is first order. In this situation the string-string interaction is attractive, and the string density jumps at the transition. Expansion of the free energy in terms of the string density is therefore only possible close to the tricritical point, where the jump in the string density is relatively small.
Close to the first-order transition the effective fermion degrees of freedom are long lived, since decay of fermions is associated with the, essentially negligible, presence of defect triangles in the 2D classical model. As a result, even though the fermions are strongly interacting, the interactions just renormalise the free-fermion terms in the Hamiltonian, but don't generate new terms. Thus one can still think in terms of an effective free-fermion dispersion, ω q . However, the deeper one goes into the first-order region the larger the jump in the string density and the more (even) powers of q have to be retained in the expansion of ω q . This is because in order to generate a jump in n string a finite region of q values must have a flat dispersion with ω q = 0 at the transition, and the larger this region, the more powers of q are required to capture it effectively (strictly all powers of q are required for a region of ω q = 0, but close to the tricritical point this effect can still be essentially captured by a finite expansion). The predictive power of the phenomenological theory thus reduces away from the tricritical point due to the rapid increase in the number of coefficients.
Taking all the results of this section together, one can see that with relatively few parameters, it is possible to understand the full gamut of critical behaviour in the dipolar TLIAF. The important parameters are the reduced temperature, (T − T c )/T c , which changes sign at a second-order phase transition, the distortion-dependent parameter b, which changes sign at the tricritical point, and the ratio E def /T c , which determines the temperature window of Ising criticality at the transition.

Correlations between spins
Next we turn to the spin correlations, the nature of which can be used to attain a more detailed understanding of the phase diagram. These can be probed via the spin structure factor, S(q) or S(r) [Eq. 4].
The correlations are very simple in the stripe-ordered phase, which has Bragg peaks in S(q) at q = q stripe = (0, 2π/ √ 3(1 − δ)), and symmetry-related wavevectors (see Fig. 5 and Fig. 7). Since at low temperatures the stripe phase is essentially fluctuationless, virtually all the spectral weight is contained in the Bragg peak, and in real space there is no decay of the spin correlations with separation. For larger values of δ the stripe-ordered phase survives to higher temperatures, and for T ∼ J 1 local fluctuations around the stripe ground state associated with pair creation of defect triangles become significant, resulting in some diffuse scattering surrounding the Bragg peak.

Domain wall network
String Luttinger liquid Paramagnet Figure 12: Typical string configurations in the domain-wall-network, string-Luttinger-liquid and paramagnetic regions. In the domain-wall-network region strings typically wind the system and attractive string-string interactions cause them to bind together. In the string-Luttinger-liquid region the strings also tend to wind the system, but repulsive interactions result in a grill-like superstructure with strings avoiding one another as far as possible. In the paramagnet there are many defect triangles that act as sources and sinks of pairs of strings, resulting in the strings being floppy and forming short closed loops.
More interesting is the disordered phase, which shows three qualitatively different regimes of spin correlations. At high T there is a weakly-correlated paramagnet in which the spin correlations are short ranged, at lower T and for δ < 0.15 there is a strongly-correlated regime with longer-range correlations that we name a string Luttinger liquid and a for T < T tri and δ < δ tri there is a different strongly-correlated regime that we call a domain-wall network (see Fig. 11).
It is instructive to discuss each of these regimes in more detail, and first we turn to the string Luttinger liquid (see also Appendix E). In this regime the density of defect triangles is low, and there is a repulsive interaction between the strings. Since the strings repel one another they form a (disordered) grill-like superstructure where the average spacing between the strings depends on the string density, n string (T, δ) (see Fig. 12). As a result of this superstructure, the structure factor is peaked at q = q string (T, δ) = (πn string (T, δ), 2π/ √ 3(1 − δ)) and related wavevectors (see Fig. 7). However, since the strings are fluctuating, the peaks are not Bragg peaks, and in real space spin correlations decay to zero for large enough separations.
In the absence of defect triangles spin correlations in real space decay algebraically, while in the presence of defect triangles the decay is exponential at large enough distances. We make the ansatz that the real-space, spin-correlation function takes the asymptotic form, where ξ ⊥ and ξ are correlation lengths in the directions perpendicular and parallel to the strings (in the isotropic case ξ ⊥ = ξ ). The exponential nature of the decay only becomes apparent when the spin separation is comparable to the correlation length ξ ⊥ or ξ . For low defect-triangle densities, this correlation length is typically large (it becomes infinite when the defect triangle density goes to zero) and for r ≤ ξ the spin correlations are essentially algebraic.
The parameter K detemines the speed of the algebraic part of the decay, and is nothing but the Luttinger parameter familiar from 1D fermionic systems [52]. That this should appear is unsurprising given the mapping between strings and spinless fermions (see Section B.2) and due to the fact that a repulsive interaction between strings maps onto a weakly attractive interaction between fermions, one expects K > 1 (for comparison non-interacting fermions have K = 1).
In practice simulations show that the reciprocal-space structure factor in the string-Luttinger-liquid region is dominated by peaks at q = q string , but there is also spectral weight on the line of q values joining q string and −q string (see Fig. 7). In consequence it is better to fit the structure factor in reciprocal space than in real space, and Fourier transforming the asymptotic form given in Eq. 18 in the vicinity of q = q string and for ξ ⊥ = ξ = ξ gives [53], where p = q − q string , p = |p| and, where J 0 (x) is the Bessel function of the first kind, Γ(x) is the Euler Gamma function and is the hypergeometric function. The result of fitting this to simulations for δ = 0 and appropriate T is shown in Fig. 11, and it can be seen that the Luttinger parameter does indeed take values K > 1, and the correlation length can be many multiples of the lattice spacing. A more precise determination of K and ξ would require simulations on larger clusters (for a numerical determination of K in a simpler model see Appendix E). It is clear from the simulations of S(q) shown in Fig. 5 that the string-Luttinger liquid regime does not survive all the way down to the phase transition when δ < δ tri . Rather than being peaked at q = q string , the low-temperature structure factor in the disordered region has spectral weight spread around the BZ boundary, and in particular weight starts to develop at q = q stripe . As was shown in Ref. [11] this type of stucture factor is associated with the formation of sizeable domains of stripe order, with neighbouring domains having stripes along different principal axes (see Fig. 12).
The formation of large stripe-ordered domains is suggestive that within this regime the strings attract one another (see also Appendix E). It makes sense that the crossover from the string-Luttinger-liquid regime (high T , repulsive string-string interactions) to the domainwall-network regime (low T , attractive string-string interactions) should occur at approximately T = T tri , and this is consistent with the S(q) measurements (see Fig. 5). A simple way to test that this is the case is to perform simulations in a highly restricted manifold of Ising configurations containing two strings, each of which winds the system. The average separation of the strings does indeed show a significant drop starting at T ≈ T tri , indicating a shift from repulsive to attractive interactions (see Fig. 11). We find that the temperature at which this change occurs is essentially independent of δ in the relevant region (0 < δ < δ tri ), and therefore the crossover between the string-Luttinger-liquid and domain-wall-network regimes is approximately flat, as shown in Fig. 11.
In terms of fermions, the domain-wall network state can be thought of as being a fluctuating, phase-separated state, with a loose analogy to the clustering of holes in superconductors [54,55].
At high temperatures the system forms a weakly-correlated paramagnetic regime. The correlations can still be described by Eq. 18, but the correlation length is comparable to the lattice spacing, and so the correlations are exponentially decaying at all length scales. While we have defined the crossover from the Luttinger liquid regime to the weakly-correlated regime in terms of the density of defect triangles reaching 10% of its saturation value (see Fig. 6) this is roughly equivalent to defining a crossover in terms of the correlation length reducing to about 2 lattice spacings.
For δ 0.15 there is a direct transition from the stripe-ordered phase to a standard paramagnet, and as such the structure factor shows the usual features of a second-order transition, with spectral weight building up at the ordering vector as the transition is approached from above, and a diverging correlation length at the transition that results in the formation of a Bragg peak.

Triangular lattice antiferromagnets with general couplings
Here we take a step back and discuss the general features of TLIAF models with monotonically decreasing further-neighbour interactions. We have argued that a good way to understand such models is in terms of the string degrees of freedom, which can be thought of either in their 2D classical incarnation or as the worldlines of spinless fermions in 1D. As such we would like to determine which energy scales present in a given microscopic model dictate the behaviour of the strings and therefore the form of the phase diagram and the physical observables.
In general TLIAF models have many competing couplings, as is clearly true in the dipolar case. Our claim is that these can in most cases be distilled into four important energy scales (it is worth noting that other energy scales can become important if the further-neighbour interactions are not monotonically decreasing interactions [10,11]).
The first and most important energy scale is the isotropic part of the nearest-neighbour interaction; that is the part of the nearest-neighbour interaction that does not vary with the anisotropy (i.e. J 1A ). This approximately sets the energy cost of creating defect triangles, and therefore "interesting", strongly-correlated physics only occurs in the region T < J 1A .
Next are two energy scales that combine the isotropic parts of the further-neighbour couplings. The first of these, J fn , is a measure of the internal energy of a string and also sets the string-string interaction energy scale. As an example, for the TLIAF with J 1 , J 2 and J 3 couplings it is given by J fn = J 2 − 2J 3 [10]. This shows that even if the further-neighbour couplings are comparable with J 1 , their combined effect can still be small due to frustration. The second energy scale is J c , and this is related to the energy cost associated with a string changing direction, and in the case of the J 1 -J 2 -J 3 model this is given by J c = J 2 . One thing that is important to note is that we always consider J fn , J c > 0, and if this is not the case different physics can be expected [11].
The final energy scale we consider is a measure of the anisotropy and is labelled J an . For the dipolar TLIAF it clearly depends on δ, and a rough estimate is given by the difference in the nearest-neighbour interaction strengths, resulting in, The energy scales J 1A , J fn , J c and J an have been constructed with the string degrees of freedom in mind, and we now make the link more explicit. We concentrate in particular on J 1A J fn , J an , which is the requirement for the existence of spin-liquid behaviour. A particularly important quantity is the internal free energy per unit length of an isolated string, which depends on J fn , J c and J an , and is approximately given by [10,27], If string-string interactions are ignored, strings will be present in the system above a temperature T string , and this is approximately given by, While a number of approximations have been made in order to arrive at this simple expression, except in the extreme case of J c J fn , it matches well to Monte Carlo simulations of simple models [11].
In reality the strings are not isolated, and the transition temperature and the nature of the correlations in the spin liquid depend on the string-string interactions. These have two main contributions, the first of which is an entropically-driven repulsion associated with the nocrossing constraint, and in the fermion language this maps onto the Pauli exclusion principle. The second is an energetically-driven attraction due to the further-neighbour interactions, which is approximately measured by J fn , and in the fermion language it is only this second contribution that counts as an interaction.
If the attractive interaction dominates in the vicinity of T = T string then an array of strings can lower their energy by binding together, and this binding energy results in a firstorder transition with T 1 < T string . As a result the string density jumps at the transition from n string ≈ 0 to a finite value. At temperatures just above the transition the string-string interactions remain attractive and the strings loosely bind together, forming a domain-wallnetwork state. The domain-wall-network state also relies on a positive J c which penalises changes in direction of the strings. The larger the value of J c and J fn relative to T , the larger the domain size will be. The dominance of attractive interactions in the string picture corresponds to the strong-coupling regime of the fermionic model.
When T J fn the entropically-driven repulsion between strings dominates over the energetically-driven attraction. If T string J fn then the strings repel one another in the critical region, resulting in a second-order transition at T = T string . As long as J 1A T string then this transition is essentially of the Kasteleyn type, since it is driven by the sudden appearance of strings that mostly wind the system. This type of phase transition is quite different from the more usual Ising transition which is driven by the proliferation of local defects. Above the second-order transition the string density, n string , increases with increasing temperature, and, while the strings fluctuate, they on average form an equally-spaced, grilllike structure due to their mutual repulsion. In the fermionic language this corresponds to weak coupling and a 2D classical equivalent of a Luttinger liquid forms.
When the attractive and repulsive interactions balance, the phase transition is tricritical, and this occurs when J an ≈ J fn . Just above the transition the string or fermion dipersions are soft, resulting in large fluctuations in the string/fermion density.
The crossover between the spin liquid and paramagnet occurs at T ≈ J 1A and at this temperature defect triangles become common. As a result strings form short closed or longer floppy loops which typically don't wind the system. If T string ≈ J 1A , then the transition is in the Ising universality class and is driven by the proliferation and growth of local defects, resulting in a direct transition from the ordered phase to the weakly-correlated paramagnet. In the dipolar TLIAF this occurs for δ 0.15.
For the dipolar TLIAF it is possible to approximately determine the appropriate energy scales as J fn ≈ 0.02J 1A , J an ≈ 9δJ 1A /4 and J c ≈ 0.08J 1A . Here J c is determined as half the energy cost of an isolated corner, while J fn is determined so as to be consistent both with T string [Eq. 23] and with Monte-Carlo, worm-update simulations, which are found to work best with approximately this value of J worm 2 − 2J worm 3 (see Section 2.1). Despite the slowly decreasing nature of the dipolar interaction with distance, it can be seen that frustration leads to a value of J fn that is considerably smaller than J 1A , resulting in a significant window in which the spins are strongly correlated.
An obvious question raised by this analysis is how to further reduce the value of J fn and J c relative to J 1A , since this would increase the size of the spin-liquid region and give a cleaner realisation of the Kasteleyn transition. One possibility would be to find systems with local interactions such that J 1 J 2 , J 3 . . . , but we are not currently aware of any such systems. A more realistic option is to change the nature of the long-range interaction such that J ij ∝ |r i − r j | −a , where a = 3 corresponds to the dipolar case. The possibility of changing a has been realised experimentally using trapped ions that naturally form a triangular lattice, and a was found to be tuneable in the range 0 < a < 3 [26]. Estimating the relationship between J fn , J c and a is complicated, due to the competition between the further-neighbour interactions, but it seems most likely that suppression of J fn would require the further-neighbour interactions to fall off faster than in the dipolar case, and therefore a > 3.
Another possibility is to add a small transverse magnetic field. This would tend to act in opposition to the further-neighbour interactions, since quantum fluctuations favour nearestneighbour-flippable configurations of Ising spins, while the stripe configuration is maximally unflippable. Therefore a transverse field would be likely to reduce the critical temperature by suppressing J fn .

Conclusion
We have shown that the dipolar TLIAF shows a variety of behaviours, with stripe-ordered, spin-liquid and paramagnetic phases. Furthermore, the nature of the spin-liquid region can be tuned by temperature between a "strongly-coupled" domain-wall network and a "weaklycoupled" string Luttinger liquid, where the strength of the coupling refers to a mapping to a 1D fermionic model. The addition of a small anisotropy allows the nature of the spin liquid to be further tuned, and this in turn changes the critical behaviour from first order to Kasteleyn-like, via a tricritical point with mixed tricritical-Ising and tricritical-Pokrovsky-Talapov characteristics.
We end with the hope that the physics we have described will soon be explored experimentally in artificial spin systems. In such a setting the physics of the isotropic dipolar TLIAF may be even richer, since it is likely that the dynamics will be too local to reliably find the stripe-ordered phase at low temperature, and instead the domain-wall network state will likely freeze to form a glassy state.
Funding information We thank the Swiss National Science Foundation and its SINERGIA network "Mott physics beyond the Heisenberg model" for financial support.

A Defect triangles in the dipolar TLIAF
In this appendix we construct a crude model for the density of defect triangles, n def , in the low-temperature paramagnetic state of the dipolar TLIAF. The aim is to justify the simple functional form of n def used to fit the Monte Carlo simulations in Fig. 4.
Defect triangles are constrained to occur in pairs, and can be considered to appear on top of microstates of the constrained manifold (configurations without defect triangles). We make the crude assumption that the energy cost of these defect triangles is only weakly dependent on position and has an average value E def . In this approximation, the total energy due to the defect triangles is given by, where N def is the number of defect triangles and interactions between defect triangles have been ignored. The number of ways N def defect triangles can be placed in the system with N plq triangular plaquettes is simply given by the binomial coefficient, and therefore the associated partition function is, where β = 1/T . The average number of defect triangles is given by, In the limit T E def / log N plq the defect triangle density is given by, in agreement with an exact calculation for the nearest-neighbour TLIAF. In the opposite limit of T E def / log N plq then for T E def one finds, While the above analysis is clearly highly simplified with respect to the true situation in the dipolar TLIAF, it suggests that at low temperature and on finite-size systems one should expect the density of defect triangles to obey the relationship, with A and E def fitting parameters. The result of fitting this to Monte Carlo simulations is shown in Fig. 4, and we find E def = 1.60J 1 in the isotropic, dipolar TLIAF. This can be compared with the nearest-neighbour TLIAF, where the energy cost per defect triangle is 2J 1 .

B Mappings and winding numbers
There are a number of possible mappings from Ising configurations of the TLIAF to dimer and string representations. Here we review the mappings used in this article and the links between them. In order to do this it is useful to define two different manifolds of Ising configurations, the unconstrained manifold that contains all possible configurations and the constrained manifold that contains only those configurations that are ground states of the nearest-neigbour TLIAF. The constrained manifold is clearly smaller than the unconstrained one, but is itself extensive [2].

B.1 Mapping to dimer coverings of the dual lattice
One useful mapping is from Ising configurations on the triangular lattice to dimer configurations on the dual honeycomb lattice [46]. We use this when constructing Monte Carlo worm updates [11]. The dual honeycomb lattice is constructed such that its bonds cut exactly one bond of the original triangular lattice (and vice versa), as shown in Fig. 13. If the triangular-lattice bond has two equivalent spins, then the honeycomb-lattice bond is covered by a dimer, while if the spins are inequivalent the honeycomb-lattice bond is left empty. The mapping between spin and dimer configurations is 2 → 1, since the dimer configuration is unaffected by a global flip of all the Ising spins.
Configurations within the constrained manifold (i.e. ground states of the nearest-neighbour TLIAF model) have one ferromagnetic bond per triangle, and therefore the number of dimers is fixed and equal to the number of triangular lattice sites, N . It follows that sites on the honeycomb lattice respect the usual dimer model constraint of being covered by exactly one dimer, as shown in Fig. 13(a).
In the unconstrained manifold, for each pair of defect triangles there are two additional dimers, and therefore the number of dimers is not fixed. The honeycomb-lattice site at the centre of a defect triangle is covered by three dimers, and therefore does not respect the usual dimer model constraint (see Fig. 13(b)).
For the unconstrained manifold of Ising configurations an alternative dimer mapping is possible, which is constructed such that the number of dimers is fixed and each site obeys the usual dimer-model constraint of being covered by exactly one dimer [56]. This involves extending the honeycomb lattice such that every original site is replaced by three new sites arranged in a triangle (see Fig. 13). Dimers are then placed on the original honeycomb lattice bonds in the same way as before, leaving a unique way of dimer covering the remaining sites of the extended honeycomb lattice such that every site is covered exactly once.

Reference
Actual

B.2 Mapping to string configurations on the dual lattice
The main mapping used throughout the article is onto string configurations on the dual honeycomb lattice [37,38]. Here we show how this is related to the dimer mapping described above. This proceeds by comparing a given dimer configuration to a reference configuration in which all the dimers are parallel (see Fig. 14). Any honeycomb-lattice bonds on which there is a discrepancy between the actual dimer configuration and the reference configuration is assigned to be part of a string. The chosen reference configuration consists of alternating horizontal stripes of aligned Ising spins, and this corresponds to all vertical bonds of the honeycomb lattice being covered by dimers (see Fig. 14(b)). This choice of reference configuration results in a number of useful properties of the strings, the most important of which is that strings never touch or cross. For periodic boundary conditions there is the additional constraint that the number of strings crossing an arbitrary reference line that winds the system has to be even, meaning that the string parity is conserved. If the Ising configurations are restricted to be in the constrained manifold the strings are directed, in the sense that they cannot turn back on themselves, and therefore have to wind the system, as shown in Fig. 14(c). In the unconstrained manifold defect triangles act as sources and sinks of pairs of strings, resulting in (non-winding) closed loops of strings as well as strings that turn back on themselves, as shown in Fig. 14(d).

B.3 Winding number sectors
In the presence of periodic boundary conditions, Ising configurations within the constrained manifold can be labelled by a pair of winding numbers.
One way to define the winding numbers, W = (W 1 , W 2 ), is to consider a pair of reference lines, as shown in Fig. 15. For each dimer crossing the horizontal part of the reference line the winding number is augmented by +1, and for each dimer crossing the angled part of the reference line it is augmented by −1.  In the string picture, the winding number is simply given by, and it follows that the density of strings in the constrained manifold can be written as, The string vacuum is therefore equivalent to the winding number sector W = (L, L) and the sector W = (0, 0) corresponds to a density n string = 2/3. The winding numbers split the constrained manifold into topological sectors, in the sense that it is not possible to move between configurations with different winding numbers by making a series of local spin flips. Instead it is necessary to flip clusters of spins that wind the system.
In the unconstrained manifold (defect triangles allowed) W remains a useful quantity, but is no longer strictly a winding number, since the creation of a pair of defect triangles on the reference line is a local move that alters W. Nevertheless, it remains a useful concept when the defect-triangle density, n def , is low.
In order to isolate and study some of the important features of general TLIAF's, we consider a number of simple models, in which the interactions are local and can be varied at will. The subject of this appendix is the simplest of these models, the TLIAF with anisotropic nearest-neighbour interactions and a constraint forbidding defect triangles. The purpose of studying such a model is to understand the Pokrovsky-Talapov critical behaviour [39,40] and the correlations within the spin-liquid phase in a simple setting. In terms of the anisotropic, dipolar TLIAF studied in the main text, the ideas will be particularly relevant to the phase transition in the region δ tri < δ 0.1 and to the correlations in the string-Luttinger liquid phase for T J 1A .
The solution of this model is already well known due to the fact it can be mapped to free fermions, and was studied by Wannier in the case of isotropic interactions [2], and can be transformed to the Kasteleyn model for anisotropic interactions [46]. The Hamiltonian is given by, where ij α denotes nearest-neighbour bonds in the α direction (see Fig. 1 for the definition of bond directions). An alternative parametrisation can be achieved by writing, and we consider the case δJ > 0 (equivalently J 1A < J 1B ). We also impose the constraint that defect triangles are forbidden, which corresponds to taking the limit J 1A /δJ → ∞.

C.1 Dimer mapping
H ABB [Eq. 32] can be mapped onto the Kasteleyn model of dimer coverings of the honeycomb lattice, which has an exact solution [46]. The mapping from Ising spins on the triangular lattice to dimers on the dual honeycomb lattice is described in Appendix B.1, and the energy of a dimer configuration is given by, where N bond = 3N is the total number of bonds (this is the same for the triangular and dual honeycomb lattices) and N α dim is the number of dimers covering α-type bonds. Since defect triangles are forbidden, the total number of dimers is fixed as N A dim + N B dim + N C dim = N bond /3. In the ground state N A dim = N bond /3 and N B dim = N C dim = 0, and therefore the energy of a given configuration relative to the ground state energy is, It follows that the partition function can be written, up to a configuration independent prefactor, as, where the sum is over all dimer coverings of the honeycomb lattice. It can be seen from Eq. 35 that the weight associated with dimer covering a B or C bond is given by,

C.2 Evaluation of the partition function
It has been known how to evaluate partition functions of the type Z hon [Eq. 36] for many years [46,6]. Here we will briefly sketch the solution, since it will prove a useful basis from which to consider more complicated models. The starting point is to introduce a real, anticommuting Grassmann variable at each site of the honeycomb lattice [6,57,58]. These variables obey the usual rules: a i a j = −a j a i , da i = 0 and da i a i = 1. Since the honeycomb lattice has a 2-site basis, it is useful to label Grassmann variables as a and b on the two sublattices, and the partition function is therefore given by, where i labels unit cells and S 2 [a, b] is the Kastelyn action. Here K is a signed adjacency matrix, known as the Kasteleyn matrix [46], and contains the weights z [Eq. 37]. The reason for the appearance of det K rather than the more usual Pfaffian is that the matrix connects sites on different sublattices, but not those on the same sublattice.
To simplify the geometry the honeycomb lattice is distorted into the brick lattice, as shown in Fig. 16. Bonds are assigned a direction in accordance with the Kasteleyn theorem, which states that transition cycles should have an odd number of arrows in each sense [46]. The bond weights are assigned according to Eq. 36, with weight 1 on A bonds and weight z on B and C bonds. It follows that the action can be written as, where the coordinate system is defined by the unit vectorsê x andê y , as shown in Fig. 16.
The action is simply diagonalised by taking the Fourier Transform, to give, Finally the partition function can be evaluated as, where * k = −k has been used.

C.3 Physical properties
In order to understand better the physical properties of H ABB [Eq. 32] it is useful to notice that the free energy, which is given by, is typically dominated by the minimal values of | k |. The physical characteristcs of the system are therefore determined predominantly by the "low-energy" part of | k |. At low temperature, the spectrum of | k | is gapped at all k, as shown in Fig. 16, and this corresponds to the stripe-ordered state. The gap closes at k = (π, 0) at the temperature, and this corresponds to the temperature at which the free energy of strings goes to zero. At T = T K strings condense into the system, and there is a Kasteleyn transition out of the stripe-ordered phase and into the spin liquid. This is second order due to the non-crossing constraint of the strings, which results in an entropically-driven string-string repulsion. The transition is in the Pokrovsky-Talapov universality class [39,40].
For all T > T K the spectrum of | k | is gapless. The position of the gapless point moves from k = (π, 0) at T = T K to k = (2π/3, −π/3) at T → ∞ and the position of this point is simply related to the string density, which smoothly increases with increasing temperature. We show below that the gaplessness of the spectrum is associated with algebraic decay of the spin-spin correlations [4,5].
In order to detect the transition between the stripe-ordered phase and the paramagnet, one possibility is to measure the local stripe order parameter m stripe [Eq. 2]. However, this is somewhat unsatisfactory, since m stripe = 1 in the ordered phase, and there is a discontinuous jump to m stripe = 0 at T = T K (see Fig. 17). Thus m stripe does not show critical behaviour, and this is due to the fact that the transition is not driven by the proliferation of local defects, but by strings that wind the system.
A more useful physical quantity is the density of strings, n string , and this does show critical behaviour. However, unlike a conventional order parameter, n string = 0 in the ordered phase, and only takes a finite value for T > T K . It can most simply be calculated in terms of dimer densities, according to, where the normalisation is such that 0 ≤ n string ≤ 1. In the case of the constrained manifold, the total number of dimers is fixed (see Appendix C.1) and this leads to the simplified expression [37], where the second equality follows from the expression for Z hon [Eq. 36]. When working in the constrained manifold, the string density can be calculated in a simple closed form. Substituting the expression for | k | [Eq. 42] into n string [Eq. 46], making the change of variables p x = (k x + k y − π)/2 and p y = (k x − k y − π)/2, and taking the thermodynamic limit results in, n string = 1 π 2 π 0 dp x π 0 dp y 4z 2 cos 2 p x − 2z cos p x cos p y 1 − 4z cos p x cos p y + 4z 2 cos 2 p x = 1 π 2 π 0 dp x u 2 ∂ ∂u π 0 dp y log[1 − 2u cos p y + u 2 ], where u = 2z cos p x . The integral over p y is tablulated and given by [2], As a result one finds, and this is plotted in Fig. 17. It can be seen that for T → ∞ (z → 1) the string density saturates at n string = 2/3, while for T ≈ T K it shows Pokrovsky-Talapov critical behaviour with n string ∝ (T − T K ) β and β = 1/2.

C.4 Correlations
In order to better understand the correlations it is useful to study the spin-spin structure factor [defined in Eq. 4]. When doing this care should be taken not to confuse q, which denotes a reciprocal vector in the Brillouin zone of the triangular lattice, with k, which lives in the Brillouin zone of the brick lattice. The structure factor can be calculated in the thermodynamic limit using the Grassmann path integral approach, following the general method proposed in Ref. [59]. For δJ = 0 this reproduces the results of Ref. [4,5]. A detailed summary of the calculation is given in Appendix F, both for pairs of spins separated by an arbitrary number of A bonds (i.e. in the direction perpendicular to the strings, denoted r x ) and for separations orthogonal to A bonds (i.e. in the direction parallel to the strings, denoted r y ). In both cases the structure factor can be written as the determinant of a Toeplitz matrix, whose dimension is proportional to the separation between the spins. Exact expressions can be written for the matrix elements, but we find it necessary to calculate the determinant numerically.
In the case of isotropic interactions the structure factor takes the asymptotic form [4,5], where q = (±2π/3, 2π/ √ 3), as can be seen in Fig. 18. The algebraic decay of correlations shows that the T = 0 nearest-neighbour TLIAF is critical, and is on the verge of forming 3-sublattice order. The combination of long-range disorder and local correlation means that the system forms a classical spin liquid.
For δJ = 0 and T > T K the structure factor retains the long-distance functional form given in Eq. 50, but the wavevector becomes temperature dependent and is given by q = ±q string (T ) = (±πn string (T ), 2π/ √ 3). This is clearly physically sensible, since the strings separate regions in which Ising spins have opposite sign, and the oscillation of the correlation function in the direction perpendicular to the strings should therefore have a period given by the average string separation. Some typical examples are shown in Fig. 18. Also shown is S(q), which has pairs of algebraically diverging peaks at q = ±q string (T ). In the vicinity of these peaks S(q string + δq) ∝ |δq| −3/2 in agreement with the algebraic decay of S(r) [Eq. 50]. It can be seen that the critical nature of the correlations is not broken by a non-zero δJ as long as T > T K , and this is due to the constraint forbidding defect triangles.
Having stated that the structure factor has the functional form given in Eq. 50 in the long distance limit, it is useful to be more precise over what counts as long distance. This has been considered in the closely related field of adsorption of a gas onto a substrate, where there exist domain walls with similar properties to the strings of the TLIAF [47]. A correlation length can be defined beyond which the long-distance algebraic correlation function given in Eq. 50 applies. In the direction perpendicular to the strings it is intuitively obvious that this is given by the average string-string separation, and therefore ζ ⊥ ∼ 1/n string . In the direction parallel to the strings it can be argued that ζ ∼ 1/n 2 string [47]. In the case of δJ = 0 (or T δJ), where n string = 2/3, the correlation lengths, ζ ⊥ and ζ , are not much longer than a single lattice spacing, and the long distance asymptotic form of the structure factor is recovered for spins separated by only a few lattice spacings. On the other hand, for δJ = 0 and T ≈ T K the string density is low, and the correlation lengths become very large, especially in the direction parallel to the strings. For T → T K an expansion of n string [Eq. 49] shows that the correlation lengths diverge as ζ ⊥ ∼ (T − T K ) −ν ⊥ and ζ ∼ (T − T K ) −ν , with ν ⊥ = 1/2 and ν = 1.
For T < T K the spin structure factor clearly does not follow the functional form of Eq. 50. Instead it is constant in real space, since the stripe state admits no fluctuations, and has Bragg peaks in reciprocal space. On cooling through T = T K the pair of algebraically-diverging peaks in S(q) coalesce to form a single Bragg peak at q = q stripe = (0, 2π/ √ 3) (see Fig. 18).

C.5 Mapping to 1D quantum model
A slightly different perspective on the nearest-neighbour TLIAF is achieved by making a mapping onto a 1D quantum model of spinless fermions. The idea is that the strings can be viewed as the worldlines of spinless fermions, and the spatial direction parallel to the strings interpreted as imaginary time. The non-crossing constraint of the strings corresponds to the Pauli exclusion principle, and periodic boundary conditions in the 2D classical model enforce periodicity in imaginary time in the 1D quantum model. This type of mapping has been frequently used for related 2D classical models with non-crossing domain walls [39,40,41,42,43,44]. We show in Appendix G that H ABB [Eq. 32] maps exactly onto the 1D quantum model, The mapping demonstrates that the Grassmann variables, a i and b i , describe coherent states of fermions/strings (for details see Appendix G). H 1D [Eq. 51] is simply diagonalised by Fourier transform, giving, and the phase diagram of the fermion model can be matched to that of the nearest-neighbour TLIAF. For µ/2t < −1 there are no fermions in the system and this is analagous to the stripe phase in which there are no strings. At µ/2t = −1 there is a phase transition due to the minima of the fermion band touching zero, and for µ/2t > −1 the fermion density, n f , is given by, n f = 1 − arccos[µ/2t]/π, which can be seen to be exactly equal to n string [Eq. 49].
According to the mapping given in Eq. 51, µ and t are not independent parameters, and the maximum value of their ratio is given by µ/2t = 1/2. This corresponds to z = 1 (equivalently T → ∞) and at this point n f = n string = 2/3 as expected. Other physical quantities, such as the heat capacity or the spin-spin structure factor can be calculated within the 1D fermion picture, and in some cases this simplifies the procedure.
It should be noted that if the 2D classical model has periodic boundary conditions, then the number of strings in the system is constrained to be even. The 1D quantum model is therefore restricted to the even-parity fermion subsector. If the 2D model is instead defined on a cylinder with the periodic direction parallel to the strings, then this restriction is lifted.
One of the utilities of the 2D classical to 1D quantum mapping is that for more complicated 2D models with longer range interactions it provides a good starting point for phenomenologial theories.

D J 1A -J 1B model with an unconstrained manifold
The next model we consider is the TLIAF with anisotropic nearest-neighbour interactions but now with defect triangles allowed (i.e. in the unconstrained manifold). The point is to better understand the crossover between the spin liquid and the weakly-correlated paramagnet and the crossover between Ising and Pokrovsky-Talapov criticality, which are both also features of the dipolar TLIAF. The Hamiltonian H ABB [Eq. 32] is the same as in Appendix C, except for the important difference that the manifold of Ising configurations is unconstrained, meaning defect triangles are allowed.

D.1 Dimer mapping
H ABB [Eq. 32] with an unconstrained manifold can be mapped onto a dimer model on the honeycomb lattice, but there is no longer a hardcore constraint, since vertices at the centre of defect triangles are covered by 3 dimers. Since the Grassmann path integral approach to determining the partition function requires the dimers to obey a hardcore constraint, it is necessary to instead consider the mapping onto a dimer model on the extended honeycomb lattice (described in Appendix B.1 and Fig. 13). This type of mapping was suggested in a more general context in [56], and makes possible an exact evaluation of the partition function.
The dimers can be categorised as those covering A, B or C bonds of the original triangular lattice (see Fig. 1 for bond labelling) or "extra" dimers, covering the bonds introduced in the act of extending the honeycomb lattice. The total number of dimers is fixed and given by, where N bond = 3N refers to the number of bonds of the original triangular lattice and N ext dim is the number of dimers on "extra" bonds. The energy of a given configuration relative to that of the ground state can be written as, and therefore, where the sum is over all dimer coverings of the extended honeycomb lattice, and the factor z 2N A ensures that Z exhon is equal to Z hon [Eq. 36] in the limit where z A → 0 and z B → 0 with z finite (i.e. the condition for being in the constrained manifold of Ising configurations).

D.2 Evaluation of the partition function
The evaluation of the partition function proceeds as in Appendix C, with the main difference being that there are 6 rather than 2 lattice sites in the unit cell (see also Ref. [60] for slightly different way of evaluating the partition function).    Fig. 16), and these are taken to be unit length. (b) The spectrum | k | [Eq. 60] along the path k = (k, k + π) for J 1A = 1 and J 1B = 1.5. For T < Tc (black) the spectrum is gapped at all k, and this corresponds to the stripe-ordered phase. At T = Tc (red) the gap closes at k = (π, 0) and an Ising transition occurs. For T > Tc (blue) the gap reopens. For Tc < T < T Is the minimum of | k | is at k = (π, 0). For T > T Is the minimum migrates away from k = (π, 0). In the limit T → ∞ the spectrum becomes flat.
Each site of the extended honeycomb lattice is assigned a real Grassmann variable, as shown in Fig. 19, and these are labelled a l i and b l i where i labels the unit cell, and l ∈ {1, 2, 3} labels sites within the unit cell. Bond weights are z on B and C bonds, 1 on A bonds and z −1 A on "extra" bonds, in accordance with Z exhon [Eq. 55]. The partition function is given by, where, This can be diagonalised by taking the Fourier Transform of the Grassmann variables (see Eq. 40), resulting in, After rewriting the action as a matrix equation, taking the Pfaffian of the matrix and absorbing the z 2N A factor, one can show that, It can be seen that this reduces to the | k | of the constrained manifold [Eq. 42] when the limit z A → 0 and z B → 0 is taken such that z remains finite (equivalently J 1A → ∞ and J 1B → ∞ while δJ remains finite). It can also easily be checked that in the T → ∞ limit the entropy per site is S/N = log 2 as expected.

D.3 Physical properties
The spectrum | k | [Eq. 60], which is is shown in Fig. 19, determines the physical properties of H ABB [Eq. 32]. As in the case of the constrained manifold there is a phase transition between an ordered and disordered phase. However, we label the transition temperature T c rather than T K , since it is not technically a Kasteleyn transition, as will be explained below.
For T < T c the spectrum is gapped at all k, and this corresponds to the stripe-ordered phase. The main difference from the case of the constrained manifold is that local fluctuations involving the creation of pairs of defect triangles are possible, though, depending on the value of T c , they can be highly suppressed.
There is a phase transition at T = T c associated with the closing of the gap in | k |, and this occurs at k = (π, 0). It can be seen from Eq. 60 that this requires, and the solution of this equation gives the critical temperature.
For T > T c a gap reopens in | k |, and this signifies that correlations are exponential in the paramagnetic state [47].
In order to investigate the nature of the phase transition, it is natural to define a second temperature, T Is , such that in the temperature range T c < T < T Is the minimum of | k | is at k = (π, 0), while for T > T Is the minimum of | k | is at a temperature-dependent incommensurate wavevector. T Is can be determined from the equation, and it can be seen from Eq. 60 that for T c < T < T Is the gap is given by, min | k | = |1−2z−z 2 B |. After setting T = T c + δT , with δT T c and δT < T Is − T c , one can show that the gap goes as min | k | ∝ δT . Taking the correlation length to be inversely proportional to the gap, ξ ∝ 1/ min | k |, results in ξ ∝ δT −ν with ν = 1, and this is typical of a 2D Ising transition [43]. Therefore Ising critical exponents are realised in the temperature window T c < T < T Is . However, the caveat to this is that the Ising temperature window can be exponentially small, and this is the case for δJ For T > T Is the minimum of | k | moves away from k = (π, 0) and the critical behaviour crosses over to that of the Pokrovsky-Talapov universality class for δT T Is −T c . Thus in the situation where δJ J 1A the transition is technically an Ising transition, but all practical measurements, whether in experiment or simulation, will show the features of a Kasteleyn transition. The values of T c and T Is are shown as a function of δJ/J 1A in Fig. 22, and it can be seen that the Ising temperature window only starts to be significant for δJ/J 1A 0.3.
Further increases in T increase the size of the gap and in the limit T → ∞ the spectrum, | k |, becomes completely flat, corresponding to an uncorrelated paramagnet where all configurations are equally likely. The density of strings, n string , can be calculated using Eq. 45 and the result is shown in Fig. 20. In the stripe-ordered phase n string is low, but not fixed to zero, as it is possible to create bound pairs of defect triangles, connected by a pair of strings. Its value increases rapidly at T = T c , since the defect triangles unbind, and therefore strings can wind the system. On further increasing T the density of strings passes through n string = 2/3 (the value realised in the constrained manifold) before saturating at n string = 3/4.
Since n string is not zero in the stripe phase, it is not, strictly speaking, an order parameter. However, it remains a useful indicator of where the transition occurs, since the derivative dn string /dT diverges logarithmically, as can be seen in Fig. 20.

D.4 Ising to Pokrovsky-Talapov crossover
The nearest-neighbour TLIAF provides a good setting in which to study the crossover from Ising to Pokrovsky-Talapov critical behaviour, since physical quantities can be calculated directly in the thermodynamic limit. For T c < T < T Is the system shows Ising critical exponents, while for T − T c T Is − T c it shows Pokrovsky-Talapov criticality. The crossover between these two limiting cases can be understood by studying the density of defect triangles, n def (see Ref. [51] for a similar analysis in terms of monopoles in spin ice), and we postulate a scaling ansatz, where z def = exp[−E def /T ], E def = 2J 1B = 2(J 1A + δJ) and g φ is an unknown function. The exponent α is the usual heat capacity exponent, and is expected to take the value α = 1/2 [48], while φ is the crossover exponent. By calculating n def in the thermodynamic limit and performing scaling according to Eq. 63 we find a convincing data collapse for φ = 1, as shown in Fig. 21.

D.5 Phase diagram and correlations
The phase diagram for H ABB [Eq. 32] in the unconstrained manifold can be calculated exactly, and is shown in Fig. 22. The nature of the correlations can be explored via the spin structure factor [Eq. 4], and this is shown in the same figure for a representative set of parameters. The phase diagram shows three regions, a stripe-ordered phase, a strongly-correlated spinliquid region and a weakly correlated paramagnet. The stripe-ordered phase is separated from the disordered region by a phase transition at T c [Eq. 61], while we take the crossover between the spin-liquid and paramagnetic regions to occur when the density of defect triangles, n def reaches 10% of its saturation value (i.e. n def = 0.025). As δJ is increased, the transition temperature T c increases faster than the crossover temperature, and therefore the spin-liquid region shrinks.
The nature of the correlations in the disordered phase changes significantly with varying T and δJ, and this can be seen from studying the structure factor [Eq. 4]. S(r) can be calculated in the thermodynamic limit via the Grassmann path integral approach (see Appendix F), and some examples are shown in Fig. 22. Also shown is S(q), which for simplicity is calculated using Monte Carlo simulation. We make the ansatz that in the disordered regions S(r) takes the long-distance asymptotic form [5] (see Appendix F), where in the case of δJ = 0 the correlation length perpendicular to the strings, ξ ⊥ , can be different from that parallel to the strings, ξ . This is found to give good fits to the calculated values of S(r) after taking into account the definition of long distance given in Appendix C.4. In the spin-liquid region the correlation length is considerably larger than the lattice spacing, and the system approximately realises the algebraically decaying correlation function studied in Appendix C.4 for the constrained manifold. In particular for δJ = 0 and T J 1A the correlation length diverges as ξ ⊥ = ξ ∝ exp[2J 1A /T ] [53]. At the crossover to the paramagnetic region, the correlation length is approximately ξ ⊥ ∼ 5, with ξ ≥ ξ ⊥ . Since the density of defect triangles is by definition low within the spin-liquid region, most of the strings wind the system, and therefore the relationship q ≈ ±q string (T ) = (±πn string (T ), 2π/ √ 3) holds to a good approximation.
In the paramagnetic region the correlation lengths become comparable with the lattice spacing, and the structure factor has a very different form to the algebraic decay found for the constrained manifold. In this region the strings mostly form short closed loops, and therefore the relationship between q and n string breaks down. For δJ 0.5 there is a direct transition from the stripe-ordered phase to the paramagnet. In the Ising critical region close to the transition the correlation function shows the usual 2D Ising scaling and is peaked in reciprocal space at q stripe .
In the stripe-ordered phase the asympotic form of S(r) given in Eq. 64 is no longer relevant, and Bragg peaks form in S(q) at the ordering vector q stripe = (0, 2π/ √ 3). Fluctuations around the ground state are not strictly forbidden, but are rare unless T ∼ J 1A , which is only possible for large anisotropies. Eq. 32] in the unconstrained manifold can be exactly mapped onto a 1D quantum model of spinless fermions, as was found to be the case for the constrained manifold in Appendix C.5. The main difference is that in the unconstrained manifold defect triangles act as sources and sinks of pairs of strings. In consequence, pair creation and annihilation terms appear in the 1D quantum model.

D.6 Mapping to 1D quantum model
Following a similar logic to that of Appendix G, there is an exact mapping of H ABB [Eq. 32] onto, where, and the evolution of these parameters with the temperature of the classical model is shown in Fig. 23. where, Physical properties of the classical TLIAF can be calculated directly from the quantum model. For example the classical quantity n string [Eq. 45] is equal to the fermion density [43]. The parameters of the 1D quantum model depend on those of the classical model according to Eq. 66, and the relationship is shown for J 1B /J 1A = 1.5. As T → 0 then t → 0, µ → −1, ∆ → 0, µ/2t → −∞ and ∆/t → 0. The phase transition occurs when µ/2t = −1. In the limit T → ∞ then t → 2, µ → 0, ∆ → 2, µ/2t → 0 and ∆/t → 1.
The last simplified model we study consists of a TLIAF with first and second-neighbour interactions and a constraint forbidding defect triangles. The motivation is that this is the simplest form of further-neighbour interactions, and can be used to study a number of features of the more general TLIAF in a simplified setting. In particular we will consider the crossover of the phase transition into the stripe-ordered phase from first to second order via a Pokrovsky-Talapov tricritical point, the string Luttinger liquid (as opposed to the free-fermion spin liquid studied in Appendix C) and its crossover into a domain-wall network state. Since the further-neighbour interactions destroy the mapping onto a free-fermion model, we rely on a combination of Monte Carlo and perturbation theory.
The Hamiltonian is given by, where the second-neighbour bonds are labelled ij 2 and we consider the constrained manifold of Ising configurations (i.e. no defect triangles).

E.1 General considerations
Before turning to detailed calculations, it is worth considering some of the qualitative features of H ABB2 [Eq. 69], both in terms of the nature of the phase transition and of the correlations in the spin-liquid phase (there is no paramagnetic region due being in the constrained manifold). The second-neighbour interaction, J 2 , and the nearest-neighbour anisotropy, δJ, act in concert with one another, in the sense that they both favour a stripe-ordered ground state. However, they act in opposition in the sense that J 2 favours a first-order phase transition, while δJ favours a second-order transition.
This can be seen by comparing the δJ = 0 case to that with δJ J 2 . At δJ = 0 the J 2 interaction selects a 6-fold degenerate, stripe-ordered ground state from the manifold of constrained Ising configurations [10]. The J 1 -J 2 TLIAF has been extensively studied, both analytically and by Monte Carlo simulation, and it is known that there is a first-order phase transition into the stripe phase [7,8,10,9,11]. In the limit of J 1 → ∞ the transition occurs at T 1 = 6.39J 2 [11]. Therefore we expect that, in the region where J 2 δJ, the transition between the paramagnet and stripe-ordered state will be first order.
In contrast, the first-neighbour anisotropy, δJ, favours a 2-fold degenerate stripe-ordered ground state, with stripes running parallel to A bonds (see Fig. 1 for the definition of bond directions). For δJ J 2 the J 2 interaction is irrelevant, and to a good approximation the analysis of Appendix C applies, indicating that the transition is second order. One focus here will be to study the crossover between the first and second-order phase transitions, which occurs when J 2 and δJ are comparable in magnitude.
When the transition is second order it is driven by the creation of isolated strings that wind the system (in [10] this is discussed in terms of the closely related concept of double domain walls). In order for a second-order transition to occur it is necessary that there is a repulsive interaction between these strings, and this repulsion is entropically driven and associated with the no-crossing constraint [39,40]. We show below that further-neighbour interactions result in an energetically-driven attraction between the strings, and that the second to first order crossover occurs when this balances the entropically-driven repulsion.
The free energy of an isolated string can be calculated exactly, and this can be used to find the exact transition temperature in the case of a second-order transition. Relative to the ground-state energy, strings cost an energy per unit length of E string = 2δJ + 4J 2 and corners, at which the string changes direction, have an energy cost E c = 2J 2 [10,11]. It follows that the free energy per unit length of an isolated string relative to the stripe ground state is given by [10], The second order transition temperature, T K , can be calculated from solving the equation f string (T K ) = 0, and in the case of J 2 = 0 it can be seen that this reduces to Eq. 44. The behaviour in the spin-liquid state should be closely related to the nature of the phase transition, since it is also sensitive to whether the interaction between strings is attractive or repulsive. In the introduction it was argued that the associated fermionic model can be weakly or strongly coupled, and it makes intuitive sense that weakly-coupled fermions correspond to repulsive string-string interactions, while strongly-coupled fermions correspond to attractive string-string interactions. In the weak-coupling case it can be expected that the spin liquid realises a 2D classical equivalent of a Luttinger liquid. In the strong coupling case it is less clear what to expect a priori. The crossover between weak and strong coupling is controlled by the ratio J 2 /T , with weak coupling for T J 2 .

E.2 Diagrammatic perturbation theory
The first approximate method we use to better understand H ABB2 [Eq. 69] in the constrained manifold is that of perturbation theory around the high-temperature limit. This approach cannot hope to compete with Monte Carlo simulations in terms of quantitative measures of, for example, the transition temperature, but does provide useful physical insights that are not apparent in Monte Carlo. While the approach is well motivated in the "weak-coupling" regime, we find that it also gives some clues as to how the system crosses over to the "strongcoupling" regime and to the appearance of a first-order phase transition.
and this can be done using a Grassmann path integral approach, following in spirit Ref. [61]. The first step is to map the Ising model, H ABB2 [Eq. 69], onto a dimer model on the dual honeycomb lattice. For the nearest-neighbour interactions the mapping is the same as in Appendix C. The second-neighbour coupling maps onto dimer-dimer interactions, where dimers on the same hexagon interact if they are separated by one unfilled bond (see Fig. 24). It follows that the partition function can be written as, where N 2 is the number of dimer-dimer interactions (see Fig. 24). The mapping of Z hon2 [Eq. 72] onto a Grassmann path integral does not result in a purely quadratic action, and therefore it is not exactly solvable by this method. Instead the mapping results in an action including terms with 2, 4, 6 . . . 2N Grassmann variables, and one can write, where the quadratic term S 2 [a, b] is given in Eq. 39 and contains products of 2 Grassmann variables, the quartic term, S 4 [a, b], contains products of 4 Grassmann variables and similarly for higher order terms, with 2N the number of honeycomb lattice sites. For a particular dimer configuration, one can ask which terms in the expansion of the action are required to correctly assign the weight. The answer depends on the size of the largest cluster of dimers connected by pairwise interactions (see Fig. 24). If the largest cluster contains n dimers, then it is necessary to consider the terms S 2m [a, b] with m ≤ n. Since clusters that include a sizeable fraction of all the dimers are common, many dimer configurations require one to consider terms up to n ∼ N .
For an infinite lattice it is necessary to truncate the expansion of the action in order to be able to perform calculations. This can be done systematically by considering z 2 − 1 to be a small parameter, which is valid for T 2J 2 . The reason that this is a useful expansion parameter is due to the fact that S 2n [a, b] has a lowest order contribution proportional to (z 2 − 1) n−1 . Thus for a chosen value of n, it is only necessary to consider terms in the action up to S 2n [a, b]. A simple worked example on a finite-size lattice is given in Appendix H to show how this type of expansion works in detail. Here we will consider n = 2, and therefore only retain the S 2 [a, b] and S 4 [a, b] terms in the action, thus working at first order in the small parameter |z 2 − 1|.
The quartic term in the action can be determined by observing that for a 2-site unit cell there are 6 terms containing 4 Grassmann variables, and these are shown schematically in Fig. 25. Thus one finds, and taking the Fourier transform using Eq. 40 results in, where, is the interaction vertex symmetrised over the pairs {k 1 , k 3 } and {k 2 , k 4 } and, The truncated action, which is given by the sum of S 2 [a, b] [Eq. 41] and S 4 [a, b] [Eq. 75], has a quartic interaction term, and therefore it is not possible to perform the path integral exactly. Instead a perturbative diagrammatic approach can be used, as is standard in quantum field theory [61]. It is important to note that the expansion order of the perturbation theory is set by the truncation of the action, and only diagrams consistent with this order should be considered.
The first step in the construction of a diagrammatic perturbation theory is the calculation of the free Green's function, and this is given by, This can be used to perturbatively construct the interacting Green's function, which is given by, where˜ k = k + Σ k and Σ k is the self energy.
In the case we are considering, the anomalous Green's functions a k 1 a k 2 and b k 1 b k 2 vanish at all orders of perturbation theory, and this is related to the absence of defect triangles. In consequence the effective quadratic action takes the simple form, and it follows that the partition function can be written as, In order to be consistent with the expansion of the action to first order in the small parameter |z 2 − 1|, we consider the Hartree-Fock diagrams, and therefore approximate the self energy as, where, At this level of approximation the partition function is given by, . On the right-hand plot the spin-liquid is split into a string Luttinger liquid for T > T tri = 9.14J2 and a domainwall network for T < T tri (separated by dashed yellow line). (Middle) Structure factor S(q) calculated by Monte Carlo simulation of an L = 72 hexagonal cluster (letters correspond to those on the phase diagram). (Bottom) Cuts through both S(q) and S(r). (a) For δJ = 0 and close to the first-order phase transition S(q) has significant spectral weight around the perimeter of the triangular-lattice Brillouin zone, as is typical of a domainwall network. (b) On increasing the temperature spectral weight rapidly accumulates at q = (±2π/3, 2π/ √ 3), as is typical for a string Luttinger liquid. In the whole string Luttinger liquid region the asymptotic form of S(r) follows Eq. 89 with a parameter-dependent Luttinger parameter, K ≥ 1. (c) At temperatures just above the Pokrovsky-Talapov tricritical point, there is a near-degeneracy between string sectors and the structure factor therefore shows extended spectral weight in the qx direction. (d) Further increasing the temperature breaks this quasi-degeneracy and sharp peaks form at q string (T ) = (±πn string , 2π/ √ 3). (e) At temperatures just above the second-order transition the structure factor is sharply peaked at q string (T ). (f) For T J2 the behaviour of the nearest-neighbour TLIAF in the constrained manifold is recovered, with K = 1.
The effective action S 2,eff [a, b] [Eq. 80] can be used to study the physical properties of the system, as in Appendix C. In particular we focus on the crossover between a second and first-order phase transition, which corresponds to the crossover from the weak to the strong coupling regimes (in fermionic language). Information about the nature of the phase transition can be extracted from the spectrum, HF k [Eq. 84], and it can be seen in Fig. 25 that this undergoes a change of structure at δJ/J 2 = 0.56.
For δJ/J 2 > 0.56 the spectrum, HF k [Eq. 84], shows the characteristic features of a second-order transition. In the disordered phase it has a gapless point at a temperature-dependent and incommensurate wavevector. As the temperature is lowered towards the critical point the gapless point migrates towards the wavevector k = (π, 0), and the critical temperature can be found from solving the equation HF (π,0) (T ) = 0. Below the transition the spectrum is gapped at all wavevectors, and the minimum is at k = (π, 0).
The second-order transition temperature is known exactly from Eq. 70, and Fig. 26 shows a comparison between the exact value and the estimate from first-order perturbation theory. First-order perturbation theory seems to work well even approaching the tricritical point, where T tri ≈ 9J 2 (the tricritical temperature will be determined more accurately by Monte Carlo simulations in the next section). At this temperature the small parameter is 1 − z 2 (T tri ) ≈ 0.2, and so the perturbation expansion is reasonably well controlled. The discrepancy in the critical temperature between zeroth and first-order perturbation theory can be seen from expanding the exact second-order transition temperature as, where it can be seen that for δJ ≈ J 2 the J 2 /δJ term is larger than the leading term. In fact further expansion of the transition temperature reveals that at δJ ≈ J 2 higher order terms are not small, but do cancel one another. However, it is important to remember that J 2 /δJ is not the expansion parameter. Exactly at the critical temperature the spectrum has qualitatively the same behaviour as the J 2 = 0 case (see Appendix C) close to the gapless point. Along the path k = (k, k + π) the spectrum grows as (k − π) 2 . This behaviour is typical of a Pokrovsky-Talapov transition [39,40].
At δJ/J 2 = 0.56 the spectrum shows a change of character. The coefficient in front of the quadratic term goes to zero, and the spectrum grows as (k − π) 4 around the gapless point. We refer to this point as a Pokrovsky-Talapov tricritical point, as the critical exponents are different from the standard ones of the Kasteleyn transition. This change of behaviour is not just an artifact of first-order perturbation theory, since its effects can be observed in Monte Carlo simulation (albeit at δJ/J 2 = 0.7 -see Appendix E.3).

δJ=0
T/J 2 =6.3 T/J 2 =6.72 T/J 2 =8 Figure 27: The spectrum HF k [Eq. 84] at δJ = 0 and for varying temperature. The path through 2D reciprocal space is parametrised by k = (k, k + π). In the paramagnet (red) there is a single gapless point in the region k > 0 and this occurs at an incommensurate wavevector. At T /J2 = 6.72 (orange) the gap at k = (π, 0) closes, but there remains a gapless point at an incommensurate wavevector. At lower T (blue) the gapless points approach one another, and the spectrum is very flat in their vicinity. While it is clear that the perturbative approach has broken down at such a small value of δJ/J2, the results are suggestive that lines of zeros appear in the spectrum, and this would be consistent with a first-order phase transition.
For δJ/J 2 < 0.56 the perturbative approach breaks down, but can be used to find some clues as to the true situation. In the paramagnet there is a gapless point at an incommensurate wavevector, as shown for the case of δJ = 0 in Fig. 27. As the temperature is reduced this migrates towards k = (π, 0), as is the case for a second-order transition. However, before the gapless point reaches k = (π, 0) the gap at k = (π, 0) closes, resulting in a pair of gapless points. This situation is not physical, and does not obviously correspond to the expected first-order phase transition. However, it can be seen that the spectrum is very flat between the two gapless points. We suggest that in reality gapless lines should develop in this region, and this would correspond to a first-order phase transition. This type of behaviour can never be exactly recovered using a perturbative approach, since a gapless line relies on the correct relationship between all coefficents in the expansion of the free energy.

E.3 Phase diagram determined from Monte Carlo simulations
As a complement to the perturbation theory approach, we also study H ABB2 [Eq. 69] using Monte Carlo simulation.
The simulations are carried out using a worm algorithm very similar to that presented in Ref. [11]. This works in the dimer representation (see Appendix B.1), and creates loops of alternating dimer-filled and empty bonds, which are then flipped, resulting in the reversal of all the Ising spins contained within the loop. The loop creation is carefully controlled such that detailed balance is maintained, and the absence of rejection results in an efficient algorithm. Hexagonal shaped clusters with periodic boundary conditions are used, containing N = 3L 2 Ising spins, where L measures the length of one side. Simulations are performed using system sizes from L = 24 up to L = 192.
The phase diagram of H ABB2 [Eq. 69], as determined by Monte Carlo simulation, is shown in Fig. 26. The phase transitions can be located either from measuring the triangular average of the winding number and associated susceptibility, defined in Eq. 3, or by measuring the heat capacity, and the results are consistent.
In the region where the phase transition is second order, the critical temperature is found from finite-size scaling analysis. We use the standard relation for a Kasteleyn transition [62], where c is a constant, L is the linear dimension of the system and ν = 1 is the critical exponent of the correlation length in the direction parallel to the double domain walls, below which the algebraic scaling of spin correlations breaks down [40,48,47]. We consider the parallel correlation length, ν , rather than the perpendicular correlation length, ν ⊥ , due to the anisotropy of the system which results in ν ⊥ = 1/2 = ν . Since the clusters used in the simulations are hexagonal in shape, and therefore isotropic, the growth of correlations parallel to the strings dominates the finite size effects. The exact second-order transition temperature is known from solving Eq. 70, and it can be seen in Fig. 26 that the finite-size-scaled Monte Carlo results are in good agreement with this.
As an example of such data collapse one can consider H ABB2 [Eq. 69] in the constrained manifold. We set δJ/J 2 = 1.5, since this is far enough from the tricritical point that deviations from the Pokrovsky-Talapov universality class are expected to be negligible. The results are shown in Fig. 28, and a convincing data collapse is found for β = 0.47 ± 0.04 and ν = 1.05 ± 0.09, which is consistent with the expected β = 1/2 and ν = 1. The line of second order transitions ends at a Pokrovsky-Talapov tricritical point, which is found to be at δJ/J 2 = 0.7 and T = T tri = 9.14J 2 .
In order to test for the presence of a Pokrovsky-Talapov tricritical point in Monte Carlo simulations one can use the scaling hypothesis, If the data for different system sizes can be collapsed using β = 1/4, then this provides good evidence of the presence of a Pokrovsky-Talapov tricritical point. We apply this scaling hypothesis to H ABB2 [Eq. 69] in the constrained manifold in Fig. 29, and find that for δJ = 0.7 the data can be convincingly collapsed using β = 0.21 ± 0.04 and ν = 0.91 ± 0.25. The findings from Monte Carlo can be seen to be in reasonable agreement with first order perturbation theory (see Appendix E.2), where a Pokrovsky-Talapov tricritical point was found at δJ/J 2 = 0.56.
For δJ/J 2 < 0.7 the transition is first order, and it is typically possible to simulate large-enough systems that the finite-size effects are small. In consequence the transition temperatures plotted in Fig. 26 are taken from the largest simulated systems. For 0.5 < δJ/J 2 < 0.7 this is L = 192, while for δJ/J 2 < 0.5 it is sufficient to consider L = 48. In order to gain physical insight into the crossover between a second and first-order phase transition, which in the spin-liquid region corresponds to the crossover between weak and strong coupling, we perform Monte Carlo simulations in a reduced manifold of states. The number of strings is fixed to be two, and the idea is to study the interaction between a pair of strings.
Monte Carlo simulations are performed on a square cluster with periodic boundary conditions, linear dimension L and total number of sites L 2 . In the 2-string manifold, allowed Ising configurations are distinguished by their J 2 energy, but all have the same energy in terms of δJ, since there are a fixed number of dimers occupying B and C bonds. In consequence it is not necessary to vary δJ/J 2 , but only T /J 2 . This shows that the string-string interactions are independent of δJ, and therefore the temperature at which weak coupling crosses over to strong coupling is also δJ independent.
For each considered temperature we measure the average separation between the strings, taking into account the periodic boundary conditions, and this is given by, where x 1 and x 2 are the positions of the strings along the x axis at a given height. It can be seen in Fig. 30 that as the temperature is reduced there is a change in ∆x starting at about T = T tri . For T > T tri the strings repel one another, and ∆x/(L/2) ≈ 1/2. This repulsion is entropically driven, and is due to the no-crossing constraint obeyed by the strings, which reduces the available fluctuations of a string if it is in close proximity to another string. This type of pairwise repulsion is crucial for the existence of a second-order phase transition out of the stripe phase, since it limits the number of strings that are condensed into the system when the free energy of an isolated, f string [Eq. 70], goes to zero.
For T < T tri the pair of strings start to approach one another, showing that the energeticallydriven attractive interaction starts to dominate over the repulsive interaction. The strings gain some binding free energy by being, on average, proximate to one another, and this is consistent with the crossover from a second to a first-order phase transition at T = T tri seen in the full Monte-Carlo simulations (see Fig. 26). The lower the temperature the more tightly the strings bind, suggesting that the transition should become more first-order as the temperature is decreased, and this is also consistent with the full simulations.
Since the string-string interactions are independent of δJ, the spin liquid region should have attractive string-string interactions in the temperature window T 1 < T < T tri , and we will discuss the implications of this in the next section.

E.5 Mapping to 1D quantum model and correlations
The nature of the correlations in H ABB2 [Eq. 69] can be used to understand the behaviour of the spin liquid, and can be determined by combining Monte Carlo simulation of the spin structure factor with insights from fermionic mappings.
It is useful to first consider at a qualitative level how the mapping to a 1D quantum model of spinless fermions is altered by the further-neighbour interactions (see Appendix C.5 and Appendix G for the nearest-neighbour case). At the level of an isolated string the J 2 interaction both increases the internal energy, and adds an energy penalty to "corners" where the string changes direction [10,11]. In the fermion model this alters the values of µ and t and adds a history dependence to the motion of the fermion, such that the passage from the imaginary timestep τ to τ + ∆τ depends not only on the fermion configuration at τ but also on the configuration at τ − ∆τ .
A second effect of the J 2 coupling is to drive an attractive interaction between strings. When strings neighbour one another their J 2 energy is reduced, and therefore the fermionic model also has an attractive interaction of the form V (z 2 )c † i c i c † i+1 c i+1 . In the string picture this attractive interaction is energetically-driven and competes with the entropically-driven repulsive interaction arising from the string non-crossing contraint. In the fermionic language the entropic repulsion maps onto the Pauli exclusion principle, which is a property of free fermions, and therefore the fermionic model is always attractive.
One advantage of mapping onto a fermion model is that it is known that fermions with weak attractive interactions form a Luttinger liquid, with Luttinger parameter K > 1 (K = 1 for free fermions) [52]. We therefore make the ansatz that the spin structure factor in the 2D classical model takes the asymptotic form [48,52], where K > 1. This corresponds to a reciprocal space structure factor with algebraically sharp peaks at q = q string . The asymptotic form given in Eq. 89 can be tested against Monte Carlo simulations, and we find that it gives a good fit to the simulations for T T tri , and some examples are shown in Fig. 26. The value of K can be extracted from the fits to the simulations, and the result of doing this for T > T tri and δJ = 0 is shown in Fig. 31. It can be seen that close to T = T tri the Luttinger parameter, K, becomes significantly different from the free fermion case of K = 1, while in the limit T /J 2 → ∞ the free fermion case is recovered, corresponding to 1/ |r| spin correlations (see Eq. 50). As a result of these findings we label the region of the spin liquid with T > T tri as a string Luttinger liquid. At T = T tri the entropic repulsion and energetic attraction between strings becomes comparable (see Fig. 30) and the strings start to bind together. At this temperature the distribution of spectral weight in the structure factor starts to rearrange itself such that S(q) is no longer dominated by a single q value, and Eq. 89 is inapplicable. Instead the weight is distributed around the perimeter of the triangular-lattice Brillouin zone (see Fig. 26), and this is typical of a domain-wall network (see supplementary material of Ref. [11]). Neighbouring parallel strings form domains in which Ising stripes are parallel to either B or C bonds, while domains with stripes parallel to A bonds correspond to an absence of strings, and an example of this is shown in Fig. 31. We find that the spin-liquid region of H ABB2 [Eq. 69] is best described as a domain-wall network in the region T 1 < T < T tri , as shown in Fig. 26. The domain-wall network state can be thought of as being a fluctuating, phase-separated state, with a loose analogy to the clustering of holes in superconductors [54,55].
The more the energetically-driven attraction between strings dominates over the entropic repulsion, the more tightly bound the strings and the larger the average domain size. In the case of H ABB2 [Eq. 69] and for δJ = 0 a first-order phase transition into the stripe phase occurs while the average domain size is relatively small. The addition of a third-neighbour interaction with 0 < J 3 < J 2 /2 suppresses the transition temperature, and therefore allows the average domain size to become larger since the attractive interaction becomes more important at low temperature [10,11].
Increasing δJ causes domains with stripes parallel to A bonds to grow, which corresponds to decreasing the string density, n string . At the tricritical point the A-domains coalesce and cover the whole system and there is a continous transition into the stripe phase.

F The spin-spin correlation function for the nearest-neighbour TLIAF
Here we show how to calculate the real-space, spin-spin correlation function, S(r) [Eq. 4], for the nearest-neighbour TLIAF, working in both the constrained manifold (i.e. without defect triangles) and the full, unconstrained manifold. Integral expressions for the correlation function can be derived in the thermodynamic limit, and numerical evaluation results in exact results up to numerical error. In the isotropic case these calculations just show how to derive the long-established results of Ref. [4,5] within the Grassmann variable approach [59]. The point of showing the calculations here is that the Grassmann approach makes it simple to extend the old results to the case of anisotropic interactions, for which it is necessary to separately consider correlations parallel and perpendicular to the string direction.
In this Appendix we show the mechanical steps used to calculate the correlation functions, while a physical discussion of the results is given in Appendix C, Appendix D and the main text. We consider separation vectors, r = r j − r i , that are either perpendicular or parallel to the average direction of the strings, and label the corresponding correlation functions as S ⊥ (r) and S (r) (corresponding to theê x andê y direction in Fig. 32).

F.1 Spin-spin correlations in the constrained manifold
First we consider the nearest-neighbour TLIAF with a constrained manifold. The calculations are slightly simplified by using a unit cell that contains 2 triangular lattice sites and 4 honeycomb/brick lattice sites, as shown in Fig. 32 (as opposed to the minimal unit cell with 1 triangular and 2 honeycomb/brick sites used in Appendix C). The two spins contained within the ith unit cell are labelled σ 1,i and σ 2,i and the perpendicular and parallel spin-spin correlation functions are, S ⊥ (rê x ) = σ 1,i σ 1,i+êx , S (rê y ) = σ 2,i σ 2,i+êy .
(90) The unit cell also contains 4 Grassmann variables, labelled a 1,i , b 1,i , a 2,i and b 2,i . These can be used to determine the partition function as in Appendix C, resulting in, where the action is, Fourier transforming the Grassmann variables results in, and this is diagonalised to give,

F.1.1 Correlations perpendicular to the strings
In order to calculate the spin-spin correlation function, it is necessary to express products of spins in terms of Grassmann variables. Before considering the general case, it is useful to first consider a pair of spins, σ 1,i and σ 1,i+êx , separated by a single honeycomb/brick lattice bond (see Fig. 32). If this bond is covered by a dimer then the spins are equivalent and σ 1,i σ 1,i+êx = 1, while if it is not dimer-covered σ 1,i σ 1,i+êx = −1. The expectation value is therefore given by, where P dim b 1,i ;a 1,i is the probability of finding a dimer on the bond connecting the Grassmann variables b 1,i and a 1,i . In order to determine P dim b 1,i ;a 1,i one can calculate a reduced partition function in which the sites b 1,i and a 1,i are excluded. Exclusion of these sites effectively fixes a dimer on the bond between them, and therefore P dim b 1,i ;a 1,i is given by the ratio of the reduced partition function to the original partition function, Z hon [Eq. 94]. In order to exclude the two sites, it is simply necessary to place b 1,i and a 1,i inside the partition function integral, using the properties of Grassmann variables (a 2 = 0). In consequence one finds, and it is clear that P dim b 1,i ;a 1,i is just the thermodynamic average of b 1,i a 1,i . In consequence, The thermodynamic average of two Grassmann variables can be calculated using, b 1,i a 1,j = 2 N k b 1,−k a 1,k e ik·(r j −r i ) e iky/2 , b 1,−k a 1,k = e iky/2 In the isotropic case (z = 1) integration yields P dim b 1,i ;a 1,i = 1/3, as expected, and therefore σ 1,i σ 1,i+êx = −1/3.
More generally, the correlation between a pair of spins with a separation vector parallel toê x is given by, This can be expanded using Wick's theorem, and rewritten as the deteminant of an r × rdimensional Toeplitz matrix, resulting in, with components, In the thermodynamic limit the sum can be converted into an integral giving, where u = 2z 2 (1 − cos k x ). The integral over k y is given by,
It follows that, where, is the Fermi wavevector of the quantum model H 1D [Eq. 51]. It can be seen that (M ⊥ ) mn = (M ⊥ ) nm and thus the Toeplitz matrix is symmetric. It is also worth noting that the matrix elements could have been calculated by making use of the exact mapping onto the 1D quantum model given in Appendix C.5.

F.1.2 Correlations parallel to the strings
The correlation between spins parallel toê y can be calculated by an analagous method. The difference is that a pair of spins are separated by not one but two dimers (see Fig. 32). As such the correlation function is given by, (2za 2,i+lêy b 1,i+lêy − 1)(2zb 2,i+lêy a 1,i+lêy − 1) .
where the z's take into account the weights of the excluded dimers. Wick's theorem allows this to be rewritten as the determinant of a 2r × 2r-dimensional Toeplitz matrix, with components, where m, n ∈ {1 . . . r}. The matrix elements can be calculated from, where,

F.2 Spin-spin correlations in the unconstrained manifold
Calculation of the spin-spin correlation function in the nearest-neighbour TLIAF with an unconstrained manifold follows a very similar pattern to that of the constrained manifold. However it is complicated by having to work with 6 or 12 Grassmann variables in the unit cell, as well as the fact that the extended brick lattice is not bipartite. The two-spin unit cell is shown in Fig. 33 and contains 12 sites of the extended brick lattice, and therefore 12 Grassmann variables, which are labelled a 1 ...a 6 and b 1 ...b 6 . The partition function can be calculated as in the constrained case, and this results in, with, The spin-spin correlation function in the direction perpendicular to the strings (parallel toê x ) is given by, and as in the constrained manifold case the correlation function can be written as the determinant of an r × r-dimensional Toeplitz matrix, M ⊥ , with matrix elements, The correlation function in the direction parallel to the strings (parallel toê y ) is given by, The matrix elements of interest can be determined from, where b 2,−k a 2,k etc. are relatively simple to calculate within the Grassmann variable approach, but result in very length expressions. Finally, we note that when calculating the determinant of the Toeplitz matrices, there is typically a finely-tuned cancellation between different terms. In the calculations presented above, this is unproblematic, since the matrix elements can be computed exactly (at least up to numerical accuracy). However, this limits the utility of the Toeplitz matrix approach for models with further-neighbour interactions where only a perturbative solution of the Grassmann variable spectrum is available. Small, unavoidable errors in the calculation of the matrix elements quickly have a significant effect on the determinant, resulting in unphysical results. For this reason we do not use this method for the J 1A -J 1B -J 2 model considered in Appendix E, but instead rely on finite-size Monte Carlo simulations of the correlation function.

G Mapping between the 2D Kasteleyn partition function and a 1D fermionic coherent-state path integral
Here we demonstrate the correspondence between the 2D, classical, nearest-neighbour TLIAF, H ABB [Eq. 32], and a 1D quantum model of spinless fermions, H 1D [Eq. 51]. For simplicity, we work in the constrained manifold of Ising configurations, but an analagous calculation can be carried out in the unconstrained manifold. In particular, we show that the Kasteleyn formulation of the partition function naturally maps onto a quantum path integral written in terms of fermionic coherent states. This makes clear that the Grassmann variables introduced in the Kasteleyn method describe coherent states of strings. The standard way to map between a d dimensional classical theory and a d−1 dimensional quantum theory is by equating the transfer matrix and the Hamiltonian according to T ≡ e −δτ H , where δτ is a small step in imaginary time. The partition function can then be viewed either as a matrix product of transfer matrices, or as a quantum path integral. Detailed examples of how to carry out this procedure in the cases of spin-ice and the cubic dimer model are presented in Ref. [63,64].
In the case of the TLIAF, the transfer matrix acts on states of strings, as shown in Fig. 34, and a translation across 4 bonds is necessary before the lattice structure repeats. These string states can be reinterpreted as the fermionic state of a 1D quantum model at a given The transfer matrix translates strings (shown in purple) by four honeycomb bonds in the y direction (from one blue dashed line to the next). At each translation the string can follow one of four possible routes, with the result that its x coordinate either remains the same or gets translated one step to the left or right. The string state on one of the blue dashed lines can be specified by the presence (1) or absence (0) of a string at each site on the line. This can be re-interpreted as the presence or absence of a spinless fermion at a particular imaginary timestep in a 1D quantum model. imaginary time coordinate, and the classical partition function sum is therefore equivalent to the fermionic path integral.
While it is possible to solve the nearest-neighbour TLIAF using a transfer matrix approach [2,3], the solution is considerably more compact using the Kasteleyn formulation expressed as a multiple integral over Grassmann variables (see Appendix C and Appendix D). Furthermore, the Kasteleyn appoach provides a good starting point for perturbative studies of more complicated models (see Appendix E.2). As such it would be useful to know how to link the Kasteleyn action to that of the fermionic path integral.
We demonstrate below that the Kasteleyn action naturally maps onto the quantum action when written in terms of fermionic coherent states. To do this we first re-examine the classical partition function using a non-minimal, 4-site unit cell, motivated by the fact that the string states shown in Fig. 34 involve a translation across 4 sites. We then examine the coherentstate fermionic path integral, and show that the action can be brought to the same form as the Kasteleyn action by introducing and summing over extra degrees of freedom that take into account the intermediate sites present in the honeycomb/brick lattice (see Fig. 34). Finally we link the Kasteleyn spectrum to that of the quantum model.
It should be noted that the mapping could just as well have been performed in the other direction, by starting from the Kasteleyn action and performing a Gaussian integral over half of the Grassmann variables to arrive at the fermionic coherent-state path integral. While this alternative method is probably slightly more direct, we feel that the method we present makes clearer the physical link between the two.

G.1 Kasteleyn action in 4-site basis
The procedure for determining the Kasteleyn action in terms of Grassmann variables was developed in [6] and is reviewed in Appendix C. To ease the comparison with the 1D quantum model, H 1D [Eq. 51], we here re-determine the Kasteleyn action of the nearest-neighbour TLIAF in the constrained manifold, using a 4-site unit cell and a different set of bond orientations compared to the main text.
The reason for using a 4-site cell is that it is natural to identify a string traversing 4 sites of the honeycomb/brick lattice with a single imaginary timestep in the quantum model (see Fig. 34). The bond orientations are shown in Fig. 35 and the reason they are different from those in the main text is just to simplify the mapping. They are of course chosen in accordance with Kasteleyn's theorem [46] and therefore there is no effect on the physical properties.

G.3 Matching the fermionic and Kasteleyn actions
The action, S 1D [η,η], exactly reproduces the partition function of the nearest-neighbour TLIAF, but it does this by averaging over some of the degrees of freedom of the Kasteleyn action. In order to make the mapping explicit, it is necessary to introduce these extra degrees of freedom into the quantum path integral. The microscopic relation between the quantum and classical partition functions requires a correspondence between an imaginary timestep and a translation of the classical system across 4 bonds in the y direction (see Fig. 36). In the classical set-up the string can hop by a maximum of one 1D lattice site per imaginary timestep, and the expansion of the quantum time-translation operator can therefore be truncated to first order without approximation, giving However, it can be seen in Fig. 36 that, if the string configuration is only known every 4 bonds (i.e. on the blue dashed lines in Fig. 36), there is an ambiguity, since a string can take two possible routes that leave its x coordinate invariant. In terms of the original Ising spins, these two possible routes describe different configurations. In order to explicitly describe all where, S 1D [η,η, θ,θ] = l,m −η l,m η i,m −θ l,m θ l,m + z(η l,m+1 θ l,m +η l+1,m+1 θ l,m +θ l,m η l,m +θ l−1,m η l,m ) .

(140)
The fermionic action is now in a form that can be directly compared with the Kasteleyn action S 2 [a 1 , b 1 , a 2 , b 2 ] [Eq. 119]. It can be seen that these can be brought to the same form simply by identifying, η l,m → a 1,i ,η l,m → b 1,i , θ l,m → a 2,i ,θ l,m → b 2,i , where there is an equivalence between i, which labels the unit cells in the 2D lattice, and (l, m), which labels the sites and timeslices in the 1D quantum problem. This justifies the mapping between the quantum and classical coefficients given in Eq. 51.

G.4 Matching the classical and quantum spectrums
As well as showing how the Kasteleyn and fermionic coherent state actions can be brought to the same form, it is also useful to show the link between the spectrums. Since the fermions/strings are free, the partition functions can be simply evaluated by Fourier transform, and the spectrums compared. The Kasteleyn spectrum is given by | K p | [Eq. 123], and the partition function, Z hon [Eq. 121] is the product of this spectrum.
In conclusion, there is an exact mapping between the Kasteleyn action and that of the fermionic, coherent-state path integral with Hamiltonian H 1D [Eq. 124]. This mapping also makes it clear that the Grassmann variables introduced in the Kasteleyn formulation of the classical partition function describe coherent states of strings, and therefore demonstrates the link between the Kasteleyn and transfer matrix approach to solving the nearest-neighbour TLIAF.
H Perturbative expansion of the action: a simple example The Grassmann path integral representation of interacting dimer problems, as used in Appendix E.2, is not unknown [61,66] but has not been widely explored in the literature. As an aid to the interested reader, we here consider Z hon2 [Eq. 72], and show a worked example of how to evaluate this partition function via Grassmann path integration on the simplest, nontrivial lattice: the hexagonal plaquette. This provides useful insights into the construction of a perturbation theory for the infinite lattice, as presented in Appendix E.2. We consider the dimer covering of a hexagonal plaquette with a dimer weight of 1 on A bonds, z on B and C bonds and a dimer interaction with weight z 2 between dimers separated by one unfilled bond (see Fig. 37). For a single plaquette there are only two possible dimer coverings, shown in Fig. 37, and each of these has a weight z 2 z 3 2 . The partition function, Z hon2 [Eq. 72], is therefore given by, Z hon2 = 2z 2 z 3 2 = 2z 2 (z 2 − 1) 3 + 3(z 2 − 1) 2 + 3(z 2 − 1) + 1 , where the second equality is an exact rewriting that will prove useful below. While in such a simple case the partition function can be calculated exactly just by inspection, it is instructive to perform the calculation via the Grassmann path integral representation. On a finite lattice the highest order term in the action is S 2N [a, b], where 2N is the number of honeycomb lattice sites, and for the 6-site plaquette the partition function can therefore be rewritten as, where i = {1, 2, 3}. The quadratic term in the action does not take into account the z 2 interaction, and is given by, If the action is truncated at quadratic order, then the usual rules of Grassmann integration can be used to find, Comparison with the exact value of the partition function [Eq. 148] shows that this quadratic approximation becomes exact in the limit z 2 − 1 → 0. The quartic term in the action takes into account pairwise interactions of the dimers in isolation from other pairwise interactions, and is given by, The factor z 2 − 1 is chosen such that if there were a configuration with only 1 dimer-dimer interaction, the -1 would remove the contribution from the purely quadratic action, while the z 2 would replace this with a contribution that takes the interaction into account. In the case of the hexagonal plaquette the allowed dimer configurations contain 3 mutually interacting dimers, and this mutual interaction is not fully taken into account by the quartic term. Direct evaluation results in, reproducing the exact partition function [Eq. 148] to first order in z 2 − 1. Finally, the hexatic term takes into account the fact that the dimers are not interacting in isolation, but are all mutually interacting, and is given by, where in the first line the −1 removes the contribution from the quadratic action, the 3(z 2 −1) removes the contribution trom the quartic action and the z 3 2 replaces these with a contribution that correctly reproduces the weight of three mutually interacting dimers. Direct calculation including quadratic, quartic and hexatic terms correctly reproduces Eq. 148 for the partition function, i da i db i e S 2 [a,b]+S 4 [a,b]+S 6 [a,b] = 2z 2 (z 2 − 1) 3 + 3(z 2 − 1) 2 + 3(z 2 − 1) + 1 = 2z 2 z 3 2 .
As the lattice size is increased, it rapidly becomes impossible to determine the partition function by inspection. A full expansion of the partition function in terms of Grassmann variables also becomes complicated due to the increase in the number of terms in the action. However, the advantage of this method is that it provides a way of systematically carrying out perturbation theory around the non-interacting limit |z 2 −1| → 0. For large or infinite lattices direct evaluation of Grassmann actions with quartic and higher order interacting terms is not possible, but approximate diagrammatic methods can be used, and the order of expansion matched to that of the truncation of the action.