Existence of N\'eel order in the S=1 bilinear-biquadratic Heisenberg model via random loops

We consider the general spin-1 SU(2) invariant Heisenberg model with a two-body interaction. A random loop model is introduced and relations to quantum spin systems is proved. Using this relation it is shown that for dimensions 3 and above N\'eel order occurs for a large range of values of the relative strength of the bilinear ($-J_1$) and biquadratic ($-J_2$) interaction terms. The proof uses the method of reflection positivity and infrared bounds. Links between spin correlations and loop correlations are proved.


Historical Setting
In this work properties of the spin-1 Heisenberg model are deduced using a random loop model first introduced in the work of Nachtergaele [18]. Random loop models have been around since the work of Tóth [21] and Aizenman and Nachtergaele [1]. In [21] a lower bound was obtained on the pressure of the spin- 1 2 Heisenberg ferromagnet; this improved the bound of Conlon and Solovej [5]. Sharp bounds have recently been found [6]. The loop model presented in [1] applies to the spin- 1 2 Heisenberg antiferromagnet. Both spin models can be applied to higher spins; for a review of these models we refer, for example, to [14]. The work of Ueltschi [23] combined and extended these loop models. It has recently seen attention for its usefulness in several aspects of quantum spin systems. In [23] it is shown that there is long-range order in various spin systems, including nematic order in the spin-1 system. The work of Crawford, Ng and Starr [7] on emptiness formation also makes use of the model, as does the work of Björnberg and Ueltschi [4] on decay of correlations in the presence of a transverse magnetic field. The loop model presented here comes from [18], it is similar in flavour to the Aizenman-Nachtergaele-Tóth-Ueltschi representation. See Refs. [1,21,23,18,19] and references therein.
Quantum spin systems are currently a very active area of research. The growth of popularity of probabilistic representations has allowed new methods to be applied to these systems with many interesting results. This work looks at the general SU (2) invariant spin-1 Heisenberg model with a two-body interaction (1.1) Here we will have x ∈ Λ ⊂ Z d and E the set of nearest neighbour edges. The operators S = (S 1 , S 2 , S 3 ) are the spin-1 matrices, see section 4 for details of the model. The work in [23] shows that in the region 0 ≤ J 1 ≤ 1 2 J 2 the system exhibits nematic order in the thermodynamic limit if the temperature is low enough and the dimension is high enough. Nematic order was also shown independently using different methods in [20]. It is also shown that if Λ is bipartite there will be Néel order for J 1 = 0 ≤ J 2 at low temperature. This corresponds to the occurrence of infinite loops in the related loop model. Alternatively in d ≤ 2 infinite loops should not occur, it is proved in [13] that this is the case for J 2 = 0, the extension to J 2 > 0 should be straightforward. The first proof of continuous symmetry breaking was shown by Fröhlich, Simon and Spencer [12] for the classical Heisenberg ferromagnet (and hence antiferromagnet). This result was extended by Dyson, Lieb and Simon [8] to the quantum antiferromagnet. The result excluded the case d = 3 and S = 1 2 , it was extended to this case in the work of Kennedy, Lieb and Shastry [15]. These works all used the method of reflection positivity and infrared bounds. For information on reflection positivity see Refs. [2,3,10,11] and references therein. The Heisenberg ferromagnet is not reflection positive and hence does not benefit from these methods.

