Dynamical Gibbs-non-Gibbs transitions in the Curie-Weiss Potts model in the regime beta<3

We consider the Curie-Weiss Potts model in zero external field under independent symmetric spin-flip dynamics. We investigate dynamical Gibbs-non-Gibbs transitions for a range of initial inverse temperatures beta<3, which covers the phase transition point beta=4 log 2 [8]. We show that finitely many types of trajectories of bad empirical measures appear, depending on the parameter beta, with a possibility of re-entrance into the Gibbsian regime, of which we provide a full description.


Research context
The past years have seen progress from various directions in the understanding of Gibbs -non-Gibbs transitions for trajectories of measures under time-evolution, and also more general transforms of measures. The Gibbs property of a measure describing the state of a large system in statistical mechanics is related to the continuity of single-site conditional probabilities, considered as a function of the configuration in the conditioning. If a measure becomes non-Gibbsian, there are internal mechanisms which are responsible for the creation of such discontinuous dependence. This leads to the study of hidden phase transitions, which was started in the particular context of renormalization group pathologies in van Enter, Fernández, and Sokal [33].
Such studies have been made for a variety of systems in different geometries, for different types of local degrees of freedom, and under different transformations. Let us mention here time-evolved discrete lattice spins [19,30], continuous lattice spins [24,34], time-evolved models of point particles in Euclidean space [17], and models on trees [32]. For a discussion of non-Gibbsian behavior of timeevolved lattice measures in regard to the approach to a (possibly non-unique) invariant state under dynamics, see [16], for relevance of non-Gibbsianness to the infinite-volume Gibbs variational principle (and its possible failure) see [22,25]. For recent developments for one-dimensional long-range systems, and the relation between continuity of one-sided (vs. two-sided) conditional probabilities see [2][3][4]31].
In the present paper we are aiming to contribute to the understanding of Gibbs -non-Gibbs transformations for mean-field models, in the sense of the sequential Gibbs property [6, 9-11, 14, 15, 18, 21]. Usually there is a somewhat incomplete picture for lattice models, due to the difficulty to find sharp critical parameters. Mean-field models on the other hand are often "solvable" in terms of variational principles which arise from the large deviation formalism, while the remaining model-dependent task to characterize the minimizers and understand the corresponding various bifurcations can be quite substantial. We choose to work for our problem in the so-called two-layer approach, in which one needs to understand the parameter dependence of the large-deviation functional of a conditional first-layer system. In this functional the conditioning provides an additional parameter given by an empirical measure on the second layer. This is more direct than working in the Lagrangian formalism on trajectory space, which would provide additional insights on the nature of competing histories that explain the current state of the system at a discontinuity point [9,20,28,29].
Compared to the Curie-Weiss Ising model, the Fuzzy Potts model and the Widom-Rowlinson models, we find in the present analysis of the time-evolved Curie-Weiss Potts model significantly more complex transition phenomena, see Theorem 2 and Figure 2. This has to be expected as already the behavior of the fully non-symmetric static model is subtle [23]. It forces us to make use of the computer for exact symbolic computations, in the derivation of the transition curves (BU, ACE and TPE in Figure 2, discussed in Sects. 4.4, 4.5 and 4.6), along with some numerics for our bifurcation analysis. We believe that these tools (see page 42) may also be useful elsewhere. Now, our approach rests on singularity theory [1,5,12,13,27] for the appropriate conditional rate functional of the dynamical model. This provides us with a four-parameter family of potentials, for a two-dimensional statevariable taking values in a simplex. It turns out that the understanding of the parameter dependence of the dynamical model is necessarily based on the good understanding of the bifurcation geometry of the free energy landscape of the static case for general vector-valued fields [23]. In that paper, which generalizes the results of Ellis and Wang [8] and Wang [35], we lay out the basic methodology. Therein we also explain the phenomenology of transitions (umbilics, butterflies, beak-to-beak) from which we need to build here for the dynamical problem.
As a result of the present paper we show that the unfoldings of the static model indeed reappear in the dynamical setup, and acquire new relevance as hidden phase transitions. It is important to note that, in order for this to be true, we have to restrict to mid-range inverse temperatures β < 3. More work has still to be done to treat the full range of inverse temperatures for the dynamical model, where more general transitions seem to appear for very low temperatures. For the scope of the present paper, it is this close connection between the static model [23] in fully non-symmetric external fields, and the symmetrically time-evolved symmetric model in intermediate β range, which is really crucial to unravel the types of trajectories of bad empirical measures of Theorem 2. It would be challenging to exploit whether an analogous non-trivial connection, that we observe for our particular model, holds for more general classes of models. This clearly asks for more research.

