Symmetry Breaking and the Geometry of Reduced Density Matrices

The concept of symmetry breaking and the emergence of corresponding local order parameters constitute the pillars of modern day many body physics. The theory of quantum entanglement is currently leading to a paradigm shift in understanding quantum correlations in many body systems and in this work we show how symmetry breaking can be understood from this wavefunction centered point of view. We demonstrate that the existence of symmetry breaking is a consequence of the geometric structure of the convex set of reduced density matrices of all possible many body wavefunctions. The surfaces of those convex bodies exhibit non-analytic behavior in the form of ruled surfaces, which turn out to be the defining signatures for the emergence of symmetry breaking and of an associated order parameter. We illustrate this by plotting the convex sets arising in the context of three paradigmatic examples of many body systems exhibiting symmetry breaking: the quantum Ising model in transverse magnetic field, exhibiting a second order quantum phase transition; the classical Ising model at finite temperature in two dimensions, which orders below a critical temperature $T_c$; and a system of free bosons at finite temperature in three dimensions, exhibiting the phenomenon of Bose-Einstein condensation together with an associated order parameter $\langle\psi\rangle$. Remarkably, these convex sets look all very much alike. We believe that this wavefunction based way of looking at phase transitions demystifies the emergence of order parameters and provides a unique novel tool for studying exotic quantum phenomena.


Introduction
In a series of ground breaking papers in the late 19th 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 nonanalyticities on the surface, which arise by considering convex hulls of analytic functions [4].
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 this work we demonstrate 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.
We show the generality of this geometric interpretation by explicitly extending Gibbs' construction to quantum and classical lattice spin models, as well as quantum field theories. Our work considerably extends and clarifies the convex set picture in the context of the N-representability problem in quantum chemistry [5][6][7][8], and provides a very nice connection to the discussions on the monogamy of entanglement and mean-field theory in [9].