Main result
In this article we use the method of reflection positivity and infrared bounds on a random loop model. Links between correlations in the spin model and probabilities of events in the loop model are also derived in section 5. We focus on the quadrant J 1 ≤ 0 ≤ J 2 , see Fig. 1 (in fact we need only work on the unit circle J 2 1 + J 2 2 = 1 but the quadrant is more convenient pictorially). The following result concerning Néel order follows from Proposition 5.3 (a), Theorem 6.1 and the discussion that follows. For the precise statements see Section 6.
Theorem. For Λ ⊂ Z d a box of even side length and d ≥ 3 there exists α = α(d) > 0 such that for J 1 ≤ 0 ≤ J 2 if −J 1 /J 2 > α then there exists c = c(α, d) > 0 such that Note that the "liminf" version of this theorem will be proved. The limits exist but this will not be proved here. It is shown in the discussion after Theorem 6.1 that this sum is positive if  have rigorous proofs of the relevant phases. The line J 1 < 0, J 2 = 0 is the Heisenberg antiferromagnet where antiferromagnetic order has been proven [8], Néel order extends into the dark yellow region. The dark blue region 0 ≤ J 1 ≤ 1 2 J 2 has nematic order at low temperatures [23], with Néel order on the line J 2 > 0, J 1 = 0. The adjacent dark yellow region has been proved to exhibit nematic order in high enough dimension [16]. Antiferromagnet order is expected here but is not yet proved. It can be shown [8,15] that I d → 0 and K d → 1 as d → ∞ and that both are decreasing in d. This means we can prove that the region where Néel order occurs will increase to the entire quadrant J 1 ≤ 0 ≤ J 2 as d → ∞ i.e., the ratio α(d) is decreasing. In d = 3 there is Néel order at low temperature in the spin system for −J 1 /J 2 < 0.46, this is a triangular region of angle 65 • measured from the J 1 axis.
Reflection positivity for this quadrant is already known; for J 1 < 0 = J 2 it was shown in [8] and for J 1 = 0 < J 2 one can see, for example, [16] lemma 3.4 for an explicit proof. It was proved in [8] that Néel order occurs for J 1 < 0 = J 2 , it is clear the result extends to a neighbourhood of the axis J 1 < 0 < J 2 with J 2 sufficiently small. However it is impossible to extend the result concerning Néel order any significant amount without some new results. This is where the loop model has been essential. Indeed in [8] an infrared bound is obtained of the form where (A, B) Duh is the Duhamel correlation function and ε(k) = 2 d i=1 (1−cos k i ) for k ∈ Λ * . Notice that this bound becomes weaker as |J 1 | decreases (equivalently on the unit Figure 2: Events of the process ρ J 1 ,J 2 , a) represents single bars, b) represents double bars and c) represents the uniform measure on vertical segments being either parallel or crossing.
circle as |J 1 |/|J 2 | decreases). Transferring this bound to S 3 0 S 3 x β (k) requires the Falk-Bruch inequality, which would involve dealing with the term After some calculation one obtains correlations in Proposition 5.3 such as S 1 x S 1 y S 3 x S 3 y β . Hence to work directly in the quantum system using the methods of [8] one must obtain good bounds on these correlations. Simple bounds such as taking the operator norm are not sufficient due to the weakening of (1.4) as |J 1 |/|J 2 | decreases. Without using the loop model it is not clear how to obtain such bounds currently. The random loop model is presented in sections 2 and 3. The spin-1 Heisenberg model is introduced in Section 4. In Section 5 the connection between the loop model and the quantum system is proved. In particular it is shown how to write various correlation functions in terms of probabilities of events in the loop model; some of these correlations are also presented in [22]. In Section 6 the main result concerning Néel order is presented and proved.

