Symmetry Breaking and Convex Set Phase Diagrams for the q-state Potts Model

We demonstrate that the occurrence of symmetry breaking phase transitions together with the emergence of a local order parameter in classical statistical physics is a consequence of the geometrical structure of probability space. To this end we investigate convex sets generated by expectation values of certain observables with respect to all possible probability distributions of classical q-state spins on a two-dimensional lattice, for several values of q. The extreme points of these sets are then given by thermal Gibbs states of the classical q-state Potts model. As symmetry breaking phase transitions and the emergence of associated order parameters are signaled by the appearance ruled surfaces on these sets, this implies that symmetry breaking is ultimately a consequence of the geometrical structure of probability space. In particular we identify the different features arising for continuous and first order phase transitions and show how to obtain critical exponents and susceptibilities from the geometrical shape of the surface set. Such convex sets thus also constitute a novel and very intuitive way of constructing phase diagrams for many body systems, as all thermodynamically relevant quantities can be very naturally read off from these sets.


Introduction
In a series of ground breaking papers in the late 19 th century, Gibbs [1,2,3] elegantly derived the thermodynamic stable state of a given substance through the minimization of some thermodynamic potential (later known as the free energy), in fact by means of a geometric construction. In particular, Gibbs considered a surface given by the possible values of the thermodynamic extensive quantities (such as e.g. energy, volume and entropy) of a system of interest and realized that points on this surface with tangent planes of equal orientation correspond to possible stable states of the substance at a temperature and pressure given by the orientation of the tangent plane. If two (or more) points belong to the same tangent plane, the corresponding states can coexist in equilibrium, characteristic for first order phase transitions. If two or more points have tangent planes with equal orientation but different distance to the origin, the state whose tangent plane is closer to the origin is metastable, corresponding to a supercritical system [3].
This geometrical construction can be interpreted as identifying the thermodynamically stable states as the extreme points of a convex set consisting of all possible realizable values of the thermodynamic extensive quantities of a given system. In the case of Gibbs' construction the relevant convex set is essentially the convex hull of the thermodynamic surface, termed "secondary surface" by Maxwell (who also produced a plaster clay model of the surface for water as a present to Gibbs in 1874). All thermodynamic properties of a system of interest can then be read off from the geometric features of this set and phase transitions correspond to non-analyticities on the surface, which arise by considering convex hulls of analytic functions [4].
In this paper we extend previous work [5] and construct in full detail convex set thermodynamic surfaces for a paradigmatic model of classical statistical mechanics on a lattice, namely the Ashkin-Teller-Potts model [6,7,8,9]. As in the case of Gibbs' original thermodynamic surface, the extensive quantities of the system are in competition with each other and stable states, which constitute the thermodynamic surface, are again those that minimize the free energy. While Gibbs' original construction is capable of detecting regions of phase coexistence at first order phase transitions, they however show no signatures at second order phase transitions, as there the thermodynamic extensive quantities vary continuously across the critical point. In [5] it is demonstrated that by including the order parameter corresponding to such a phase transition as an extensive quantity into these sets, phase transitions are signaled through the appearance of characteristic geometrical features in the form of ruled surfaces. As these sets exist as a collection of all possible realizable states of a given system without any prior reference to any Hamiltonian which generates dynamics, the reason for the occurrence of symmetry breaking phase transitions thus lies in the geometry of the space of all possible realizable states.
Similar convex set pictures have been studied in the context of the N -representability problem in quantum chemistry [10,11,12,13], but without including order parameters. In the following we will construct convex set thermodynamic surfaces for the qstate Potts model and study its geometrical features. In section 2 we construct and discuss these sets, in particular in section 2.1 we demonstrate how symmetry breaking leads to characteristic ruled surfaces and flat parts, which are a signature of symmetry breaking phase transitions. In sections 2.2 and 2.3 we describe additional features of the surface where the model at zero temperature becomes equivalent to coloring problems or hard-square lattice-gas models with nearest neighbor exclusion. We further describe in section 2.4 how to obtain critical exponents and susceptibilities from a given convex set surface. We conclude with final remarks and outlooks in section 3. We additionally give information about the tensor network representations used to obtain numerical data in Appendix A and show scatter plots generated by drawing random probability distributions of spin configurations in the Potts model in Appendix B.