Convex sets for quantum lattice systems
Ground states of quantum many body Hamiltonians composed of local interactions are very special: in order to minimize the energy expectation value they have extremal local correlations, but those correlations must be compatible with the global symmetry of the many body Hamiltonian such as translation invariance. The competition between those two requirements is responsible for the emergence of the typical long-range properties as exhibited in strongly correlated materials. This is best illustrated by the = S 1 2 Heisenberg quantum antiferromagnet: the energy density would be minimized if all nearest neighbor reduced density matrices (RDMs) were singlets, but due to the monogamy properties of entanglement [10,11], this is not compatible with the translation invariance of the ground state. Hence the symmetry requirements smear the entanglement out into a globally entangled state with algebraically decaying correlations and (quasi) long range order.
This competition is nicely reflected in the convex set of the RDMs of all possible pure and mixed quantum many body states of the entire system 3 .
Let us for example consider a lattice of spin-1/2 quantum systems and arbitrary spatial dimension d, take random quantum states r a , and make a scatter plot of the expectation values r á ñ º å , where X and Z are Pauli matrices, N is the number of sites, K is the coordination number of the lattice and the sum is over nearest neighbors only. This is equivalent to restricting to translation invariant states r a TI and measuring the 1 or 2-site observables Z and Ä X X , for which only the 2-site RDM is needed. As these terms do not commute, a large expectation value á Ä ñ X X will lead to a small expectation value á ñ Z , giving rise to a curved boundary of the generated set. Due to the convexity of the set of 2-site RDMs, the generated body is also convex and corresponds to a two-dimensional projection of the full 15dimensional set of all possible 2-site RDMs. The extreme points of this set correspond to ground states of a family of quantum Ising Hamiltonians of the form Indeed, surfaces of constant energy are represented by lines = -á Ä ñ -á ñ E J X X B Z z in this plot, where the orientation of these lines is given by the parameters J, B z and their distance to the origin is proportional to (minus) the energy. Hence the expectation values of the states with minimal energy must correspond to extreme points for which the lines are tangent to the convex set and thus at maximum distance from the origin, or equivalently, every point on the boundary of the generated set corresponds to the ground state of (1) with parameters given by the orientation of the tangent line through that point.
The situation becomes much more interesting and the presence of symmetry breaking becomes immediately evident when adding an extra axis corresponding to the expectation value of to the scatter plot. The extreme points of the resulting convex set now correspond to ground states of the quantum Ising model including a longitudinal field which explicitly breaks the  2 spin flip symmetry In figure 1 we show the surface of this set in zero, one and infinite spatial dimensions (see appendix B.1 for a scatter plot and appendix D for enlarged versions of these surface plots).
For an infinite system in  d 1 spatial dimensions we witness the emergence of a ruled surface with all lines parallel to the new axis, which turns out to be the defining signature for symmetry breaking. Indeed, all points on such a line are ground states to the same instance of (2)-with parameters given by the orientation of the tangent plane 4 -but with different values of á ñ X . This implies that the ground state is not unique and there is symmetry breaking, as an infinitesimal perturbation of the form of a longitudinal magnetic field å  X i i to the Hamiltonian 3 When restricting oneself to RDMs of pure states, this problem is related to identifying the joint numerical range of a set of operators, which is not necessarily convex anymore [12,13]. 4 Here, all tangent planes with normal vectors = ( ) will touch the set on a line of the ruled surface, instead of a single point.
(1) will break the symmetry, and make sure that the magnetization of the ground state will be polarized in the xdirection with a magnitude given by the extreme points of the convex set lying on the border of the ruled surface. á ñ X is then obviously the order parameter, and the shape of the border of the ruled surface encodes all the information about the ground state expectation values such as the order parameter as a function of the (transverse) magnetic field B z .
Furthermore, from figure 1 we observe three additional remarkable but obvious facts. (i) The convex set of the zero-dimensional case completely contains the one-dimensional set, which in turn completely contains the infinite-dimensional case. This reflects the fact that more and more symmetry constraints (e.g. translation invariance along an increasing number of spatial dimensions) restrict the convex set of possible 2-site RDMs further and further. (ii) The ruled surfaces only arise in the thermodynamic limit, hence demonstrating the need for the well known fact that the order parameter is obtained by first taking the limit of the system size to infinity and only then the longitudinal magnetic field to zero. Any set obtained for a finite one-dimensional chain of N spins would look similar to (c). Put differently, with increasing N the surfaces of these sets will be gradually deformed to asymptotically yield the surface of set (b) as  ¥ N , i.e. only in this limit will the green ruled surface and thus symmetry breaking emerge. In [7,8], the concept of speed was introduced to describe the curvature of convex sets of systems on finite lattices, and observed to diverge when doing finite size scaling. (iii) We can extract critical exponents by investigating the geometry of the convex set around the critical point 5 . More generally, any thermodynamic property of the system such as susceptibilities can be read off from this convex set and the properties of its surface, hence demonstrating the power of such convex set plots (for a detailed analysis in the case of classical lattice systems, see 3.4). Figure 1. Convex sets generated by nearest neighbor correlation á Ä ñ X X , transverse magnetization á ñ Z and longitudinal magnetization á ñ X for all possible translation invariant states on an infinite lattice of = S 1 2 spins in (a)  ¥ d and (b) d=1 spatial dimensions, as well as (c) d=0 spatial dimensions (all possible states of two spins). We plot the surfaces of extreme points of these sets, corresponding to ground states of (2) for various values of J, B z and > B 0 x (due to symmetry we only plot the upper half). Blue lines represent points with constant =  J 1 and B x and varying B z . In (b), the black line corresponds to the exact solution for = B 0 x [14] and points A and A′ mark critical points, where at A a ruled surface (green) emerges, signaling a degeneracy of the ground state, which leads to symmetry breaking and a finite value of the order parameter á ñ X . The corresponding order parameter for A′ is the staggered magnetization á -ñ ( )X 1 i i and we would therefore need to add another axis corresponding to that order parameter to observe the corresponding ruled surface. In all three cases, the red line corresponds to J=0 and thus separates regimes of ferromagnetic ( > J 0) and antiferromagnetic ( < J 0) coupling. As J=0 corresponds to the decoupled case the red line consists of product states. It is therefore independent of the spatial dimension d and thus common to all three sets. For enlarged versions of these plots see appendix D. 5 In the one-dimensional case, the dependence of á ñ X on the magnetic field B z in the symmetry broken phase slightly below the critical point e = -B J 1 z with e 1  can be recovered from the orientation of the tangent plane = -á Ä ñ -á ñ E J X X h Z as e = - ¶á Ä ñ ¶á ñ 1 X X Z to find that indeed in this regime á ñ µ + ¶á Ä ñ ¶á ñ