The random loop model
We now introduce the loop model presented in [18]. To begin we take a finite set of vertices, Λ, with a set of edges, E ⊂ {{x, y}|x, y ∈ Λ, x y}. We associate to this lattice a new lattice,Λ, and edge set,Ẽ: There are two lattice sites inΛ for every site in Λ and four edge inẼ for each edge in E. We will write x 0 , x 1 in place of (x, 0), (x, 1).
For β > 0 consider a process, ρ J 1 ,J 2 , consisting of a Poisson point process on E × [0, β] and a uniform measure on segments of Λ × [0, β] between events of the Poisson point process. The Poisson point process has two events that we will refer to as 'single bars' and 'double bars'. Note that this process is on the edge set E, the events define corresponding events on the edge setẼ. The single bars will occur at rate −2J 1 and double bars at rate J 2 for J 1 ≤ 0 ≤ J 2 . The rate for the single bars is written in this way to be consistent with the connection to the quantum spin system that will be introduced in Section 4. The interval [0, β] will be referred to as a time interval. The uniform measure is on two possibilities, "crossing" and "parallel". How to build loops from these events is described in detail below, see Fig. 2 for pictorial representations of the events. We first define the single and double bars. Single bars occur at a point (x, y, t) for {x, y} ∈ E. We define the corresponding geometric event onẼ as a bar joining x 1 and y 0 at time t. Double bars occur at a point (x, y, t) and the corresponding event oñ E is a bar joining x 1 and y 0 and a bar joining x 0 and y 1 , both at time t. A loop of length l is then a map γ : [0, βl] per →Λ × [0, β] per such that γ(s) γ(t) if s t, γ is piecewise differentiable with derivate ±1 where it exists. If s is a point of nondifferentiability then {γ(s−), γ(s+)} ∈Ẽ. Loops with the same support and different parameterisations are identified. For a realisation ω of ρ J 1 ,J 2 we associate a set of loops as follows: Starting at a point (x i , s) ∈Λ × [0, β] we move upwards (i.e. in direction of increasing s). If a bar is met at time t it is crossed and we then continue in the opposite direction from (y j , t), where y j is the other site associated to the bar. Each maximal vertical segment between bars (x 0 , x 1 ) × [s, t] (i.e. bars involving the site x occur at times s and t and no bar involving x occurs for u such that s < u < t) is either parallel (nothing happens) or crossing (the sites x 0 and x 1 are exchanged). If time β is reached the periodic time conditions mean we continue in the same direction starting from time 0. We denote by L(ω) the set of all loops associated to a realisation ω. Loops are most easily understood pictorially, see Fig. 3. Note that the loops could be defined via a Poisson point process onẼ × [0, β] where bars can occur between x i and y j with each (i, j) being equally likely. However one would still need to introduce the crossing or parallel events so that it is still possible to have x 0 and x 1 in the same loop even when there is no bar occurring on any edge containing x.

For this loop model we have partition function
Here θ > 0 is a parameter and ρ J 1 ,J 2 is the probability measure corresponding to a Poisson point process of intensity −2J 1 for single bars and J 2 for double bars. The relevant probability measure is then We are interested in sets of realisations, ω, where certain points ofΛ × [0, β] are in the same loop. Probabilities of these events are connected to correlations in the spin-1 quantum system presented in Section 4, they will be required in the proof of Néel order in Section 6. Particular events of interest will be denoted pictorially , see Fig. 4. These events are defined and denoted as follows.
a) The event that sites x i and y j are connected (in the same loop). Note that the probability of x i and y j being connected is independent of i and j. Denoted E[x i y j ].
b) The event that x 0 and x 1 are connected, y 0 and y 1 are connected but there is no connection from any x i to any y j . Denoted E x 0 c) The event that x 0 and y 0 are connected, x 1 and y 1 are connected but x 0 and x 1 are not . We can also have x 0 and y 1 connected and x 1 and y 0 connected but x 0 and x 1 not connected and denote the event in the analogous way. These events both have the same probability.
d) The event that all four sites The definition of bars means that if a loop is followed starting from a point x i ∈Λ (by moving in either the up or down direction) then the direction it is travelling upon arriving at a point y j ∈Λ in the same loop is determined only by the number of bars the loop has encountered between the sites. For example on a bipartite lattice defined by sublattices Λ A and Λ B such that {x, y} ∈ E ⇐⇒ x ∈ Λ A , y ∈ Λ B the direction that x i is left and y j is entered will be the same if x and y are in the same sublattice and different if they are in different sublattices.
Sometimes the order in which sites are encountered along the loop will be important. In this case arrows will indicate the order that sites will be encountered in on following the loop (up to parameterisation). The events E x 0 T c , means that upon leaving site x 0 if we encounter y 0 before encountering x 1 then we will then encounter y 1 and then x 1 before closing the loop. As this notation is potentially confusing (but also seemingly unavoidable) the reader will be told explicitly when the order is important. When wanting the probability of these events we will drop the E from the notation, as below.
It is intuitively clear that P(x 0 y 0 ) decays exponentially fast with respect to x − y for β small. Hence P