The Potts Model and its Convex Set Representations
The q-state Potts model [6,7,8,9] is a generalization of the ubiquitous Z 2 -symmetric Ising model [14,15] to Z q -symmetry. It has been shown to be correspond to a Z q lattice gauge theory of matter [16,17] and in certain parameter regimes to coloring problems [18,19] and hard-square lattice-gas models with nearest neighbor exclusion (1NN) [20].
The Potts model in a magnetic field is defined by the Hamiltonian where z i = 1, . . . , q is a q-state spin on site i, ij denotes nearest neighbors and δ is the Kronecker delta function. We consider the model in two spatial dimensions on a square lattice. At zero field, where the model possesses Z q -symmetry, it undergoes a symmetry breaking phase transition at some finite critical inverse temperature β c = log( √ q + 1) [7,21], where for β > β c the Z q -symmetry is spontaneously broken. For q = 2 the Potts model is equivalent to the classical Ising model [14] and can thus be solved exactly in zero field for all temperatures [9,22]. For general q > 2 and zero field the model can be mapped onto a staggered six-vertex model, which can be solved exactly only at criticality [23,24]. Other solvable cases include J < 0 at T → 0 and zero field for q = 3 on the square lattice [18], and q = 4 on the hexagonal lattice as well as q = 3 on the Kagome lattice [19]. The symmetry breaking phase transition in zero field is continuous for q ≤ 4 and of first order for q > 4 [7]. The nature of the phase transition will become apparent from the geometrical features of the corresponding convex set phase diagrams which we construct below.
Consider the space of all possible probability distributions P (z) of configurations of q-state spins z i = 1, . . . , q with i the position on a two-dimensional square lattice with N sites, which form a convex set in some high-dimensional parameter space. In particular we consider three-dimensional projections of this set in the thermodynamic limit N → ∞, parameterized by the three observables nearest neighbor interaction energy per site shifted magnetization per site and entropy per site where . . . denotes expectation values with respect to P (z). The convex set C is then given by all possible points X = [ δ(z, z ) , z , s], such that δ(z, z ) , z and s are compatible with each other, i.e. they stem from a common valid probability distribution P (z). This is an instance of the classical marginal problem [25,26,27,28]. Notice that we are using a shifted magnetization with an offset q+1 2 , such that the convex set is reflection symmetric with respect to z . The extreme points on the surface of this set are then naturally given by Gibbs states of (1).
To see this, consider (hyper)planes in this three-dimensional parameter space, which are defined as families of points X ∈ C, related by a plane equation of the form n · X = n x δ(z, z ) + n y z + n z s = n d, where n is the normal vector of the plane and d is the distance of the hyperplane to the origin. Setting n x = 2J, n y = h and n z = T , this yields exactly the (negative of the) free energy per site of (1) where the factor 2 comes from the fact that every site has 4 nearest neighbors on a two-dimensional square lattice. ‡. For a given set of parameters (i.e. normal vector) the hyperplane tangent to the convex set has maximum possible distance from the origin and thus also minimizes the free energy, which is the definition of a Gibbs state. Every point on the surface thus corresponds to a state of thermodynamic equilibrium, at parameters given by the orientation of the tangent plane and free energy proportional to the distance of the tangent plane to the origin. Conversely, every point inside the convex set corresponds to a possible non-equilibrium state of the system.
If the tangent plane touches the convex set at a unique point only, then the thermodynamic stable state is unique and exactly given by a Gibbs state which yields the observables given by the tangent point for the parameters (J, h, T ) defined by the ‡ On a general isotropic lattice the free energy is given by −f = JK orientation of the tangent plane, i.e. its normal vector n. If however the tangent plane touches the set on an entire line or even a plane, then the state which minimizes the free energy for these parameters is not unique, which is a prerequisite of symmetry breaking. The set of valid states can then be parameterized by one (or more) real parameters. Such ruled surfaces (continuous sets of tangent lines) or planes are thus the geometrical signatures that will enable us to detect symmetry breaking and the emergence of a connected order parameter.
We show the surfaces of the resulting convex sets for the Potts model for q = 3 and q = 5 in figures 1 and 2 respectively (for the special case of the Ising model, corresponding to q = 2, see [5]). These sets show interesting geometrical features from which a wealth of other information, such as the nature of phase transitions, locations of critical points, critical exponents, susceptibilities, etc. can be extracted. The numerical data for plotting these surfaces has been obtained by means of tensor network techniques described in Appendix A. For scatter plots of points obtained from random probability distributions, which approximate the convex set from the inside, see Appendix B. Figure 1. Convex set generated by nearest-neighbor interaction energy δ(z, z ) , shifted magnetization z and entropy per site s of all possible probability distributions of 3-state spins on a two-dimensional square lattice. We plot the surface of this set, corresponding to Gibbs-states of (1) for q = 3. Due to reflection symmetry we only plot the upper half of the set. Blue lines denote points of constant J = ±1 and h and varying temperature T . The red line denotes the exactly solvable decoupled case J = 0 and thus separates regions of ferromagnetic and antiferromagnetic coupling. At the critical point A the emergence of a (green) ruled surface signals a non-uniqueness of the thermal equilibrium state at zero field and thus symmetry breaking. As a guide to the eye we have plotted a few vertical lines on the ruled surface, along which the tangent plane touches the convex set. Point B marks the end point of the bifurcation line of J = −1, h = 4 and T → 0, leading up to the (blue) top plane where the lowest energy state is exponentially degenerate, resulting in a finite residual entropy as described in section 2.2. A similar situation arises at point C, corresponding to the end point of the line J = −1, h = 0, T → 0. There again the lowest energy state is exponentially degenerate, resulting in a finite residual entropy as described in section 2.3. This plane is only present for q > 2 and does therefore not appear in the convex set drawn for the Ising model in [5]. As a guide to the eye we have drawn two-dimensional grids onto the top and left plane, emphasizing the fact that there the tangent plane touches the set on the entire respective planes.  figure 1 for the case of 5-state spins where the surface of this set is given by Gibbs-states of (1) for q = 5. For q > 4 the phase transition is of first order and thus comes with a discontinuity of the three observables at the critical point. This results in a coexistence region of the ordered and disordered phases and the critical point A gets stretched out into a (gray) flat triangular surface, where any mixture of the two phases is a valid state, i.e. the two phases coexist. This flat part then smoothly connects to the (symmetry broken) ordered phase represented by the green ruled surface. As a guide to the eye we have drawn a two-dimensional grid onto the flat triangular surface, emphasizing the fact that there the tangent plane touches the set on the entire triangular surface and we have also plotted a few vertical lines on the green ruled surface, along which the tangent plane touches the convex set. The flat surfaces emerging from points B and C are the same as described in figure 1.