( )
Notice that the depicted convex sets with their ruled surfaces and non-analyticities exist prior to any reference to any underlying Hamiltonian. It is rather the choice of a collection of plotted observables (i.e. a choice of projection of the full convex set of RDMs) that enables access to thermodynamic properties of a corresponding Hamiltonian defined by this collection. The occurrence of symmetry breaking is therefore encoded in the geometrical structure of a certain projection of the convex set of all possible RDMs, and quantum phase transitions are a consequence of the geometry of such convex sets.
In the case where one would like to learn about symmetry breaking phase transitions of a particular Hamiltonian without knowing the order parameter, one can still use a random observable as additional axis. In general a random observable will have a finite contribution of the true order parameter, thus still generating ruled surfaces, albeit not with maximal extent. In order to find the true order parameter one can then optimize over random observables to obtain the maximal extent of the emerged ruled surfaces (see appendix C).
The general features of the plots in figure 1 are clearly generic for second order phase transitions. Note that in the case of first order phase transitions, one would not need to add an extra axis (corresponding to the order parameter) to witness symmetry breaking, and this would already be present at finite size (see section 3 and [9] for examples).
In figure 1 the surface of set (a) for  ¥ d was numerically obtained using semidefinite programming: as a consequence of the monogamy properties of entanglement [10,11] and the quantum de Finetti theorem [15], this set is equivalent to the convex set generated by all separable states [9,16,17]. For the particular case of two = S 1 2 spins, separability is completely determined by semidefinite constraints [18,19] and the surface of the set can be obtained by minimizing the energy = -á ñ -á ñ -á ñ with respect to all separable density matrices of two spins. Set (b) for d=1 was obtained by doing extensive variational matrix product ground state calculations [20], while set (c) for d=0 was obtained by exact diagonalization of a system of 2 spins.

Top plane
. One can easily see that the ground state is exponentially degenerate with growing system size and for d=1 the degeneracy is given by the Fibonacci sequence + F N 1 where N is the number of lattice sites. The ground space is then given by all possible product states such that no two neighboring spins are in the -ñ | x state (see also section 3.2). Any linear combination of these states is a valid ground state and the entirety of all such possible states is given by the top blue plane.
To determine the edge of this plane for d=1 we consider an infinitesimal perturbation a b = å + å H X Z j j j j 1 away from this point, with a b , 1  and project this perturbation onto the Fibonacci The points on the edge of the top blue plane then correspond to ground states of this projected Hamiltonian in the case of d=1 spatial dimensions. In other words, for all values of a b , we seek linear combinations of states within the degenerate ground state subspace which maximize the magnetization along the direction , 0, . As the normal vector of the tangent plane is given by x this plane has the same orientation = -[ ] n 1, 0, 2 but slightly different boundaries in all three cases.