Overview and organization of the paper
In the present paper we study the simplest model which is, together with its time-evolution, invariant under the permutation group with three elements: We consider the 3-state Curie-Weiss Potts model in zero external field, under an independent symmetric stochastic spin-flip dynamics. Based on previous examples [21], one may expect loss without recovery of the Gibbs property for all initial temperatures lower than a critical one (which then may or may not coincide with the critical temperature of the initial model), and Gibbsian behavior for all times above the same critical temperature. We show that this is not the case for our model, and the behavior is much more complicated:  The trajectories of the model show a much greater variety, depending on the initial temperature. We find a regime of Gibbs forever (I), a regime of loss with recovery (II) and a regime of loss without recovery (III). Figure 1 shows the non-Gibbs region in the two-dimensional space of initial temperature and time. The boundary of this non-Gibbs region consists of three different curves which correspond to exit scenarios of different types of bad empirical measures. Bad empirical measures are points of discontinuity of the limiting conditional probabilities as defined in Definition 1. Under the time evolution t ↑ ∞ (or equivalently g t ↓ 0 given by (4)) the system moves along vertical lines of fixed β towards the temperature axis. Intersections with a finite number of lines occur along this way, which are responsible for the transitions described in our main theorem, Theorem 2. These additional relevant lines are shown in Figure 2. Theorem 2 rests on the understanding of the structure of stationary points of the time-dependent conditional rate function given in Formula (9) via singularity theory. It turns out that the bifurcations we encounter for general values of the four-dimensional parameter (α, β, t) ∈ ∆ 2 × (0, ∞) × (0, ∞) (see (6)) are of the same types as for the static model depending on a three-dimensional parameter. However, this holds only if we restrict to mid-range inverse temperatures β < 3 and to endconditionings α taking values in the unit simplex (and not in the full hyperplane spanned by the simplex). Nevertheless, in order to understand the relevant singularities, the analysis is best done by first relaxing the probability measure constraint on the parameter α and allow it to take values in the hyperplane. The analysis proceeds with a description of the bifurcation set, where the structure of stationary points of the conditional rate function changes, and the Maxwell set, where multiple global minimizers appear. To pick from these transitions the ones which are relevant to the problem of sequential Gibbsianness and visible on the level of bad empirical measures, we have to take the probability measure constraint for α into account. This step is neither necessary in the static Potts nor in the dynamical symmetric Ising model. The lines Symmetric cusp exit (SCE), Asymmetric cusp exit (ACE), Triple point exit (TPE) and Maxwell triangle exit (MTE) depicted in the full phase diagram in Figure 2 are examples of such exit scenarios. For those lines there is an exit of a certain particular critical value of α from the unit simplex (observation window). The detailed dynamical phase diagram in Figure 2 shows more information about the transitions during time evolution. Preliminary investigations show that the structural similarity with the static case may no longer be valid in the regime β > 3. Therefore we leave the region of very low temperatures for future research.
We describe the model we are considering together with its time-evolution in Sect. 1.3 where we also define what we mean by Gibbsianness (or the sequential Gibbs property). In Sect. 2 we present our main theorem and describe the transitions of the sets of bad empirical measures as a function of the parameters β and t. We will establish the connection between the analysis of the potential function G α,β,t and the Gibbs property of the time-evolved model in Sect. 3. The analysis of the potential function using the methods of singularity theory is then carried out in the Sects. 4 and 5.

The model and sequential Gibbsianness
We consider the mean-field Potts model with three states in vanishing external field under an independent symmetric spin-flip dynamics. The space of configurations in finite-volume n ≥ 2 is defined as Ω n = {1, 2, 3} n and the Hamiltonian of the initial model is So at time t = 0 the distribution of the model is given by We consider a rate-one symmetric spin-flip time-evolution in terms of independent Markov chains on the sites with transition probabilities from state a to b where We are interested in the Gibbsian behavior of the time-evolved measure The unit simplex contains the empirical distributions of spins. By Gibbsian behavior we mean the existence of limiting conditional probabilities in the following sense.
exists for every family η n,k ∈ {1, 2, 3} with n ≥ 2 and 2 ≤ k ≤ n such that We call α bad, if it is not good. The model µ β,t is called sequentially Gibbs if all α in the unit simplex ∆ 2 are good points.

Dynamical Gibbs -non-Gibbs transitions: main result
Our main result on the dynamical Gibbs -non-Gibbs transitions in the high-tointermediate temperature regime for the initial inverse temperature β < 3 is as follows. This temperature regime ranges from high temperature, covering the phase transition temperature (Ellis-Wang inverse temperature β = 4 log 2), up to the elliptic umbilic point β = 3 (where the central stationary point of the time-zero rate function in zero external field changes from minimum to maximum). Essential parts of the structure of the trajectories of dynamical transitions as a function of time t in the regime β < 3 remain unchanged over the three inverse-temperature intervals I, II and III, which were already visualized in Figure 1. The type of transitions can be understood as deformations of the sequences of transitions found in the static Potts model in general vector-valued fields analyzed in [23], where in that case only the one-dimensional parameter β was varied. Observe that however, the dynamical transitions we describe here, do not necessarily occur in a monotonic order with respect to what is seen in the static model under temperature variation. This is for instance (but not only) apparent in the phenomenon of recovery of Gibbsianness. At very low temperatures (β > 3) different bifurcations seem to occur which will be left for future research. While reading the following theorem it is useful to have Figure 2 in mind as the inverse temperatures and transition times are related to the lines depicted in the dynamical phase diagram. (1 -2) in zero external field, for initial inverse temperature β > 0 and at time t > 0 under the symmetric spin-flip dynamics (3 -5). Then the following holds. (ii) For β BE < β < 8 3 the bad empirical measures are given by three symmetric straight lines in a first time interval t NG (β) < t < t BU (β). For a second time interval t BU (β) < t < t TPE (β), the set of bad empirical measures consists of three symmetric Y-shaped sets not touching. For t TPE (β) < t < t ACE (β) the set of bad empirical measures consists of six disconnected arcs. For t > t ACE (β) the system is Gibbsian again.