Symmetry Breaking and the Ruled Surface
For zero field, J > 0 and T < T c the thermodynamic state that minimizes the free energy is q-fold degenerate and the Z q -symmetry can be spontaneously broken, such that z = 0. The maximum possible value z max can then be taken as the order parameter associated to this phase transition § . For a given set of parameters any state within this q-fold degenerate space thus minimizes the free energy and is characterized by the same values for δ(z, z ) and s, but different z This is nicely reflected in the convex sets through the emergence of a (green) ruled surface at the critical point. Zero field implies tangent planes with normal vectors lying in the z = 0 plane, i.e. n = [2J, 0, T ]. The tangent plane touches the convex set on a unique point in the z = 0 plane everywhere except for J > 0 and T < T c , where the tangent plane in fact touches the convex set along a whole line for each J and T , given by X(t) = [ δ(z, z ) , t z max , s] with t ∈ [−1, 1] and z max > 0 the maximum value of the order parameter. An infinitesimal value of h = 0 then immediately explicitly breaks the symmetry and causes the tangent plane to touch the set on a unique point of the set infinitesimally close to the edge of the ruled surface. Or equivalently, the curve of tangent points of a tangent plane given by n = [2J, h = 0, T ] as h → 0 ± will end in a point with z = ± z max = 0 for T < T c . This nicely reflects the fact that the order parameter can be obtained by first taking the thermodynamic limit at nonzero field before letting the field go to zero.
The nature of the phase transition changes from continuous to first order for q > 4, where a first order phase transition is characterized by a latent heat and a discontinuity of first derivatives of the free energy at the critical point. The internal energy and all other expectation values that can be written as a derivative of the free energy, such as the order parameter and also the entropy per site s therefore have a discontinuity at the critical point. In the convex set we can thus detect first order phase transitions through the appearance of flat hyperplanes at the boundary that arise even without additionally plotting the order parameter. At the critical point the thermal equilibrium state is not unique and any point on this hyperplane is a valid state of the system at the critical temperature. This corresponds to the coexistence of phases at the critical point which is characteristic for first order phase transitions. In the case of the Potts model, this flat hyperplane then smoothly connects to the ruled surface representing the symmetry § Given z max the shifted magnetization is then z = (k−(q+1)/2) z max with k = 1, . . . , q the integer enumerating the maximally symmetry broken states, characterized by one-site marginal distributions given by p(z) = 1/q + p(2δ z,k − 1), where p < 1/q is a function of T . Other order parameters for the Potts model have also been proposed. One possibility for defining an observable whose expectation value in the symmetry broken phase is independent of k is e.g. given by defining O(z) = exp(2πiz/q) and measuring | O | = pq ∈ [0, 1].
Mixtures of maximally symmetry broken states generally do not correspond to physically realizable states as they cannot be converted into each other by means of local modifications. Mathematically they are elements of disjoint Hilbert space sectors [29,30]. A hint towards this fact is given by the peculiar structure of the random scatter plots for quantum and classical systems is shown in section Appendix B. Figure 3. Surfaces of the two-dimensional convex sets generated by nearest neighbor interaction δ(z, z ) and entropy per site s for the zero field Potts model for q = 3 and q = 5. The phase transition is continuous for q = 3 and cannot be detected from the convex set without adding an additional axis corresponding to the order parameter z . For q = 5 the phase transition is however of first order and can thus be detected through the discontinuities of δ(z, z ) and s across the critical point A, which gets stretched into a straight (red) line where the two phases can coexist. As a guide to the eye we have extended this line to both sides to see that there is (albeit very small) curvature to both sides of the phase coexistence part.
broken phase T < T c (c.f. figure 2). For continuous phase transitions the thermodynamic state at the critical point is still unique and there is no such additional hyperplane. We can thus already detect first order phase transitions in the lower dimensional convex set that does not include the order parameter. In the case of the Potts model, a two-dimensional convex set parameterized by δ(z, z ) and s thus already suffices to detect the phase transition for q > 4, it will however show no signature of the phase transition for q ≤ 4 (see figure 3), for which adding an additional axis corresponding to the order parameter z is necessary.
We want to emphasize here that these convex sets and thus also the ruled surfaces exist prior to making any references to any model Hamiltonian, we just consider finite dimensional projections of the convex set of all possible probability distributions of a system of physical degrees of freedom. This means that the reason for the occurrence of symmetry breaking phase transitions ultimately lies in the geometrical structure of the space of all possible probability distributions. It would therefore be interesting to investigate all possible projections of this set and classify all possible ruled surfaces that can arise on such projections.