Convex sets for classical lattice systems
A natural question is whether a similar picture emerges in the case of classical statistical physics. As opposed to the competition between non-commuting terms in a quantum Hamiltonian, classical phase transitions emerge as a consequence of the competition between the internal energy E and the entropy S in the free energy = -F E TS. As the free energy is a linear function of energy and entropy, we expect similar convex sets as for the quantum case when making a scatter plot of the expectation values of energy, entropy and the order parameter with respect to all possible probability distributions (see figure 2 for the classical two-dimensional Ising model). Remarkably, we obtain a very similar picture as for the quantum case. The extreme points of the convex set now correspond to expectation values for Gibbs states which minimize the free energy. Note also that pictures involving intensive quantities such as temperature and/or pressure on the axes would not make sense in this setting, as those quantities are not defined for general probability distributions out of equilibrium.
In the following we work out the case for a paradigmatic example of a classical lattice spin model in full detail: the q-state Potts model [22][23][24][25] is a generalization of the ubiquitous  2 -symmetric Ising model [26,27] to  q -symmetry. It has been shown to be correspond to a  q lattice gauge theory of matter [28,29] and in certain parameter regimes to coloring problems [30,31] and hard-square lattice-gas models with nearest neighbor exclusion (1NN) [32].
The Potts model in a magnetic field is defined by the Hamiltonian is a q-state classical 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  q -symmetry, it undergoes a symmetry breaking phase transition at finite critical inverse temperature b = + ( ) q log 1 c [23,33], above which the  q -symmetry is spontaneously broken. For q=2 the Potts model is equivalent to the classical Ising model [26] and can thus be solved exactly in zero field for all temperatures [21,25]. 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 [34,35]. Other solvable cases include < J 0 at  T 0 and zero field for q=3 on the square lattice [30], and q=4 on the hexagonal lattice as well as q=3 on the Kagome lattice [31].
The symmetry breaking phase transition in zero field is continuous for  q 4 and of first order for > q 4 [23]. The nature of the phase transition will be apparent from the geometrical features of the corresponding convex set phase diagrams which we construct below.
Consider the space of all possible probability distributions ( ) z P of configurations of q-state spins = ¼ z q 1, , i 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 Figure 2. Convex set generated by the average nearest neighbor correlation á ¢ñ zz , the entropy per site s, and the expectation value of the magnetization per site á ñ z for all possible probability distributions of classical 2-state spin configurations on an infinite 2D square lattice. The surface of extreme points corresponds to Gibbs states of the classical Ising model on a square lattice, given by and entropy per site where á¼ñ denotes expectation values with respect to ( ) z P . The convex set  is then given by all possible points , , á ñ z and s are compatible with each other, i.e. they stem from a common valid probability distribution ( ) z P . This is an instance of the classical marginal problem [36][37][38][39]. 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 (5).
To see this, consider (hyper)planes in this three-dimensional parameter space, which are defined as families of points Î  X , related by a plane equation of the form , 9 x y z   where n is the normal vector of the plane and d is the distance of the hyperplane to the origin. Setting = n J 2 x , = n h y and = n T z , this yields exactly the (negative of the) free energy per site of (5) d 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 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 3 and 4 respectively (the special case of the Ising model, corresponding to q = 2, is shown in figure 2). 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.2.

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  q -symmetry can be spontaneously broken, such that á ñ 1 z 0. The maximum possible value á ñ z max can then be taken as the order parameter associated to this phase transition 6 . 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 d á ¢ ñ ( ) z z , and s, but different á ñ z 7 . 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 J T 2 , 0, . 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 1, 1 and á ñ > z 0 max the maximum value of the order parameter. 6 Given á ñ z max the shifted magnetization is then á ñ = -+ á ñ ( ( ) )z k q z 1 2 max with k = 1,K, q the integer enumerating the maximally symmetry broken states, characterized by one-site marginal distributions given by , where < p q 1 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 pq 0, 1 . 7 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 [40,41]. A hint towards this fact is given by the peculiar structure of the random scatter plots for quantum and classical systems as shown in appendix B.
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 = 0, as   h 0 will end in a point with á ñ = á ñ 1z z 0 max for < T T c . This nicely reflects the fact that the order parameter can be obtained by first taking the thermodynamic limit at non-zero 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   figure 3 for the case of 5-state spins where the surface of this set is given by Gibbs-states of (5) 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 of the same nature as described in figure 3. 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 broken phase < T T c (see figure 4). 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 d á ¢ ñ ( ) 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 5), 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 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 q j , the magnetic field term is minimized, whereas the interaction part has a positive energy contribution, resulting in a net energy ofq 2 4 per site. If we now flip one spin at an arbitrary position from q toq 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 - , which is only zero for = N 1 f . The two lowest energy states with the smallest magnetization are thus the two Néel states between q andq 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 q q , 1 j such that every =z q 1 i is completely surrounded by = z q j (see also figure 6). This restricted space is equivalent to the configuration space for the nearest-neighbor exclusion latticegas model (1NN) [32] 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 Figure 5. Surfaces of the two-dimensional convex sets generated by nearest neighbor interaction d á ¢ ñ ( ) 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 d á ¢ ñ ( ) 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.
with equal free energy =f q 2 4 . 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  t given by the top plane and m Î . The entropy per site s is then given by is the partition function per site. The other observables d á ¢ ñ ( ) z z , and á ñ z are computed as usual but with respect to (12). Note that entropy and d á ¢ ñ ( ) 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 (12) is equivalent to the 1NN model in a chemical potential μ [42][43][44], where states q andq 1 correspond to an empty and occupied site respectively. The limits m  ¥ correspond to the the completely polarized and the Néel states respectively (or equivalently the completely empty and maximally filled lattice  [32] up to machine precision. The tensor network we use to simulate (12) is described in appendix A.1.

Left plane
The left (red) plane with d á ¢ ñ = ( ) z z , 0corresponds 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 q 1, j , 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 [45,46] 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 this problem is trivial as only two valid configurations exist (the two Néel states) and the left plane is absent (see figure 2).
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 [30], for which the value is known exactly as = ( ) s 3 2 log 4 3 res [47,48]. 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 3.2 and simulate where the sum is now over all proper vertex colorings  c . The entropy is again given by (13).
We have calculated s res for several values of q (see figure 7), where we can reproduce the exact value for q=3 up to an accuracy of -( )  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 (14) 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 (9) and (10) 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 , . 17 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 0 t z 3 . We can then write where we have extracted the critical temperature from the orientation of the tangent plane at the critical point A z log we expect a linear relation near A and we can read off b from the slope 8 .
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 versus model parameters (i.e. a logarithmic fit of (16)).
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 dT along n 2 , what are the resulting dt , ds and dz. The relation between these changes is of course dictated by the function ( ) According to (17) the Jacobian of this function is then proportional to the Hessian of ( ) t s z , via so we can express it purely in terms of the observables. The infinitesimal change in the normal vector when moving infinitesimally on the surface is then given by d d We are however interested in the converse direction, i.e. the derivatives which are the elements of the Jacobian of the inverse function = The inverse function theorem then gives the elements of this Jacobian by inverting (21) and we can thus obtain the susceptibilities from the second derivatives of ( ) t s z , , i.e. we can obtain d d With the 8 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 (19) 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 3. Alternatively we could have formulated (19) 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. determinant of (21) given by we get e.g.

Convex sets for quantum field theories
As a final example, let us consider a 3D quantum system of free bosons in the continuum at finite temperature where we expect to witness Bose-Einstein condensation. Motivated by the above findings we plot the expectation value of the kinetic energy, the entropy, and a U(1) symmetry breaking order parameter y á ñ, with respect to all possible quantum states to obtain a meaningful convex structure (see figure 8). The extreme points of this set correspond to Gibbs states of the Hamiltonian ò y y y y at fixed density r = 1, where again a symmetry breaking term has been added explicitly. A ruled surface emerging below the critical temperature for v = 0 beautifully signals the onset of Bose-Einstein condensation, where the equilibrium state is not unique and can be parameterized by a finite value of y á ñ. Again, the critical exponent can be extracted from the change of the orientation of the tangent plane around the critical point. , the entropy density s and (the absolute values of) the order parameter y á ñfor all possible states of free bosons in three-dimensions with fixed particle density. We show the surface of this set corresponding to Gibbs states of the ideal Bose gas for particle density r = 1 and various values of the temperature T and Î  v . The blue lines correspond to RG flows through parameter space as indicated by the black arrows. These flows can be performed exactly and are given by (33) and (34). The black line corresponds to v=0 and displays a critical point A at a critical temperature T c . Beyond this point a ruled surface emerges, witnessing symmetry breaking and the occurrence of Bose-Einstein condensation.
The system of an ideal Bose gas in the presence of a U(1)-symmetry breaking term can be solved analytically [49] and the thermodynamic extensive quantities plotted in figure 8 are readily found to be with s the entropy density and the chemical potential μ always chosen such that for Î -¥ ¥ ( ) s , .

Conclusion
In conclusion, we investigated the convex structure of RDMs and marginal probability distributions of many body systems and illustrated how the concept of symmetry breaking emerges very naturally through the appearance of ruled surfaces at the boundaries of these sets. As these sets exist without any prior notion of an underlying Hamiltonian, this shows that the reason for the occurrence of symmetry breaking lies in the geometrical structures of the convex sets of all possible RDMs or marginal probability distributions. This picture seems to capture all the thermodynamically relevant features of many body systems in an extremely concise way. It would therefore be very interesting to classify all possible ruled surfaces that can arise on such convex sets. Our work is very close in spirit to the original groundbreaking papers of Gibbs [2,51], which clarified that phase transitions and the coexistence of different phases can be understood in terms of non-analyticities in the parametrization of the surface of thermodynamic diagrams. It is also related to ideas developed in the context of N-representability [5][6][7][8] for describing quantum phase transitions in fermionic systems. It provides an explicit construction of the famous thermodynamic surface of Maxwell [3] for the case of classical and quantum spin systems, and illustrates very concisely the mathematical physics point of view of symmetry breaking as a breakdown of ergodicity [41]. It also complements the ideas developed in [52,53], where a systematic procedure for finding order parameters was developed by contrasting the RDMs of the low-lying excited states in finite size quantum many body systems. Note however that the starting point of our work is very different: we make no a priori reference to a Hamiltonian, and just consider scatter plots with respect to all possible many body wavefunctions and/or probability distributions. Only the choice of plotted observables relates the obtained convex set to the ground/equilibrium state properties of a corresponding Hamiltonian. Finally, the works [16,17] reported on the convex structure of expectation values of separable density matrices; in retrospect, those are the convex sets obtained in the mean-field regime of an infinite-dimensional lattice, as illustrated in figure 1(a).
One remaining open question is how phase transitions of different types, such as e.g. Berezinskii-Kosterlitz-Thouless phase transitions [54,55] or especially phase transitions in the ground states of two-and higherdimensional quantum systems with topological order fit in this description, as these cannot be characterized in terms of local order parameters. The tensor network description of quantum states might yield one possible resolution, as the topological order induces certain symmetries onto the virtual boundary theory of the tensor network [56][57][58]. Topological phase transitions then correspond to symmetry-breaking phase transitions in the virtual boundary theory [59], i.e. in the structure of the fixed-point subspace of the transfer matrix of the tensor network. These transitions can thus be characterized in terms of a local order parameter at the virtual level of the tensor network [60]. By bringing this virtual operator back to the physical level, it can be associated to the nonlocal string order parameters that characterize the topological phase [61]. When considering systems on a torus, the natural approach would hence be to plot the expectation value of such a Wilson loop around the torus; the different ground states in the topological phase can then be distinguished by different values of this (nonextensive and non-local) order parameter, and hence a ruled surface should emerge at the topological phase transition. Convex sets for SPT phases have been constructed in [62], which appeared as a follow up to the preprint [63] of the current paper.
In order to evaluate the partition function per site = = z Z z D N L 1 1 2 and the local observables in the thermodynamic limit  ¥ L we obtain the dominant eigenvector of the DTM by means of matrix product state (MPS) [65] techniques. More specifically, we use a modification of the algorithm presented in [20] for MPOs in the thermodynamic limit [66] 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 d á ¢ ñ ( ) 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 A.4 with the internal energy per site . For graphical representations of the tensor network and related quantities see figure A1.
In section 3.2 it is also mentioned that (A.5) is equivalent to the 1NN model in a chemical potential [32,[42][43][44] by interpreting = z q q , 1 j as empty and occupied sites of a lattice gas with nearest neighbor exclusion. We can therefore arrive at a formulation of (A.5) where the entropy per site s and the interaction d á ¢ ñ ( ) z z , are independent of q and á ñ z for different q are related by an offset.
with z the partition function per site. In the special case of a Gibbs distribution, á ñ ( ) P log is nothing but the internal energy times the inverse temperature bE, which is a local observable (i.e. a sum of local terms). The entropy is then given by the familiar formula b = -( ) s e f , with e the internal energy per site and f the free energy per site.
For arbitrary ( ) z P the quantity á ñ P log N 1 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 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 These surfaces are asymptotically obtained by taking the convex hull of more and more random points generated that way with  ¥ R . Figure B2 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.

Appendix C. Order parameter optimization
It is argued in section 2 that in the case where the order parameter with maximal symmetry breaking is not known, one can still use a random observable, which will in general have a finite overlap with the true order parameter. Using this observable as an additional axis will thus yield a convex body showing a ruled surface, albeit not with maximal extent. One can then optimize over all possible observables of that type to find the observable which maximizes the extent of the ruled surface and thus corresponds to the order parameter with maximal symmetry breaking.
We demonstrate this procedure in the case of the quantum Ising model (1). We assume the order parameter to be a local one-site observable and thus optimize over all possible linear combinations of Pauli-matrices X, Y and Z, with fixed spectral radius. If the order parameter is not a one-site observable one would have to optimize over multi-site observables, i.e. over linear combinations of products of Pauli-matrices in the case of spin-1/2 systems. Along that line it would be very interesting to know if it is possible to a priori determine (possibly by different means) just from a random scatter plot if there is a symmetry breaking order parameter at all.
As the observable Z however is already present in the Hamiltonian itself (i.e. we already plot á ñ Z ), it is enough to consider only linear combinations of X and Y. We thus define the one-parameter family of observables Figure B2. Scatter plot of observables d á ¢ ñ ( ) 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 3 and 4. 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. Figure C1. Scatter plot of observables á Ä ñ X X , á ñ Z and á Q ñ ( ) O (as defined in (C.1)) generated from the same set of random states as was also used for figure B1, for p p p Q Î { } 0, 0.35 , 0.75 , . We also draw the convex hulls of the respective sets in green. It is clearly visible that the ruled surface is most prominent in the case Q = 0, where = ( ) O X 0 corresponds to the true maximum symmetry breaking order parameter associated to the quantum phase transition of (1).
And construct convex sets from random states, using á Q ñ ( ) O as a third axis, and vary Θ. Figure C1 shows instances of these sets, together with their convex hulls, for certain selected values of Θ, all generated from the same collection of random states that was also used to generate figure B1. It is clearly visible that the ruled surface in the form of a distinct lobe structure is most prominent for Q = = ( ) O X 0 , whereas there is absolutely no signature of a ruled surface for