Theorem 2. Consider the time-evolved Curie-Weiss Potts model given by
(iii) For 8 3 < β < β * and t NG (β) < t < t BU (β) the bad empirical measures consist of three symmetric straight lines. For t BU (β) < t < t TPE (β), the set of bad empirical measures consists of three Y-shaped sets not touching. For t TPE (β) < t < t B2B (β) the set of bad empirical measures consists of six disconnected arcs. For t B2B (β) < t < t MTE (β) the set of bad empirical measures consists of three disconnected arcs. For t > t MTE (β) the system is Gibbsian again. The inverse temperature β * is given by the intersection point of the two lines B2B and TPE in Figure 2.
(iv) For β * < β < 4 log 2 and t NG (β) < t < t BU  The meaning and computation of these lines are discussed in Sects. 4 and 5. While only the three lines SCE, ACE and MTE appear as part of the boundary line of the non-Gibbs region, the other lines are relevant for structural changes of the set of bad empirical measures. There are lines which are explicit in the sense that they are given in terms of zeros of one-dimensional non-linear functions, for example, the entry time t NG (β) (formula (60)) or the butterfly unfolding time t BU (β) (Formula (72)). The least explicit lines are the MTE and TPE lines which involve a Maxwell set computation, the most explicit line is SCE which is given in parametric form s → (β(s), g t (s)) as described in Proposition 8. Figure 3 gives a graphical overview of the possible types of sequences of bad empirical measures with increasing time for the different temperature regimes. There is an even more detailed graphic that illustrates all the transitions involved in the bifurcation set as well as in the Maxwell set. You can find this graphic in the electronic supplemental material (ESM) under the filename detailed_overview.pdf.

II Recovery
III Loss without recovery ii) straight lines enter the simplex, become non-touching Y-shaped sets at the butterfly transition time t BU (β) and move out of the simplex. The midpoints of the Y-shaped sets exit at t TPE (β) and the set leaves the simplex completely at t ACE (β). In (II.iii) the midpoints of the Y-shaped sets leave the unit simplex at t TPE (β) and the two respective arcs connect at the beak-to-beak transition time t B2B (β). The remaining three arcs move towards the corners and leave the unit simplex at t MTE (β). The exit of the midpoints of the Y-shaped sets and the connection of the six arcs occurs in reversed order in the next row (II.iv). In (III) the central triangle shrinks to a point and forms the star-like set that remains in the simplex forever.

Infinite-volume limit of conditional probabilities
The existence of the infinite-volume limit of the conditional probabilities, that is, the question of sequential Gibbsianness, can be transformed into an optimization problem of a certain potential function. As the parameters (β, t) are fixed throughout this section let us write µ n for the measure µ n,β,t .
has a unique global minimizer, then α is a good point, that is, the infinite-volume limit of the conditional probabilities µ n (·|α n ) with α n → α exists independently of the choice of (α n ).
The idea of the proof goes as follows: We can rewrite the conditional probabilities µ n (·|α n ) in terms of an expected value with respect to a disordered mean-field Potts modelμ n (see Lemma 4). Thus, we have to study the weak convergence of L n , where L n is the empirical distribution of the spins σ 2 , . . . , σ n . Note that this is equivalent to the weak convergence of can prove the theorem by an asymptotic analysis of integrals of the form as was done by Ellis and Wang [8]. So it suffices to prove the Lemmata 4 and 5.
A point is good if the respective random field model shows no phase transition, that is, the law of large numbers holds. To be precise, we have the following representation: Lemma 4. The finite-volume conditional probabilities are given by andμ n is a quenched random field Potts model Proof. The proof follows from explicit computations with conditional probabilities.
This representation of the conditional probabilities transforms the problem of understanding bad points to the analysis of disordered mean-field models and their phase transitions. This analysis is done using the Hubbard-Stratonovič transformation which is successfully used for many models [7,8,21].

Lemma 5. Write
for the empirical measure of n − 1 spins with lawμ n [η 2 , . . . , η n ] • L −1 n . Furthermore, let W be a standard normal random vector independent of L n . The distribution of W/ β(n − 1) + L n has a density proportional to e −(n−1)G αn,β,t with respect to Lebesgue measure.
Proof. Denote by σ 2 , . . . , σ n independent {1, 2, 3}-valued random variables each distributed according to p t (dσ i , η i ) with a fixed boundary configuration η 2 , . . . , η n with empirical measure α n . We denote the expectation with respect to this distribution by E. Then in order to calculate the distribution of we calculate for every bounded continuous function f the expectation Now we apply the transformation m = w/ β(n − 1) + L n and obtain In order to complete the proof, we have to calculate the expectation Now we take the logarithm to raise the expression back into the exponent again. So the expected value (16) of the bounded continuous function f is equal to the following up to a normalizing constant: We can now identify G αn,β,t in the exponent using that