Space-time spin configurations
In order to make the connection with spin systems we need the notion of a space-time spin configuration. The spin system we shall connect to is the spin-1 Heisenberg model, we shall make this connection via an intertwining that merges two spin-1 2 models. For this reason we will take θ = 2 from Section 2 (2S + 1 for S = 1 2 ). This is also the reason the latticeΛ has two sites for every site in Λ. It is also possible to represent the spin-S model for general S by merging 2S spin-1 2 models, this will meanΛ will have 2S sites for every site in Λ. See [18] for more details. This generalisation together with some results analogous to the ones presented here should be straightforward once the spin-1 model is understood. It is not immediately clear which results will still hold however, investigation is required.
From now on we take the cubic lattice in Z d with side length L, denoted Λ L , with periodic boundary conditions. The edge set, E L , will consist of pairs of nearest neighbour lattice points. Precisely Where x − y is the graph distance between x and y. A space-time spin configuration is a function σ x i ,t is piecewise constant in t for any x i . We further define Σ to be the set of all such functions with a finite number of discontinuities. For a realisation of the process ω we consider σ that are constant on the vertical segments of each loop in L(ω) and that change value on crossing a bar. This restriction on configurations will allow to make the link with spin systems. We call such configurations compatible with ω and denote by Σ (1) (ω) the set of all compatible configurations. The following relation holds as we work on a bipartite lattice, meaning fixing a configuration's value at some (x i , t) determines the configuration on the entire loop containing (x i , t): from which we can obtain We further define the set Σ (1) x i ,y j (ω) ⊃ Σ (1) (ω) to be compatible configurations along with configurations that flip spin at points (x i , 0) or (y j , 0) (or both) but are otherwise compatible.
In the remainder of this section it will be convenient to ignore the condition that compatible configurations flip value on crossing a bars as it adds unnecessary extra complication. Hence we further define Σ (2) (ω) to be configurations that are constant on loops (and hence do not flip value at bars). Σ (2) x i ,y j (ω) ⊃ Σ (2) (ω) denotes the set of configurations in Σ (2) (ω) along with configurations that flip spin at points (x i , 0) or (y j , 0) (or both) but are otherwise consistent with the definition of Σ (2) (ω). The reader should bear in mind in the sequel that the connection with the spin-1 quantum system requires spin flips at bars, with this in mind the modifications to the remainder of this section required to incorporate the spin flips are easy.
As in [23] we will later need a more general setting for the measure on space-time spin configurations. We consider a Poisson point process onẼ × [0, β] with events being specifications of the local spin configuration. We will consider discontinuities involving two pairs of sites (x 0 , x 1 , y 0 , y 1 ). The objects of the process will be a set of allowed configurations at these sites immediately before and after t. We can denote these events as Implicit here is an ordering on Λ with x < y. An event A is a subset of {−1/2, 1/2} 8 and occurs with intensity ι(A). More precisely we let ι : and σ x i ,t is otherwise constant in t. The measure is then given by ρ ι with the counting measure on compatible configurations. We note that different intensities can give the same measure as in [23], for ι and ι intensities it is shown in [23] that We want to write the Poisson point process involving bars in terms of intensities of specifications of spins. We require that specifications corresponding to single and double sets of bars have intensity −2J 1 and J 2 respectively. If we naively defineι by For any a, a , b, c, c ∈ {1/2, −1/2}. Now each specification is disjoint from the other two and single and double sets of bars have intensities −2J 1 and J 2 respectively, as required. We also have ι(A) = 0 for any other specification. Then This representation can be used to show reflection positivity of the loop model, however in light of the famous paper of Dyson, Lieb and Simon [8] showing reflection positivity for J 1 < 0 = J 2 and the recent paper [16] showing reflection positivity for J 1 = 0 < J 2 we have reflection positivity for J 1 < 0 < J 2 more directly, as we shall see.