Top Plane
The top (blue) plane corresponds to parameters J = −1, h = 4 and T = 0 where the tangent plane touches the convex set on the entire top plane, meaning that the thermal Flip one spin from q to q-2 ΔE = +4 Flip two adjacent spins from q to q-1 Flip one spin from q to q-1 ΔE = 0 -4q + 4 Figure 4. Construction of the degenerate space of lowest energy configurations for the top plane, corresponding to J = −1, h = 4 and T = 0. Starting from the fully polarized state z i = q with lowest possible energy, flipping single spins from q to q − 1 leaves the overall energy invariant. Flipping two or more adjacent spins however results in a net energy increase, as does flipping from q to any z < q − 1. The resulting space consists of all configurations where z j ∈ [q, q−1] such that every z i = q−1 is completely surrounded by z j = q.
equilibrium state is not unique and in fact all states on this plane are valid equilibrium states for these parameters. Exactly at this point, the two terms in the Hamiltonian become "equally strong" in the following sense. If we start from the completely polarized state z j = q, the magnetic field term is minimized, whereas the interaction part has a positive energy contribution, resulting in a net energy of 2−4q per site. If we now flip one spin at an arbitrary position from q to q − 1, we gain exactly the same amount of energy from the interaction term as we lose from the magnetic field term and the overall energy stays the same. We can now continue flipping spins that way without changing the energy, as long as we never flip any spins next to an already flipped one, which would result in a net energy increase of +2. In general, a cluster of N f flipped spins and a boundary of length N b results in a net energy change of 4N f − N b ≥ 0, which is only zero for N f = 1. The two lowest energy states with the smallest magnetization are thus the two Néel states between q and q − 1. Similarly, flipping from q to any z < q − 1 always results in a net energy increase and the restricted space of lowest energy states is thus given by all configurations z j ∈ [q, q − 1] such that every z i = q − 1 is completely surrounded by z j = q (see also figure 4). This restricted space is equivalent to the configuration space for the nearest-neighbor exclusion lattice-gas model (1NN) [20] and grows exponentially with the system size.
At T = 0 all such configurations are equally likely; the entropy per site is therefore finite and measures the exponential growth of the space of lowest energy configurations. This symmetry of equal probability can however be spontaneously broken as any statistical mixture of such configurations is a valid state of the system with equal free energy f = 2 − 4q. The entirety of all such mixtures is exactly given by the top blue plane in the convex sets, where point B marks the state of equal probability which has maximal entropy.
To calculate the boundary of the top blue plane we consider tiny perturbations away from this point in parameter space, which immediately cause a jump onto the edge of the plane. Similar to degenerate perturbation theory we then simulate this perturbation Hamiltonian only within the restricted subspace of the top plane to lift the exponential degeneracy and determine its extreme points. The perturbation Hamiltonian is just the magnetic field term where µ is usually small. Since we however simulate this Hamiltonian in the restricted subspace only (which also makes the simulation non-trivial), µ need not be small and just controls the position along the edge of the top plane. We therefore wish to evaluate where the sum is only over the space of valid configurations Z t given by the top plane and µ ∈ R. The entropy per site s is then given by where z = Z 1/N is the partition function per site. The other observables δ(z, z ) and z are computed as usual but with respect to (8). Note that entropy and δ(z, z ) are independent of q and z for different q are related by just an offset. The top plane thus has the same shape for all q, but different vertical offset in z . Note that (8) is equivalent to the 1NN model in a chemical potential µ [31,32,33], where states q and q − 1 correspond to an empty and occupied site respectively. The limits µ → ±∞ correspond to the the completely polarized and the Néel states respectively (or equivalently the completely empty and maximally filled lattice respectively in terms of the 1NN model) and thus have zero entropy, while µ = 0 corresponds to point B with maximal residual entropy s res . Our calculated value at this point reproduces the (log of the) value κ(1) given in section 1.1 of [20] up to machine precision. The tensor network we use to simulate (8) is described in Appendix A.1.