Recovery of the Gibbs property
The regime β < 8 3 is split into three parts given by the intervals (0, β NG ], (β NG , β BE ] and (β BE , 8 3 ). In the first part we find that the model is sequentially Gibbs for all times t > 0 whereas in the other two parts the system recovers from a state of non-Gibbsian behavior. The driving mechanism in this "recovery regime" is due to the butterfly singularity which is already found in the static model [see 23, Sect. 2.4.1]. However, in contrast to the static model the bifurcation set might leave the unit simplex so that in order to answer the Gibbs -non-Gibbs question the location of this set (and the contained Maxwell set) with respect to the unit simplex is also important.

Elements from singularity theory
In order to investigate the Gibbs -non-Gibbs transitions we have to study the global minimizers of the potential G α,β,t (Theorem 3). We will use concepts from singularity theory to derive and explain our results.
Singularity theory allows us to understand how the stationary points of the potential change with varying parameters. This can be achieved by looking at the geometry of the so-called catastrophe manifold, which contains the information about the stationary points of the potential for every possible choice of parameter values. More precisely, it consists of the tuples (m, α, such that m is a stationary point of G α,β,t given by (9). The bifurcation set consists of those parameter values (α, β, t) in ∆ 2 × (0, ∞) × (0, ∞) such that there exists a degenerate stationary point m in R 3 , that is, a point at which the Hessian has a zero eigenvalue. The parameter values of the bifurcation set give rise to a partition of the parameter space whose cells contain parameters at which the number and nature of stationary points do not change. Although we are only interested in α that are bad empirical measures, hence probability measures, it is convenient to loosen this constraint and consider α in the hyperplane H = {m ∈ R 3 |m 1 + m 2 + m 3 = 1} into which the unit simplex is embedded. The following proposition is the basis for the analysis of the bifurcation set.
Then we have the following: (a) Let ρ be any permutation of {1, 2, 3}. Then where we interpret the permutation ρ as a 3 × 3-matrix and M as a column vector. For example, if for two distinct elements a, b of {1, 2, 3}.
(c) The catastrophe manifold of the HS-transform G α,β,t is the graph of the map (m, β, t) → α = χ(m, β, t) given by for m ∈ R 3 . In these coordinates, the β-scaled simplex β∆ 2 is an equilateral triangle in the (x, y)-plane centered at the origin. The Hessian matrix G α,β,t (m) in these coordinates is in block diagonal form: The set of degenerate stationary points is given by the solutions (m, β, t) of the following equation: Before we present the proof, let us stress the importance of this proposition. The matrix Γ naturally appears in the derivatives of G α,β,t and has the two important properties: Firstly, the rows of Γ are probability vectors and secondly the map M → Γ(M, t) is compatible with the symmetry of the model. The fact that the catastrophe manifold is given as a graph allows us to write the bifurcation set as the set of (χ(m, β, t), β, t) such that We can therefore take the same point of view as in the static case [cf. 23, Lemma 3]: We study the zeros of the Hessian determinant as a function of m with β and t fixed. This is a twodimensional problem since we only have to consider points in the unit simplex ∆ 2 . Additionally, ∆ 2 is bounded so that we can simply compute the zeros of the Hessian determinant numerically on a discretization of ∆ 2 as accurately as we want to. In this way we can get insight into the global shape of the bifurcation set. It is convenient to look at this set as composed of the bifurcation set slices B(β, t), that is, the subsets for which the parameter (β, t) is fixed. Figure 4 shows an example of the zeros of the Hessian determinant together with the respective image under the map χ(·, β, t) for a fixed pair (β, t). We now continue with the proof of the above proposition.
Proof of Proposition 6. Let us prove the claims in increasing order. Fix arbitrary M ∈ R 3 and positive t. The following equation proves (22).
We proceed with the second point. Note that the matrix Γ(M, t) can be written as the product DE of the diagonal matrix D = (D a,b ) with entries for a, b ∈ {1, 2, 3} and the matrix Since det Γ(M, t) = det(D) · det(E) and the determinant of D is clearly positive, we have to check that det(E) is positive to see that Γ(M, t) is in the general linear group. We find that the determinant of E is given by χ(·, β, t) χ(·, β, t) The branches of the degenerate points on the left and their corresponding images under χ(·, β, t) on the right are marked with the same color. Note that despite the fact that the degenerate stationary points in the left plot lie inside of ∆ 2 in the right plot we see that parts of the bifurcation set slice lie outside of the simplex. This is a major difference to the static case.
which is clearly positive for all positive g t .
To prove the formula for the inverse, let a, b and d be pairwise different elements of {1, 2, 3}. Substituting the right-hand sides of (23-24), we have the following Adding the right-hand sides of the first three equations yields zero and adding those of the last two gives one. This proves the formula for the inverse.
We now prove that the catastrophe manifold is the graph of χ. First, let us check that the range of χ is indeed the hyperplane H. Take an arbitrary point (m, β, t) in H × (0, ∞) × (0, ∞) and let α = χ(m, β, t).
Since (1, 1, 1) T is an eigenvector of Γ(βm, t) for the eigenvalue 1, it is also an eigenvector of Γ −1 (βm, t) for the same eigenvalue. Therefore, we find so α is an element of H. Next, we show that the catastrophe manifold is the graph of χ. The differential of G α,β,t is given by Since Γ(βm, t) is invertible, the equation G α,β,t (m) = 0 can be solved for α and we find α = χ(m, β, t). Assume α is in ∆ 2 , then G α,β,t (m) = 0 implies that m also lies in ∆ 2 since 0 < Γ b,a (βm, t) < 1 for all b, a in {1, 2, 3}.
To show (27) and (28) observe that the second derivative of G α,β,t is given by the matrix where Γ = Γ(βm, t). The partial derivatives of Γ c,a are elements of the tangent space of ∆ 2 for every c in {1, 2, 3}, that is, summing over a yields zero. Therefore: Since the coordinate basis of the (x, y, z)-chart is an orthogonal basis, we find Since β > 0 and α = χ(m, β, t), the condition for degenerate stationary points det G α,β,t (m) = 0 is equivalent to equation (28).

Universality hypothesis connecting the mid-range dynamical model with the static model
In our work we are guided by the following universality hypothesis, which provides a useful organizing principle to understand the transitions which appear. It is suggested by the universality seen in local bifurcation theory, and verified for our model in the full set of mid-range temperatures β < 3, by means of our analytical treatment in the sequel of the paper, aided in some parts by computer algebra and numerics. There exists a map from the two parameters temperature and time of the dynamical model to one effective temperature parameter of the static model of the form which for our model is defined on the whole subset {(β, t) | 0 < β < 3, t > 0} of the positive quadrant (and not only locally) and this map has the following property. At fixed (β, t) the bifurcation set slice B(β, t) ⊂ ∆ 2 , in the space of endconditionings α for the dynamical model, is diffeomorphic to a subset of the corresponding bifurcation set slice B st (β st ) ⊂ ∆ 2 of the static model under a smooth (β, t)-dependent map  (·, β, t), which we call the effective observation window, always contains the uniform distribution. However, it may be much smaller than ∆ 2 for some parameter values. In fact, this will happen as t ↑ ∞, as we will see. The map β st (β, t) from dynamical to static parameters is (only) uniquely defined on the critical lines EW, B2B and BU of the dynamical model (see Figure 2) which get mapped to the corresponding static values β st = 4 log 2, β st = 8 3 , and β st = 18 7 [see 23, Table 1]. The following conjecture underlies this hypothesis, as it expresses the structural similarity of dynamical and static rate functional, by means of a parameterdependent map acting on the state space ∆ 2 , compare with the definition of equivalent potentials in [27, Chapter 6, Section 1]. (c) For every (α, β, t) in D and every m in ∆ 2 the following identity holds:

Conjecture 7. There exists a set U which contains the unit simplex ∆ 2 and is open in the hyperplane H such that
where f β,α denotes the potential (5) where pr 1 denotes the projection (0, ∞)×∆ 2 → (0, ∞). In other words, the effective static inverse temperature β st does not depend on the dynamical α.
A comparison of Figure 9 with [23, Figure 5] gives evidence for the existence of the map ψ 1 as the bifurcation set slice of the static model looks structurally similar to the bifurcation set slice in a neighbourhood of the unit simplex of the dynamical model. The contour plots in the rightmost plots of the two figures support the existence of the map ψ 2 as the contour plot of the dynamical potential G α,β,t looks structurally similar to a subset of the contour plot of the static potential f βst(β,t),αst(α,β,t) . Note, however, that we are not going to construct the maps ψ 1 and ψ 2 in the following sections of the paper and we do not need to do it. Instead, we explicitly compute the critical lines from the dynamical potential following the ideas of singularity theory. This means that the lines can be found independently of the construction of the maps ψ 1 and ψ 2 . The behavior of the model in the vicinity of these lines follows from Thom's classification theorem [see 26, Section 5 of Chapter 3] and our global analysis is supported by the global numerical analysis of the relevant parts of the dynamical bifurcation set. In the following sections we now proceed with the discussion of the critical lines.

The symmetric cusp exit (SCE) line and the non-Gibbs temperature
The non-Gibbs inverse temperature β NG is defined as the supremum of all β such that µ n,β,t is sequentially Gibbsian for all positive t. It turns out to be a maximum. As the type of transitions of the dynamical model for mid-range temperatures can be understood in terms of the static case, let us remark that in the static Potts model the first type of bad magnetic fields that show up with increasing β are due to three symmetric cusp singularities, the "rockets"

(b) The solutions of the system (44-45) can be explicitly parametrized in the form
for s < 0.

(c) The non-Gibbs temperature is given via
Proof. Let us first prove item 1. A symmetric cusp point α is the image of a symmetric degenerate stationary point m under the map χ(·, β, t) at which the tangent vector of the curve of degenerate stationary points (given by vanishing Hessian determinant) is parallel to the direction of degeneracy. The partial derivatives of G α,β,t with respect to x and z vanish at m because of symmetry, so it is sufficient for a stationary point m to have a vanishing partial derivative with respect to the y-coordinate of m. Now, for the gradient we note that where we have abbreviated Γ b,a = Γ b,a (βm, t) and used the fact that α lies on the simplex edge α = (0, 1 2 , 1 2 ). This yields Equation (44). We will now derive equation (45). Note that the mixed partial derivative, which appears in the degeneracy condition (28), vanishes at partially symmetric points: Plugging in α = (0, 1 2 , 1 2 ), the right-hand side of the last equality in (52) vanishes because Γ 3,3 − Γ 3,2 = Γ 2,2 − Γ 2,3 for points m which have the partial symmetry m 2 = m 3 . Therefore the degeneracy condition (28) is in product form. We calculate the remaining partial derivatives: The partial derivative (53) is always positive for β < 8 3 . This means we only have to consider the zeros of (54). This yields equation (45).
We will now explain the parametrization of the set of solutions given in 2. First note that the variable β can be eliminated from Equation (45) using Equation (44) for all y = 0. When we set w = e gt + 1 we find that the resulting equation is a quotient of quadratic polynomials in w: Since w > 2, it suffices to consider the numerator of the left-hand side. The discriminant of this quadratic polynomial is given by D = ((3y − 1)e 3y + 12y) 2 + 8(6y + e 6y ).
It is positive for all real y. Therefore, this polynomial has two real roots. Because w > 2, we choose the larger of the two solutions where we have defined s = 3y and used the definition of F (s) in Equation (48). Furthermore, F (s) > 4 for s = 0 such that Equation (47) yields positive values for g t . Finally, the non-Gibbs inverse temperature is the minimal value of β along the curve given by the parametrization (46-47). Therefore we calculate the derivative of (46) which gives Since 4e s − F (s) is never zero for any s in (−∞, 0), we only have to consider the numerator of the fraction. We calculate the derivative of F Plugging everything together, dβ/ds = 0 is exactly fulfilled for the zero of the function defined in (50).

Lemma 9. Suppose β lies in the interval (β NG , 3). The entry time t NG (β) into the non-Gibbs region is given by
where y is the largest root in (− β 6 , 0) of y →2 β 2 + 24 βy + 72 y 2 − β 2 + 3 βy − 18 y 2 − 9 β e 6 y − 4 β 2 + 3 βy − 18 y 2 e 3 y . (61) Proof. The entry time t NG is given by the first entry of rockets into the unit simplex while increasing the time t and keeping β fixed. This is because, if the pentagrams unfold at all under increase of time, they unfold after the rockets have entered the unit simplex ∆ 2 . This will be clear in the next subsection where we compute the butterfly line. So let us consider the system (44-45) and fix any positive β < 3. Since the relation (4) between g t and t is strictly monotonically decreasing, we have to look for the maximal g t such that (β, g t , y) with negative y is a solution to the system (44-45), which defines the symmetric cusp exit line. Here, y is a magnetization-type variable. We can solve Equation (44) for w = e gt + 1 to obtain Plugging this into the left-hand side of the degeneracy condition (45), we arrive at 2e −6y 3β 2 2β 2 + 24βy + 72y 2 − (β 2 + 3βy − 18y 2 − 9β)e 6y − 4(β 2 + 3βy − 18y 2 )e 3y = 0. (63) This yields the expression of (61). Since the right-hand side of (62) is increasing with y, we have to pick the largest root of (61).

The butterfly unfolding (BU) line and butterfly exit temperature
The unfolding of the pentagrams is a very important mechanism since it changes the set of bad empirical measures from straight lines to Y-shaped, branching curves. This mechanism is already present in the static case, however, in contrast to the static case we have to deal with the fact that in some parameter regions the pentagrams do not fully lie inside of the unit simplex. This leads us to the definition of a butterfly exit inverse temperature β BE for which at some point in time t > 0 there is a cusp point on an edge of the simplex that is about to unfold into a pentagram. By definition, β BE lies between β NG and 8 3 . The value 8 3 is the first inverse temperature for which a beak-to-beak scenario inside of the unit simplex appears as we will see in Section 4.7. (46-47). The butterfly exit β BE is given by and γ s is the implicit function y = γ s (x) defined in a neighbourhood of (x, y) = (0, s 3 ) by the degeneracy condition (28). Note that equation (65) is explicitly computed by a computer program because its expression is very complicated. Nevertheless it is possible to plot the function (see Figure 6).

Proposition 10. Let v(m, β, t) = (ϕ β ) 2 • χ(m, β, t) be the parallel coordinate of χ(m, β, t) and let β(s) and t(s) be given by
Proof. Let us first fix β between β NG and 8 3 and a positive t. Consider a point α on the midpoint of one of the edges of ∆ 2 such that (α, β, t) belongs to the bifurcation set. Furthermore, without loss of generality by symmetry let us assume that α 2 = α 3 . To this point corresponds a degenerate stationary point m that has the same symmetry m 2 = m 3 . We can solve the degeneracy condition (28) in a neighbourhood of m in the form y = γ β,t (x) such that γ β,t (0) is the y-coordinate of m. In α-space in a neighbourhood of α = χ(m, β, t) we can now write the bifurcation set as χ(ϕ −1 β (x, γ β,t (x), 0), β, t). We know that the parallel component v of α fulfills when we follow the curve γ s through the bifurcation set. This is because it has a minimum before the pentagram unfolds and it has a maximum after the pentagram has unfolded. The curve γ of degenerate stationary points is obtained by solving equation (28) in the form y = γ(x) around (0, y * ) where y * is the parallel component of m * . Let us now compute the second derivative of the v-component of the curve: The other mixed partial derivatives of v vanish sinceγ(0) = 0 because of symmetry.
Furthermore, we computeγ(0) via implicit differentiation: Let us write f (x, y) for the left-hand side of (28) viewed as a function in the unit simplex in (x, y)-coordinates. By implicit differentiation we then find: And therefore:γ Using the symbolic calculus tools (see page 42) we can obtain an expression for (65). Using a similar approach it is possible to compute the line in the dynamical phase diagram for which we find butterfly points no matter where these points are with respect to the unit simplex. The key idea that the parallel component Then the butterfly transition time t BU (β) is given by and s * (β) < 0 is the largest zero of Proof. Using the same reasoning as in the proof of Proposition 10, we find that the point m maps under χ(·, β, t) to a point α that is about to unfold into a pentagram if where γ β,t is obtained by solving the degeneracy condition (28) in the form y = γ β,t (x) in a neighbourhood of the point m. This equation is now dependent on m, β and t, that is, we have one equation and three variables (m is onedimensional because m 2 = m 3 ). Additionally, since we know that the direction of degeneracy is the x-direction, we have the equation This equation can be solved for w = e gt + 1 which yields (70). Plugging this into (74), we are left to find the zeros of (73) for some fixed β in the interval (β BE , 8 3 ).

Reentry into Gibbs: the asymmetric cusp exit (ACE) line
In the β-regime (β NG , β BE ), three pentagrams unfold inside of the simplex at an intermediate time and leave the simplex as t increases further. Since we are interested in phase-coexistence of the first layer modelμ n (Lemma 4) and the phase-coexistence lines of the pentagram end in the asymmetric cusp points of the pentagrams, we must compute the exit time t G (β) of these points for β in the above regime. Like in the previous subsection, this is done using a combination of symbolic and numerical computation (see page 42). First, let us state the problem that we need to solve. (a) There is exactly one branch of solution with m 2 = m 3 and it is given by the graph of a map x → y = γ β,t (x).
Then the asymmetric cusps of the pentagrams are on the simplex edges if Proof. The location of the asymmetric cusps of the pentagrams on the curve x → χ(ϕ −1 β (x, γ β,t , 0), β, t) are given by the local maxima of the parallel component v(x) as a function of the curve parameter x (see Figure 8). This yields (78). Equation (77) comes from the constraint that the cusp point lies on the simplex edge because for points on the edge the parallel component equals − 1 6 β in the chart (26). Now, similarly to the case for the butterfly line, the computation ofγ β,t (x) by hand is impractical. Therefore we compute the expression symbolically with the help of the computer. This allows us to numerically determine the course of the line in the dynamical phase diagram. Now, because it is impossible to solve the degeneracy equation (28) in the form y = γ β,t (x) explicitly, we proceed as follows. Note that it is possible to solve (77) for β and plug it into equation (78). We then fix some value of g t , and numerically solve the system consisting of the degeneracy condition (28), where β is substituted from (77), and equation (78), where γ β,t is substituted by y anḋ where f denotes the left-hand side of (28) considered as a function of (x, y). This yields two equations in the two variables x and y.

The triple point exit (TPE) line
To each of the three pentagrams there belongs a special point, the triple point [see 23,Sects. 3.2]. This point is characterized by the coexistence of three global minima, that is, the functional values of all the three minimizers are equal. First, we discuss the existence of these points and then we determine for each fixed positive β the exit time t triple (β). This is the last time for which there are bad empirical measures with partial symmetry that lie inside the unit simplex.
there exists exactly one α in the hyperplane H with α 1 ≤ α 2 ≤ α 3 such that G α,β,t has precisely three global minimizers.
Since the pentagrams in the bifurcation slices leave the simplex (observation window), it is necessary for a discussion of the bad empirical measures that we find the time when the triple points leave the unit simplex. The problem that we have to solve is stated in the following proposition. The exit time t TPE (β) is then given by t TPE (β) = t(β, y (β)) where ϕ β (0, y (β)) and ϕ β (x(β), y(β)) lie in the fundamental cell m 1 ≤ m 2 ≤ m 3 and the triple (y (β), x(β), y(β)) is a solution to the following system of equations.
Note that the expressions of the equations (84 -86) are computed symbolically by the computer (see page 42 for more information). They are not displayed here because of their length. Figure 9 shows a contour plot of the HS transform G α,β,t with α = (0, 1 2 , 1 2 ) and (β, t) on the line TPE.
Proof. The system of equations mainly comes from two ingredients: equal depth of two minimizers and same end-conditioning α for these two minimizers. The triple point is characterized by a coexistence of three global minimizers and since a triple point α must fulfill the symmetry relation α 2 = α 3 , we find that it is sufficient to compare the two minimizers in the fundamental cell m 1 ≤ m 2 ≤ m 3 .  Figure 10: The beak-to-beak mechanism is characterized by the merging of two horns of two different pentagrams. This merging joins two connected components of the complement of the bifurcation set slice when crossing the red line from right to left. As can be seen in the two rightmost plots, this merging happens on the axis of symmetry. The red dots in the dynamical phase diagram on the left mark the time -temperature pairs that correspond to the bifurcation set slices from left to right. The dots in the central plot correspond to the points of the same color in Figure 11.
Because α 2 = α 3 , we always have one symmetric stationary point so that the two minimizers have the coordinates (0, y , 0) and (x, y, 0). Since we now that either minimizer is a stationary point, we can use the vanishing of the first partial derivative of G α,β,t with respect to the y-coordinate to eliminate the time variable t from the equations. This yields the function in equation (83). Using this function we can eliminate the variable t from the equal depth condition and the other two equations that require that the minimizers belong to the same end-conditioning α.

The beak-to-beak (B2B) line
The beak-to-beak point in the static model is characterized as a cusp point that lies in a segment from the center of the simplex to one vertex, that is, for example it has y > 0. The following proposition describes the line of beakto-beak points and a parametric representation in terms of roots of a cubic polynomial. Note that, despite the fact that the line continues to exist for β > 3, the structural behavior of the bifurcation set around the beak-to-beak point might change in the regime β > 3.

Proposition 15.
Fix any positive β and t, let m be a point in H with coordinates (0, y, 0).
Proof. From the analysis of the static model [see 23, Figure 2, rightmost plot of the first row and neighbouring plots for smaller or larger β] we know that the beak-to-beak point (α * , β * , t * ) is such that if we fix α = α * but change the parameters β or t we either find that α = α * is contained in a cell with two minimizers or in a cell with one minimizer. Since α * lies on the axis of symmetry, we know α * = χ(m * , β * , t * ) where m * lies on the axis of symmetry as well, and we find in coordinates ϕ β (α * ) = (0, v(m * , β * , t * ), 0), so it suffices to study as a function of the y-coordinate of m. As before substitute w = e gt + 1. In Figure 11 you see a minimum and a maximum collide and form a saddle point. This is exactly the beak-to-beak behavior. The point (β, t) for which this collision has just happened is given by the vanishing of the first and second derivatives of v(m, β, t) with respect to the y-coordinate of m. Now, the derivatives are given by: Since w > 2, it suffices to consider the numerators of the above expressions. This yields equations (87)  which is only fulfilled for y = 2 9 . However, this leads to the contradiction e gt = 1 e 4 3 −1 < 1 but g t > 0. Therefore, we can assume that we can solve (88) for β. Plugging this into equation (87) we arrive at the following fraction of polynomials in w.
We will now discuss the roots larger than 2 of this cubic polynomial. It is convenient to change variables θ = w − 2, so that we are interested in the positive roots of the following polynomial: Using Descartes' rule of signs, we know that the number of positive roots is equal to the number of sign changes among consecutive, nonzero coefficients of the polynomial or it less than it by an even number. Note that the coefficients in increasing order for s = 0 are given by (9, 12, 3, 0). Therefore we do not find any positive roots for very low positive values of s. The first sign changes appears for the coefficient of order zero which yields equation (92). All of the coefficients except the highest order coefficient eventually become negative. However, with increasing s this happens with increasing order of the coefficient so that we have only one sign change between consecutive coefficients for each s larger than s * . Thus, for all s > s * there exists only one root w * (s) larger than 2.

Reentry into Gibbs: the Maxwell triangle exit (MTE) line
For β in the interval ( 8 3 , 4 log 2) the model displays recovery as well but due to a different mechanism. After the horns of two pentagrams have touched, the Maxwell set which consisted of three connected components now has become one connected component. It consists of three straight lines on the axes of symmetry and a triangle with curved edges. The model recovers from the non-Gibbsianness when this triangle completely leaves the unit simplex which happens on another line in the dynamical phase diagram we call Maxwell triangle exit (MTE). where s > 2 log 2 and w * (s) is the unique zero in (2, ∞) of w → s 1 + (e s − 1)(we s − e s + w) (w + e s )(we s − e s + 2) + log (w + 1) 3 (w + e s ) 2 (we s − e s + 2) . (107) Proof. First, let us derive the system of equations (103-104). Since α has the full symmetry, that is, it is invariant under any permutation of S 3 , it suffices to consider the equal-depth of the central minimum m 0 with one of the three outer ones denoted by m. In the following, we assume m 2 = m 3 . The relative difference between the values is given by where Γ b,a = Γ b,a (βm, t). The partial derivative with respect to the x-coordinate of m vanishes because of symmetry. Plugging in the expressions for Γ 1,1 and Γ 2,1 yields (103). Now, let us come to the parametrization. Equation (105) follows by substituting w = e gt + 1 and s = 3y in equation (103) and solving for β which is possible since s = 0. Plugging this into equation (104) and making the same substitutions we find (107). Note that w * (s) is increasing with s and that the solution of w * (s) = 2 is s = 2 log 2. For lower values of s (107) has no zeros larger than two.

The elliptic umbilics (EU) line
In the static model there is a special point called elliptic umbilic. This catastrophe at the center of the unit simplex is responsible for the fact that the central minimum changes to a maximum. In the dynamical model -due to the additional parameter g t -we have a whole line of these points. This line we call the line of elliptic umbilics (EU).
Using the Taylor expansion of the logarithm and the exponential function, (111) follows by an elementary computation. Note that (113) is actually the HS transform of the static Potts model.
Using symbolic computation with the help of a computer, it is also possible to obtain a Taylor expansion for every pair (β, g t ) on the Elliptic umbilic line. Because of symmetry, the β-dependent coefficients of x 2 y and y 3 differ only by a factor of − 1 3 . This means that for any (β, g t ) on the Elliptic umbilic line the potential G α,β,t with α representing the uniform distribution has the following Taylor expansion up to order three around the simplex center.