The general spin-1 SU(2) invariant Heisenberg model
Let S ∈ 1 2 N. A spin-S model has local Hilbert spaces H x = C 2S +1 . Observables are Hermitian matrices built from linear combinations of tensor products of operators on ⊗ x∈Λ H x . Physically important observables are expressed in terms of spin matrices S 1 , S 2 and S 3 , operators on C 2S +1 that are the generators of a (2S +1)-dimensional irreducible unitary representation of su (2). They satisfy the commutation relations where α, β, γ ∈ {1, 2, 3} and E αβγ is the Levi-Civita symbol. Denote S = (S 1 , S 2 , S 3 ), then S · S = S (S + 1)1. The case S = 1 2 gives the Pauli spin matrices. For S = 1 there are several choices for spin matrices, we will use the following matrices: Consider a pair (Λ, E) of a lattice, Λ ⊂ Z d , and a set of edges, E, between points in Λ.
We will take Λ to be a box in Z d , hence Λ is bipartite. We denote by Λ A and Λ B the two disjoint lattices such that Λ A ∪ Λ B = Λ and every e ∈ E contains precisely one site from Λ A and one site from Λ B .
Then we take the operator S i x for i = 1, 2, 3 to be shorthand for the operator S i x ⊗ Id Λ\{x} . Recall the definition ofΛ andẼ above, we shall use these below.
The most general SU(2) invariant Hamiltonian with two-body interactions for spin-1 is We will soon drop the parameters J 1 , J 2 from H J 1 ,J 2 Λ for readability. In this article we will be concerned with the region where J 1 ≤ 0 ≤ J 2 . Associated to this Hamiltonian we have the following partition function and Gibbs states for β > 0: Again we shall drop the parameters J 1 , J 2 from the notation.
The following new definitions come from Nachtergaele [18]. We introduce an isometry V : C 3 → C 2 ⊗ C 2 with the property V D 1 (g) = (D 1 2 (g)) ⊗2 V for g ∈ S U(2) and D S the spin-S representation of SU (2). Here the representation D 1 is given by the matrices (4.2) and D 1 2 is given by the Pauli matrices. it is clear such an isometry exists as we can define it for spin matrices and then extend by linearity (recall that the spin matrices generate the representation). From this we obtain the key relation where σ i are the spin- 1 2 matrices (hence 2σ i are the Pauli matrices). Further we have where P is the projection onto the spin triplet. Hence VS i acts on C 2 ⊗ C 2 and so using the notation before V x S i x acts on ⊗ x∈Λ C 2 ⊗ C 2 . We make the following definition One can check that R i = (σ i ⊗ 1 + 1 ⊗ σ i ). To make expressions more concise we will also denote A X := ⊗ x∈X A x for X ⊂ Λ. For these new operators we have a new Hamiltonian (note we have now dropped the J 1 and J 2 parameters) and associated Gibbs states The connection with the previous Gibbs state can easily be made explicit, Λ,β . (4.12) We use Dirac notation in the following way: |a, b denotes an element of the one site Hilbert space C 2 ⊗ C 2 and |a, b ⊗ |c, d for two sites etc.
There are two operators of particular interest, both act on two sites. Firstly we define S (1) by its matrix elements Geometrically this requires spin b and c and the spins b and c to be the negative of each other and also requires a = a and d = d . This corresponds to the the single bars in the loop picture. The second operator, D (1) , is also defined via its matrix elements 14) The geometrical interpretation this time is that of the double bars. The actual operators needed are S (1) = PS (1) P and D (1) = PD (1) P in order to account for bars occurring between any x i and y j with each i and j from {0, 1} being equally likely. Note here that from this definition we see that we require the spin value to change sign on crossing a bar as was mentioned in Section 3. There are also extra factors in S (1) and D (1) of e iπa for the bottom half of a bar (denoted ) and e −iπa for the top half of a bar (denoted ) where a = ± 1 2 is the spin value on the site in Λ A associated to the bar. By direct computation of the matrix elements we can prove the relations Using these relations we can rewrite the Hamiltonian in the region J 1 ≤ 0 ≤ J 2 as