Left Plane
The left (red) plane with δ(z, z ) = 0 corresponds to parameters J = −1, h = 0 and T = 0 where the tangent plane touches the convex set on the entire left plane, meaning that the thermal equilibrium state is not unique and in fact all states on this plane are valid equilibrium states for these parameters. For these parameters the lowest energy states are given by all configurations z j ∈ [1, q], such that no nearest neighbors are in the same state. This is the famous vertex coloring problem and consequently, the partition function can be written as a chromatic polynomial in q [34,35] and counts the number of proper vertex colorings of the two-dimensional square lattice with q colors. For any q > 2 the number of valid configurations is exponentially large in system size and we are thus presented with the same situation as for the top plane in the previous section, but with a different restricted subspace. For q = 2 (i.e. the Ising model) this problem is trivial as only two valid configurations exist (the two Néel states) and the left plane is absent).
Again, at T = 0 all these configurations are equally likely, leading to a residual, non-zero entropy s res . For q = 3 this can be mapped onto the problem of residual entropy of square ice [18], for which the value is known exactly as s res = 3/2 log(4/3) ≈ 0.431523 [36,37]. For q > 3 there are no exact solutions for the square lattice. The symmetry of equal probability can again be spontaneously broken and any point on the left flat surface then corresponds to a valid statistical mixture of configurations within the restricted subspace, giving the same free energy f = 0. All these mixtures are represented by the left red plane in the convex sets, where point C corresponds to the equal probability mixture which has maximal entropy s res .
To determine the boundary of the left red plane we proceed the same way as in section 2.2 and simulate Z = z∈Zc e −µ j z j (10) where the sum is now over all proper vertex colorings Z c . The entropy is again given by (9). We have calculated s res for several values of q (see figure 5), where we can reproduce the exact value for q = 3 up to an accuracy of O(10 −10 ) with bond dimension D = 800 of the MPS-representation of the dominant eigenvector of the transfer matrix.
The tensor network used to simulate (10) is described in Appendix A.2.