The random loop representation
We can neglect the term (J 1 + J 2 )P x,y in the Hamiltonian (4.17) and instead add (2J 1 − J 2 )1, this does not change the Gibbs states. Doing this allows to use a useful lemma from [1] exp Here ρ is the measure associated to a Poisson point process on E×[0, 1] with two events occurring with intensities u and ν respectively. The product is ordered according to the times at which the events occur. C is either A or B depending on which event occurs. This is actually a slight extension of the lemma presented in [1]. From this we can obtain here each A (1) is one of S (1) or D (1) . The process has intensity −2J 1 for single bars and J 2 for double bars. Again the product is ordered by the time events occur.
We now prove the connection between the loop model and the quantum system. This will enable us to understand certain important correlation functions. After this we should have the tools we need to calculate any two point correlation (at least ones involving only spin operators). The first thing to understand is the extra factor, which we shall denote by z x i ,y j (σ, ω), the product of all factors e ±iπa from operators S (1) and D (1) corresponding to the bars in loop(s) containing x i and y j in a realisation ω of ρ J 1 ,J 2 . Again a ∈ {1/2, −1/2} is the value that σ assigns to the portions of these loop(s) in the Λ A sublattice (or if all of a loop is on the sublattice Λ B a is given by the negative of the value assigned to the loop). The value of z x i ,y j (σ, ω) is specified by the following lemma: Before the proof we should note that the lemma says that the only dependence on σ is at x i and y j at time zero. If the spin does not flip at both sites that we get total factor 1, else it depends on which sublattices the sites are in. If the spin only flips at one site then there are no compatible configurations hence the value of the total extra factor is unimportant.
Proof. To begin note that we can take (i, j) = (0, 0). The result for (i, j) (0, 0) follows as the choice of i or j does not affect which sublattice the two sites are in. Suppose σ ∈ Σ (1) (ω). Moving upwards from x 0 the first bar encountered is , the bars encountered then alternate between and . Moving downwards from x 0 we first encounter a bar then alternate between and . This means we can make a matching If σ ∈ Σ (1) x 0 ,y 0 (ω) \ Σ (1) (ω) and (−1) x−y = 1 and ω ∈ E[x 0 y 0 ], then x 0 and y 0 are in the same sublattice. We can thus deduce that the section of loop that moves upwards/downwards from x 0 crosses an even number of bars before reaching y 0 . This means that the loop containing x 0 and y 0 contains an even number of bars of each type ( or ). Hence we can make a matching of a bar in one 'half' of the loop with a bar in the other 'half' and the same with bars , with some bars left over. The factors from bars in the matching will thus be 1 as the spin flip at x 0 at time 0 means one bar in each pair has factor e ±iπa and one bar has factor e ±iπ(−a) . Here by 'half' of a loop we mean the section that connects x 0 at time 0+ with y 0 at time 0− or x 0 at time 0− with y 0 at time 0+. There are still possibly some bars left over as each half of the loop may have a different number of bars in it. A moments thought reveals that there must be an even number of bars left, half of type and half of type . As the bars have factor e −iπ(±a) and the bars have factor e iπ(±a) we have full cancellation again and have total factor 1.
For the remaining case σ ∈ Σ (1) x 0 ,y 0 (ω) \ Σ (1) (ω) and (−1) x−y = −1 and ω ∈ E[x 0 y 0 ], we have x 0 and y 0 in different sublattices. We can see as last time that the factors from the 'extra bars' (that arise from each half of the loop having a different number of bars) will cancel as again there are equal numbers of and . For the remaining bars there are an odd number in each half of the loop, this means we can make a matching for all but two of the bars. The factors from bars in the matching will cancel each other. For the remaining two bars one is a with factor e iπ(±a) and one is a with factor e −iπ(∓a) (the sign of a is opposite due to the spin flip at x 0 at time 0). This means the overall factor is (±i) 2 = −1. This completes the proof.
In light of Proposition 5.1 the following proposition can be proved in the same way as Theorem 3.2 in [23].
Proposition 5.2. The partition functions Z (1) Λ,β are given by We also have the following identity where σ z i is the value of a space time configuration, σ, at time 0 and site z i .
With the important details understood we can calculate some correlations in terms of probabilities in the loop model. The most important correlations here are the Néel and nematic correlations (Proposition 5.3 (a) and (b) respectively).
Proof. We will calculate the correlations in order. First note that each S i plays an equivalent role, hence cyclic permutations of the indices (1, 2, 3) does not alter the expectation. Using this together with (S i S j ) T = ±(S j S i ) (the sign depending on the value of i and j) means we can take i = 3 and j = 1. For each we will expand using (4.8) and (4.12).
Proof of (a). First Λ,β . (5.6) We see that due to sites z 0 and z 1 being interchangeable for z ∈ Λ each of the four terms in the sum have the same expectation. We also know from Proposition 5.2 T r(σ 3 ⊗ 1) x (σ 3 ⊗ 1) y e −βHΛ = ρ(dω) We note that the integral differs from zero only on the set where x 0 and y 0 are connected. If x and y are in different sublattices the product of spin configuration values is − 1 4 , if in the same sublattice the product is 1 4 . We can deduce that Proof of (b). For the second correlation (R 3 We see that expanding as before gives (5.10) From this and the fact that (S 3 x ) 2 Λ,β = 1 3 S x · S x Λ,β = 2 3 we can deduce that For the first term in the correlation we again note that (S 3 x ) 2 (S 3 y ) 2 Λ,β = (R 3 x ) 2 (R 3 y ) 2 1 Λ,β . We then calculate as before: x,y . (5.12) Now following through the same expansion as before we have (5.13) Using (5.11) and noting that the last term in the sum requires either two loops containing two of the sites x 0 , x 1 , y 0 , y 1 each or one loop containing all four sites to give a non-zero contribution to the sum overall (if one site is not connected to any other its spin value can be ± 1 2 independently of other sites, averaging the integral on this set to zero) we have (5.14) The probability P x 0 x 1 y 0 y 1 comes with twice the weight because there are two ways to connect both sites at x to different sites at y (but only one way both sites at x can be connected and both sites at y can be connected).
Proof of (c). For the third correlation we use the same expansion The factor 4 has come from grouping together terms such as σ 1 ⊗ σ 3 and σ 3 ⊗ σ 1 that have the same expectation. A useful observation at this stage is that σ 1 σ 3 = −i 2 σ 2 . Calculating further and noting that the two cross terms in the above product have the same expectation we see (5.16) From the symmetric roles of σ i for i = 1, 2, 3 and part a) we know the first term is − (−1) x−y 4 P(x 0 y 0 ). For the second term we need σ 2 ⊗ 1 ⊗ σ 1 ⊗ σ 3 (1) Λ,β . This is the expectation of a matrix with purely imaginary entries, due to the one appearance of σ 2 . Now we note three pieces of information that allow us to calculate this expectation. All the matrices e −βH (1) Λ , P Λ , σ 1 , σ 2 and σ 3 are Hermitian. The matrices σ i are acting on different sites inΛ and hence they commute. e −βH (1) Λ and P Λ commute and have real entries. This means taking the adjoint of the operator leaves the expectation unchanged. Because the operator is purely imaginary we should obtain the negative of what we started with on taking the adjoint. Hence the correlation must be zero.
For the last term we expand as in Proposition 5.2 and obtain σ∈Σ (1) x 0 ,y 0 (ω) z x 0 ,y 0 (σ, ω) σ ·,0+ |σ 1 ⊗σ 3 ⊗σ 1 ⊗σ 3 |σ ·,0− (5.17) Here σ ·,0± denotes the full spin configuration for some σ ∈ Σ x 0 ,y 0 (ω) at time 0± respectively. Also note that as σ 1 flips spins and σ 3 does not the set of space-time spin configurations Σ (1) x 0 ,y 0 (ω) is the correct set. We could expand the set of configurations we sum over to include configurations that flip spin at sites x 1 and y 1 at time zero but these would not be compatible with σ 3 acting at time zero at those sites hence they would not contribute. Recall that a loop that contains a site that spin flips at time zero cannot contain only one such site, hence the set of configurations that contribute to the integral is E[x 0 y 0 ]. Again the set of configurations where one of the sites x 1 or y 1 is not connected to any of the other three does not contribute to the integral. Combining these two facts we see that the only sets of configurations that contribute to the integral are those where there are two loops each containing two sites (one with x 0 and y 0 and the other with x 1 and y 1 ), or one loop containing all four sites. For the case of two loops there is one factor of z x 0 ,y 0 (σ, ω) = (−1) x−y from the loop containing x 0 and y 0 (where σ 1 acts). Another factor of (−1) x−y comes from the loop containing x 1 and y 1 and the condition that the spin flips on crossing a bar. Note that for the first loop there is no such factor coming from spin flips at bars because σ 1 | ± 1 2 = + 1 2 | ∓ 1 2 hence there is a factor of + 1 2 regardless of the spin value at the site. For the case of one loop containing all sites the order that sites occur in the loop is important; this is because both σ 1 and σ 3 are acting at sites in the loop. If, when following the loop, the site y 1 appears directly before or after the site x 1 then the section of loop between these sites follows the normal rule of flipping spins at bars (or if we follow the loop the other way we pass through two spin flips at time zero as well, these cancel each other out as far as the product of spins at sites x 1 and y 1 is concerned). This means we have a factor of (−1) x−y as before. If one of the sites x 0 or y 0 appears between sites x 1 and y 1 on the loop the effect of the extra spin flip changes the sign of the factor coming from the product of spins, giving a factor of −(−1) x−y . As before we also have the factor z x 0 ,y 0 (σ, ω) = (−1) x−y in both cases. This means the correlation is Recall that the arrows in the events show the direction that the loop is traversed. From this we can finally deduce that Now we note that the last two probabilities are equal (swap y 0 and y 1 ).
The correlations (d) and (e) follow easily using the same techniques and considerations as above.
From this we can easily obtain some bounds on these correlations that are potentially very difficult to obtain without the loop model.
Proof. All inequalities are immediate from Proposition 5.3 when we note that E x 0 Other inequalities of interest involve correlations between nearest neighbour points. Equation (29) in [20] allows us to obtain the following bound in the ground state (β → ∞) Now looking at Proposition 5.3 (b) for x − y = 1 (say x = 0, y = e 1 ) we see that if J 1 = 0 then the event P(0 0 e 1 0 ) puts us into the case of one of the last two probabilities.
Ignoring the first probability (as it is difficult to control) we obtain (for J 2 > 0 = J 1 ) This bound is positive for d ≤ 4, however it was not sufficient to deduce nematic order from a theorem analogous to 6.1 but concerning the nematic correlation function.
6 Occurrence of Néel Order