Critical Exponents and Susceptibilities
If we are given the entire convex set as a function of the extensive observables we can determine critical exponents and susceptibilities purely from the geometrical shape of its surface, i.e. completely independent from the intensive parameters J, T and h. To ease notation in this section we will write Critical exponents for q ≤ 4 can be extracted from the change of the tangent plane orientation around the critical point. For this we need the functional relation between an observable and a model parameter close to the critical point. As an example consider the shifted magnetization z for zero field slightly below the critical temperature T c . There we expect z to behave as with b the critical exponent for the magnetization. We assume the thermodynamic surface to be given e.g. by the interaction energy t as a function of the (independent) variables entropy s and shifted magnetization z, i.e. t = t(s, z). Our intention is to extract b entirely from the geometrical form of the thermodynamic surface, i.e. from the surface given by the function t(s, z). We therefore need a way to express the model parameters J, T and h in terms of the observables t, s and z. From (5) and (6) we saw that they are precisely the elements of the normal vector to the surface function t(s, z). On the other hand, the normal vector to the thermodynamic surface t(s, z) at a given point is and we can immediately identify Without loss of generality we fix J = 1 and consider the case h = 0, i.e. the path of normal vectors with n 3 = − ∂t ∂z = 0. We can then write vs. log z we expect a linear relation near A and we can read off b from the slope ¶.
Estimates for the critical exponents calculated that way from the obtained given numerical data are of the same accuracy as estimates obtained from conventional fits of observables vs. model parameters (i.e. a logarithmic fit of (12)).
Furthermore, susceptibilities defined as the derivatives of the (extensive) observables t, s and z with respect to the (intensive) model parameters J, T and h can be calculated from the curvature of the surface. Loosely speaking, we would like to know how we move on the surface if we change the orientation of the normal vector infinitesimally along one component. In other words, if we change the orientation of n by δT along n 2 , what are the resulting δt, δs and δz. The relation between these changes is of course dictated by the function t(s, z) (or in fact any other representation of the surface, e.g. as s(t, z) or z(t, s)).
With fixed J = 1 we have established the model parameters as functions purely of the observables in (14), i.e. we consider the vector-valued function According to (13) the Jacobian of this function is then proportional to the Hessian of t(s, z) via (18) ¶ As per definition of the ruled surface, z is not unique along this path and it is understood that we take the maximum of z in (15) for each s and t, i.e. the order parameter. This path is nothing but the upper boundary of the ruled surface shown e.g. in figure 1. Alternatively we could have formulated (15) in terms of derivatives of s = s(t, z). Notice however that z = z(s, t) is not a good choice as it is a highly multivalued function on the ruled surface.
The inverse function theorem then gives the elements of this Jacobian by inverting (17) and we can thus obtain the susceptibilities from the second derivatives of t(s, z), i.e. we can obtain δO = J O · δp = J −1 p · δp. With the determinant of (17) given by we get e.g.
These relations are only valid if (17) is invertible and the susceptibilities can diverge if (19) becomes zero. Consider for example the magnetic susceptibility χ T . In figure 1 along the path h = 0 we have z = 0 and ∂t ∂z = 0 (i.e. normal vectors with n 3 = 0), as t becomes maximal when z = 0. This is the case for any s along this path, the mixed derivative ∂ 2 t ∂s∂z is thus also zero everywhere. For T > T c the second derivative ∂ 2 t ∂z 2 is finite, but becomes zero as T → T + c (and is in fact zero at every point on the ruled surface per definition). With ∂ 2 t ∂s∂z = 0 and ∂ 2 t ∂z 2 → 0 the determinant of the Jacobian becomes zero as T → T c and χ T diverges.

Conclusions
We have presented an explicit construction of Gibbs' thermodynamic surface in the form of convex sets for the classical q-state Potts model on a two-dimensional square lattice. We established that points on these surfaces correspond to thermodynamically stable states of the model at parameters given by the orientation of the tangent plane going through that point. Points on the inside on the other hand correspond to nonequilibrium states. These convex sets also constitute a novel and very intuitive way of constructing phase diagrams for many body systems, as all thermodynamically relevant quantities are very naturally included in these sets.
In particular we have demonstrated that symmetry breaking phases appear in this sets in the form of ruled surfaces, where the thermodynamically stable state is not unique. Especially for first order phase transitions the critical point gets stretched out into a flat surface, corresponding to the coexistence of phases at the critical point, characteristic for first order phase transitions. As these sets exist in probability space of the physical degrees of freedom prior to any notion of a Hamiltonian, this implies that the occurrence of symmetry breaking phase transitions is purely a consequence of the geometrical structure of probability space. To further elucidate that point we have shown scatter plots of points obtained from random probability distributions, which all lie inside the convex set per construction and give further information about the internal structure of the constructed convex sets.
We have also identified two regimes, where the ground state at T → 0 is exponentially degenerate and the Potts model becomes equivalent to the vertex coloring problem and the 1NN model respectively. The corresponding flat parts in the convex sets constitute all possible states of the system in these regimes. The symmetry of equal weight superposition of these degenerate states can be spontaneously broken on and distinguished by the observables chosen to constitute the convex set, causing the emergence of these flat parts.
Additionally we have shown how thermodynamic relevant quantities such as critical exponents and susceptibilities can be extracted from the curvature of the thermodynamic surface.
In terms of projections of the set of all possible probability distributions of a physical system it remains to investigate and classify all possible ruled surfaces that can arise on such convex sets projections. Some attempts for the case of fully connected graphs have been made in [38,39]. In the context of models of classical statistical mechanics it would be interesting to obtain equivalent convex set representations in the presence of different types of phase transitions, such as e.g. Berezinskii-Kosterlitz-Thouless phase transitions [40,41]. As topological phase transitions in two-dimensional quantum many body systems appear as symmetry breaking phase transitions in the boundary theories of the entanglement degrees of freedom [42], the question remains what would be the equivalent in the context of classical mechanics. in other words, T represents a Matrix Product Operator (MPO) [43,44] decomposition of the DTM. Other tensor decompositions -e.g. yielding the row-to-row or column-tocolumn transfer matrix upon concatenation -are also possible, the advantage of (A.2) is however that the DTM is hermitian for all J, h and β.
We therefore have Z = Tr(T L D ) and in the limit N → ∞ the dominant eigenvalue of the DTM corresponds to the partition function per diagonal z D = Z 1 L of the system (c.f. e.g. [9]).
In order to evaluate the partition function per site z = Z 1 N = z 1 2L D and the local observables in the thermodynamic limit L → ∞ we obtain the dominant eigenvector of the DTM by means of Matrix Product State (MPS) [44] techniques. More specifically, we use a modification of the algorithm presented in [45] for MPOs in the thermodynamic limit [46] to calculate the partition function per site z and an MPS approximation of the dominant eigenvector of the DTM, which can be used to calculate all local observables, in particular δ(z, z ) and z = z − q+1 2 . As we have access to the partition function per site z, we can now easily evaluate the entropy per site, which is given by with the internal energy per site e = −2J δ(z, z ) − h z .
Appendix A.1. Top Plane and the 1NN As described in section 2.2 in order to determine the boundary of the top plane we simulate the trivial perturbation Hamiltonian in the restricted subspace Z t given by all configurations z j ∈ [q, q − 1] such that every z i = q − 1 is completely surrounded by z j = q, i.e. we wish to evaluate There it is also mentioned that (A.5) is equivalent to the 1NN model in a chemical potential [20,31,32,33] by interpreting z j = q, q − 1 as empty and occupied sites of a lattice gas with nearest neighbor exclusion respectively. We can therefore arrive at a formulation of (A.5) where the entropy per site s and the interaction δ(z, z ) are independent of q and z for different q are related by an offset. By substituting z j = q − s j with s j = 0, 1 we get where Z hs is the partition function of the 1NN model and S is the restricted set of all configurations s j ∈ [0, 1] such that every s i = 1 is completely surrounded by s j = 0. The partition functions per site are then related by z tp = e −µq z hs .
To evaluate Z hs we can achieve a summation over the restricted subspace only by summing over all configurations s j ∈ [0, 1], but giving configurations with neighboring s i = s j = 1 statistical weight zero. This way we obtain a MPO decomposition with bond dimension 2, with MPOs given by where the 2 × 2 matrix f is given by The magnetization then becomes giving for the entropy per site δ(z, z ) is invariant as δ(z i , z j ) = δ(s i , s j ) and the expectation value is evaluated with respect to the same probability distribution. Figure B1.
Scatter plot of observables δ(z, z ) , z and s of q-state spins on a two-dimensional square lattice for q = 3 and q = 5, together with the surface of extreme points of the convex set shown in figures 1 and 2. Here we explicitly plot both the upper and lower half of the convex sets for the sake of completeness. Black dots correspond to observables of single random probability distributions with interaction distance R = 1 (see text). Around the green ruled surface, where there is symmetry breaking, the random scatter points clearly cluster along q separate branches, whereas the other accumulation points are a finite R effect.