Setting and results
We take the cubic lattice in Z d with side length L, denoted Λ L , with periodic boundary conditions. The edge set, E L , will consist of pairs of nearest neighbour lattice points. Precisely For the main theorem we need to introduce two integrals, they come about due to similar considerations as in [15] Here (·) + denotes the positive part and ε(k) = 2 d i=1 (1 − cos(k i )). Positivity of this lower bound implies Néel order for those values of J 1 and J 2 in the spin-1 system at low enough temperature. This result is also equivalent to the occurrence of macroscopic loops in the loop model.
Of course we see that for −J 1 , J 2 > 0 the positivity of the lower bound doesn't depend on the value of J 2 1 + J 2 2 , only on the ratio −J 1 /J 2 . This means there corresponds an angle, measured from the J 1 axis, such that for angles less than this we have proved the existence of Néel order/macroscopic loops. The bound is positive if One of these is certainly satisfied if I d K d < (−4J 1 )/(−J 1 + 4J 2 ). A table of values of I d and K d for various d is presented in [23]. If J 2 1 + J 2 2 = 1 this is the case in d = 3 for J 1 < −0.42, d = 4 for J 1 < −0.28 and d = 5 for J 1 < −0.22.
A similar theorem concerning nematic order (corresponding to correlation (b) in 5.3) can be proved using the same methods. Unfortunately showing that one of the lower bounds obtained was positive proved difficult due to the seemingly unavoidable issue of bounding more complicated connection probabilities from below.

Proof of Theorem 6.1
The result can be proved directly using the loop model. However in light of the proof of reflection positivity for the J 2 interaction [16], we can complete the proof more directly, appealing to Proposition 5.3 when required. To begin we define a Hamiltoniañ denote its Gibbs states · ∼ . Note that for h = 0 this Hamiltonian is unitarily equivalent to H J 1 ,J 2 Λ L as can be seen from defining U = x∈Λ B e iπS 2 x where Λ B ⊂ Λ L is the odd sublattice. It has been shown [8,16] that this Hamiltonian is reflection positive for and Z(v) = T re −βH(v) . (6.9) By the usual methods for reflection positivity, using [16] to handle the J 2 terms we can show that for any v we have Z(v) ≤ Z(0)e −J 1 β(v,∆v) . By choosing v according to v x = η cos(k · x) for η small we can recover the infrared bound in [8]: with y = 0 and y = e 1 . We see from using unitary operator U with (6.10) and (6.14) we have after taking limits the bounds in Theorem 6.1.