Appendix B. Random Scatter Plots
In section 2 we have built on the fact that the surface of the convex sets are given by Gibbs states of (1), which can be efficiently simulated using tensor network techniques. As these convex sets however exist in probability space prior to any definition of a Hamiltonian, the occurrence of symmetry breaking is thus purely a consequence of the geometrical structure of probability space. To further elucidate this argument we show scatter plots of points from random probability distributions P (z), not necessarily being Gibbs distributions. The points generated by expectation values with respect to these distributions must therefore all lie on or within the convex set surfaces shown in figures 1 and 2.
In order to simulate random probability distributions we resort to the class of distributions representable by tensor networks consisting of 4-index tensors T P (z) = n T = T z i ,z j ,z k ,z l T z l ,zm,zn,zo T zp,zq,z j ,zr T zr,zs,zm,zt . . . , (B.1) such that the partition function, obtained by summing over all configurations z, is given by a tensor trace Z = tTr n T = z T z i ,z j ,z k ,z l T z l ,zm,zn,zo T zp,zq,z j ,zr T zr,zs,zm,zt . . . , (B.2) which can again be efficiently evaluated using tensor network techniques. The expectation values δ(z, z ) and z can then be calculated the usual way (c.f. Appendix A). For a general (unnormalized) probability distribution P (z) the entropy per site is given by with z the partition function per site. In the special case of a Gibbs distribution, log(P ) is nothing but the internal energy times the inverse temperature βE, which is a local observable (i.e. a sum of local terms). The entropy is then given by the familiar formula s = β(e − f ), with e the internal energy per site and f the free energy per site. For arbitrary P (z) the quantity 1 N log P is in general however not a local observable. With (B.1) on the other hand we essentially restrict ourselves to probability distributions, for which the entropy is given by the sum of local observables and the entropy per site can be evaluated as s = log(z) − 1 2 log T , (B.4) The factor 1 2 comes from the fact that there are half as many tensors T as there are sites on the lattice. We thus obtain random points within the convex set by sampling T (or rather log T ) from some probability distribution and measuring δ(z, z ) , z and s according to (B.4).
The class of distributions given by (B.1) contains all possible nearest-neighbor interactions as well as 4-site interactions around the face on every other plaquette.
Higher order interactions and distances can be achieved in principle by blocking sites, i.e. transforming to variables z i = ⊗ R r=1 z r . The bond dimension of T is then given by q R and only moderate values of R are computationally feasible. As a demonstration we have however resorted to R = 1 and have drawn log(T ) from a gaussian distribution with varying standard deviation σ ∈ [0.2, 1.5]. The resulting scatter plots for q = 3 and q = 5 are shown in figure B1, together with the surfaces of extreme points already plotted in figure 1 and figure 2.
These surfaces are asymptotically obtained by taking the convex hull of more and more random points generated that way with R → ∞. Figure B1 shows that R = 1 already gives a quite good qualitative approximation of the convex set. Especially around the green ruled surface where we expect spontaneous symmetry breaking it is apparent that the points cluster along q distinct branches. This can be interpreted as a signature of the existence of q disjoint probability spaces in the symmetry broken phase and thus statistical mixtures of configurations from different sectors do not correspond to physically realizable states.