Euclidean Dynamical Triangulation revisited: is the phase transition really first order?

The transition between the two phases of 4D Euclidean Dynamical Triangulation [1] was long believed to be of second order until in 1996 first order behavior was found for sufficiently large systems [3,4]. However, one may wonder if this finding was affected by the numerical methods used: to control volume fluctuations, in both studies [3,4] an artificial harmonic potential was added to the action; in [4] measurements were taken after a fixed number of accepted instead of attempted moves which introduces an additional error. Finally the simulations suffer from strong critical slowing down which may have been underestimated. In the present work, we address the above weaknesses: we allow the volume to fluctuate freely within a fixed interval; we take measurements after a fixed number of attempted moves; and we overcome critical slowing down by using an optimized parallel tempering algorithm [6]. With these improved methods, on systems of size up to 64k 4-simplices, we confirm that the phase transition is first order.


Introduction
Euclidean Dynamical Triangulation (EDT) in four dimensions, as introduced below in Sec. 1.1, was first studied by J. Ambjorn and J. Jurkiewicz back in 1992 [1]. They found that the model possesses two phases and the transition between them initially seemed to be of 2 nd order, which is necessary for a continuum limit to be defined. In 1996 then, P. Bialas, Z. Burda, A. Krzywicki and B. Petersson reported for the first time the finding of some 1 st order behavior in this phase transition for systems consisting of N 4 = 32k 4-simplices [5]. Shortly afterwards B. V. de Bakker verified this finding and extended the study to larger systems with N 4 = 64k. However, we were not completely convinced by the numerical methods used in the latter work. In particular, there were three things which disturbed us: 1. Measurements were taken after a fixed number of accepted (instead of attempted) moves, which introduces a systematic error. 2. The use of an artificial harmonic potential to control volume fluctuations also introduces a systematic error. 3. Autocorrelation and thermalization times could easily have been underestimated. Therefore we wanted to cross-check these old results with our own, hopefully correct code which satisfies detailed balance, uses a potential well instead of a harmonic potential to control volume fluctuations, and makes use of parallel tempering to cope with critical slowing down.
The paper is organized as follows: In the remainder of this section we give a brief overview of the EDT model and its phase diagram. In section 2 we describe our simulation methods while in section 3 we present our results: after having verified in part 3.1 that the phase transition is 1 st order, we address in part 3.2 the question whether we can observe a coexistence of the two phases and present therefore a local criterion to determine whether a piece of triangulation is in a crumpled or elongated state and try to identify the nature of the metastability in the Markov chain that causes the 1 st order transition. Finally in the appendix, we propose a modification of the path-integral measure, based on a counting of the number of possible moves in each triangulation, which could weaken the 1 st order nature of the phase transition.

The EDT Model
In 4-dimensional Euclidean Dynamical Triangulation (EDT) [1] the formal path integral for Euclidean (local SO (4) instead of SO (3, 1) symmetry) gravity, with the Einstein-Hilbert action S EH : is regularized by approximating the configuration space (space of all diffeomorphism inequivalent 4-metrics) with the space of simplicial piecewise linear (PL) manifolds consisting of equilateral 4-simplices with fixed edge length a (such manifolds are also called abstract triangulations). Under such a discretization, (1.2) turns into the Einstein-Regge action S ER which for equilateral 4simplices and a space-time of topology S 4 takes the simple form with V n = a n √ n+1 n! √ 2 n being the volume of a n-simplex. The partition function (1.1) can now be written as Z (κ 2 , N 4 ) e −κ4 N4 , (1.5) where after the first equality sign the sum runs over all abstract triangulations T of S 4 and C T is a symmetry or degeneracy factor (to avoid over-counting) which is assumed to be ∼ 1 for sufficiently large systems. After the second equality sign, the canonical partition function
For constant N 4 , we can define (see Fig. 1) a line κ pcr 4 (κ 2 , N 4 ) as a function of κ 2 , along which two phases are separated by a pseudo-critical point at κ 2 = κ pcr 2 (N 4 ). For κ 2 < κ pcr 2 (N 4 ) we are in the crumpled phase where a typical configuration is highly collapsed in the sense that the distance between any two 4-simplices is very short, leading to a large (infinite) Hausdorff dimension. For κ 2 > κ pcr 2 (N 4 ) we are in the elongated phase with Hausdorff dimension ∼ 2, where a typical configuration consists of a so-called baby-universe tree: the total volume is subdivided into smaller parts, the baby-universes 2 , which are pairwise connected by only a small minimal neck 3 . This structure is hierarchical in a treelike manner: consider the largest baby-universe as "mother" with outgrowing smaller "babies" which in turn give birth to their own "babies", and so on (see right-hand side of Fig. 2). The true critical point in the thermodynamic limit is obtained as (1.8) Figure 2: Representative configurations in the crumpled (left, κ2 = 1.26) and elongated (right, κ2 = 1.30) phase at system size N4 ≈ 64k: in the crumpled phase, the triangulation consists of one large, highly connected bunch with outgrowths which are at least an order of magnitude smaller. In the elongated phase on the other hand, although a largest component still exists and may be called "mother universe", it is much smaller than in the crumpled phase and some of its outgrowths (the "babies") are of comparable size.

Pachner Moves
In d dimensions, there exist (d + 1) Pachner moves. They form an ergodic set of local updates [15] in the space of abstract triangulations of fixed topology without boundary. A n-move (n ∈ {0, . . . , d}) consists of the following procedure: pick a n-simplex which is contained in (d + 1 − n) d-simplices. The complex consisting of these (d + 1 − n) d-simplices has the same boundary as a corresponding dual complex spanned by (n + 1) d-simplices that share a common (d − n)-simplex. We can therefore just replace the complex around the selected n-simplex with its dual (see Fig.  3 for an illustration in four dimensions). The only additional constraint is the so called manifold constraint, that is: the (d − n)-simplex shared by the (n + 1) newly created d-simplices of the dual complex must not already exist in the triangulation as this could result in topology changes. From now on we will consider only the 4-dimensional case.

Detailed Balance
Calling T k the current triangulation in our Markov chain, we obtain T k+1 as follows: 1. randomly choose a move type n ∈ {0, . . . , 4}, 2. randomly choose one of the N 4 4-simplices of T k and call it D, 3. randomly choose one of the 5 n+1 n-simplices contained in D and call it S, 4. perform a Metropolis test with acceptance probability p n (T k , S): • accept: T k+1 is obtained from T k by applying the n-move at S, The acceptance probability at step 4 is given by [10] p n (T, S) = p n (N 4 (T )) if n-move possible at S ∈ T 0 else , where p n (N 4 ) = min 1, N4 N4+∆N4(n) e κ2∆N2(n)−κ4∆N4(n) is the so-called reduced transition probability, ∆N i (n) labels the change of N i under a n-move, and a n-move is considered as possible at S if S is contained in (5 − n) 4-simplices and the application of the move does not violate the manifold constraint mentioned above in Sec. 2.1. Equation (2.1) can be derived from the detailed balance condition: assume first we have two valid 4-dimensional triangulations T, T , where T can be obtained from T by applying a n-move at a specific n-simplex S of T . The detailed balance equation then reads where ρ (T ) = e κ2 N2(T ) − κ4 N4(T ) and P (T n → T ) is the transition probability. The latter can be written as where p n (N 4 ) is again the reduced transition probability and the factor in front of it is the probability for selecting (with the update scheme mentioned above) the n-simplex S through which T and T can be related by applying a n-move: 1/5 is the probability for choosing the correct move type n, 5−n N4 the one for selecting a 4-simplex D which contains S and 1 ( 5 n+1 ) is the probability for selecting S out of the 5 n+1 n-simplices of D. Note that 5−n ( 5 n+1 ) is the local 4-volume of a n-simplex 4 that allows for a n-move (in units of V 4 ) and p n (N 4 ) can therefore be interpreted as the transition probability 4 The 4-volume containing all points which are closer to the n-simplex under consideration than to any other n-simplex.
for a n-move per unit volume, as the term reduced suggests. For the inverse transition probability on the right hand side of (2.2) we have analogously: which, by noting that N 4 (T ) = N 4 (T ) + ∆N 4 (n), is satisfied by setting p n (N 4 ) = min 1, This gives the upper part of (2.1). For the lower part, when T does not exist, we have ρ (T ) = 0 and in order to satisfy (2.2), we have to set (2.7) We now have a prescription for how to produce a Markov chain containing the configurations required to evaluate (1.5). As neighboring elements in our Markov chain are highly correlated, it is appropriate to take measurements only on every k th element in the chain, where k must be a constant in order to preserve the probability distribution. As already mentioned at the beginning, this was not always respected in previous work, as e.g. in [9] the measurements were separated by a fixed number of accepted moves, which turns the true separation k between measurements into a random variable whose value depends on what kind of configurations is currently sampled by the Markov chain.

Controlling the Volume
As κ pcr 4 κ 2 , N 4 is monotonically growing with N 4 , it is practically impossible to run fully grand canonical MC simulations for the EDT model. But canonical simulations are anyway better suited to investigate finite size scaling. Unfortunately, as already mentioned before, there is no set of ergodic moves known for the space of triangulations of fixed N 4 and it is therefore not possible to run canonical simulations either. The best we can do is to run quasi canonical simulations based on (1.5) but with N 4 constrained to fluctuate around some desired N 4 . In previous work [1,5,9], this was often achieved by adding a harmonic potential, to the action (1.3). This of course introduces a systematic error for all moves which change N 4 . We therefore decided to rather use a infinite potential well of some reasonable width w ≈ 2 σ (N 2 ) / 2.5, ∆N4(n) and σ (N 2 ) is the square root of the N 2 -variance.
As with such a potential well we cannot use the saddle point expansion method from [10] to tune κ 4 to its pseudo-critical value κ pcr 4 κ 2 , N 4 , we instead made use of a method mentioned in [4]: as the N 4 -histogram has to be flat around N 4 if κ 4 = κ pcr where p pcr n (N 4 ) is the reduced transition probability (2.6) with κ 4 = κ pcr 4 κ 2 , N 4 and p geo n (N 4 ) is the average geometric probability for a n-move, i.e. the fraction of the total volume that allows for a change (through one of the Pachner moves), to which a n-move could be applied 5 (see Fig. 4): where f legal n (T ) is the fraction of n-simplices in the triangulation T where a n-move is possible. One can then solve for κ pcr 4 κ 2 , N 4 which leads to Note that the 4-move is always possible and one therefore effectively only has to measure the average fraction of vertices which allow for a 0-move in order to determine κ pcr 4 κ 2 , N 4 . Furthermore, as the average geometric probabilities vary slowly with N 4 , we can, as long as the width of the potential well for N 4 is much smaller than N 4 , just use to set κ 4 to its pseudo-critical value.

Autocorrelation Time
The autocorrelation time τ is a measure for the typical distance between two uncorrelated elements in a Markov chain. It can be thought of as the time it takes for a change to propagate through the typical volume of the system over which the degrees of freedom are correlated. One therefore writes τ ∝ ξ z , where ξ is the correlation length and the dynamical critical exponent z is expected to be z ≈ 2 for a local update scheme; i.e. updated information propagates in a diffusion-like manner.
For a 2 nd order transition where the correlation length diverges when approaching the critical point, ξ is truncated by the linear system size L and we have τ ∝ L z = V z/d H with d H being the Hausdorff dimension of the system. In this case the autocorrelation time obviously diverges as a power of the system size.
For a 1 st order transition ξ remains finite at the transition point for all system sizes. Nevertheless the autocorrelation time can diverge even more dramatically as transitions between the two phases that coexist at the transition point become exponentially suppressed with increasing system size. The autocorrelation function should then consist of two parts: a relatively steep first one corresponding to the decay of autocorrelations within a single phase, as well as a second, much less steep part which reflects the fact that the system remains in one and the same phase for a rather long time. Unfortunately it is almost impossible to verify this as it would need ridiculously long simulations to obtain the required accuracy on the auto-correlation function.
In both cases, for 1 st and 2 nd order phase transitions, parallel tempering can be used to reduce autocorrelations [11,12]. The idea is to run K simulations for different values of the couplings in parallel where the couplings are chosen such that they lie on a line in coupling space which connects a region with slow relaxation with another where relaxation is fast. One now periodically attempts to swap configurations between neighboring sets (called replicas) with an acceptance probability given 5 Alternatively, one could define p geo n (T ) as the fraction of the overall total volume of a triangulation T , to which a nmove could be applied, which would lead to a different normalization: p geo n (T ) = 1 Nn (T ) f legal n (T ), such that p geo n (T ) would coincide with the probability to select a good location for a n-move within the update scheme described above in Sec. 2.2.
are the sets of couplings and configuration variables for the two neighboring replicas and S (κ i , N i ) is the action of a configuration with variables N i at couplings κ i . The advantage of this procedure is twofold: first, we no longer have just one Markov chain per simulation point in coupling space that has to build up the whole corresponding statistical ensemble but now all K chains alternately contribute to all of the statistical ensembles at different couplings. Second, if every Markov chain frequently reaches regions in coupling space where relaxation is fast before passing again through the critical region, then the configurations which the chain contributes to the statistical ensembles in the critical region, are on average much less correlated than corresponding configurations of a Markov chain that remains all the time in the critical region. For more details see [11,12] where it is also explained how this procedure can be optimized. Especially in [12] the application of parallel tempering to 1 st order transitions is discussed.
In our implementation we chose, for a fixed average volume N 4 , 24 or 48 equally spaced (w.r.t. the κ 2 direction) couplings along the pseudo-critical line κ pcr 4 κ 2 , N 4 , such that they join a region in the crumpled phase where relaxation is fast, with one in the elongated phase where relaxation is also fast (compared to the critical region), and thereby pass through the critical region around κ pcr 2 N 4 . After some runtime, the optimization procedure of [11,12] is applied which gives us a new set of couplings for which the configuration exchange between replicas is more frequent. In Fig. 5, two Monte Carlo time histories for the observable N 2 (number of triangles) are shown for comparison, both for a system of size N 4 = 32k at the pseudo-critical point κ 2 = κ pcr 2 (32k) ≈ 1.258. The left one stems from an ordinary simulation while the right one was obtained using parallel tempering.

Data Analysis
Due to the use of a potential well instead of a harmonic potential to control the system volume and due to the tuning of κ 4 to its pseudo-critical value, we have significant volume fluctuations in the data which also affect for example the N 2 distribution. To take this into account, we project the data in the (N 2 , N 4 )-plane along the "correlation direction" before evaluating any observables (see Fig.  6), i.e. instead of N 2 we use (2.14) We checked that this leads to the same results as when evaluating the observables only on data subsets corresponding to single, fixed N 4 values. As long as the fluctuations in N4 are forced to be much smaller than the average system size N 4 itself, we can project the data along the N2-N4-correlation direction (indicated by a red line) and evaluate observables as if we had a true canonical simulation at system size N 4 and fluctuating triangle number N 2 as given by (2.13) and (2.14).
After that, we use multi-histogram reweighting [13] with respect to κ 2 . The parallel tempering optimization procedure mentioned above also leads to a good distribution of simulation points for the reweighting. The errors are determined by the Jack-Knife method with 20 sets. In multi-histogram reweighting, these sets consist of the simultaneous data of all the simulations at different κ 2 values (i.e. to form the Jack-Knife sets we consider as a measurement all the measurements at different κ 2 values which correspond to the same Monte Carlo time), therefore cross correlations should automatically be taken into account.

Order of Phase Transition
As stated in [9], the pseudo-critical N 2 -distribution starts to be visibly double-Gaussian for systems consisting of more than about 32000 4-simplices. For systems containing 32k, 48k and 64k 4simplices, these distributions are shown in Figure 7, where κ 2 was tuned to produce peaks of equal height (left) or equal area (right). Either way it can be seen that the double peak structure becomes more pronounced with increasing system size and that there is so far no sign that the peaks will merge again in the thermodynamic limit. This behavior is characteristic of a 1 st order transition.  Figure 7: Normalized N2-distribution for systems of size N4 = 32k (blue), 48k (red) and 64k (green): the solid lines are double-Gaussian fits to the data. To the left, the values of κ2 were chosen such that the two Gaussian parts of the distribution function have the same height, whereas on the right-hand side the κ2 values are such that the two Gaussians have the same area, i.e. the two states are equally probable. It can be seen that the double peak structure becomes more pronounced with increasing system size and there is no sign that the peaks will merge again in the thermodynamic limit. This is characteristic of a 1 st order transition.
The reason for the double peak structure in the N 2 -distribution is, that at a 1 st order transition point, the two phases can coexist (see Fig. 8) while N 2 takes on different average values in each of these phases.

Scaling of B 4
A more quantitative method to determine the order of a phase transition is to study finite-size scaling of the 4 th order Binder cumulant (Kurtosis) of the According to [14], for large enough N 4 this quantity should scale like as a function of the average linear system size L pcr instead of volume N 4 , as a diverging correlation length will be truncated to this L pcr , which is defined as with n (r; N 4 , κ 2 ) being the average volume profile of a triangulation of size N 4 at a given value of κ 2 , i.e.: where s 1 , s 2 run over all N 4 4-simplices and d ∆ (s 1 , s 2 ) is the geodesic distance between the simplices s 1 and s 2 in the triangulation ∆ and . . . κ2 refers to the average over triangulations at κ 2 . As is typical for a weak 1 st order transition, the 2 nd order fit seems to work fine, but the obtained values  [N2] as a function of average linear system size L pcr (N4) (3.2), assuming a 2 nd order transition, together with a fit of the form (3.1) where N ω 4 = (L pcr ) 1/ν . We also included higher order corrections. It can be seen that the fit seems to work fine, but the obtained values

Coexistence of Phases
As the phase transition is 1 st order we expect for finite systems a coexistence of the elongated and crumpled phases in some neighborhood of the pseudo-critical point. Instead of speaking of crumpled and elongated configurations as in Fig. 8, we should rather speak of configurations in which the crumpled or the elongated part dominates. For example, the left-hand graph in Fig. 8 is dominated by the large bubble 6 in the middle that is in the crumpled state, but attached to that large bubble are baby-universe trees which are in the elongated state. Similarly, the right-hand graph in Fig. 8, whose strong branching indicates that the configuration is mainly in the elongated state, also contains some larger bubbles that presumably are in the crumpled state. Although criteria like "strong branching" and "large bubbles" seem to work well to decide if a piece of triangulation corresponds to the elongated or crumpled phase, there is some ambiguity in what "strong branching" should mean or what bubble sizes should be considered as small. As a first attempt, we could define the elongated phase as consisting of bubbles of size 6 and everything else as belonging to the crumpled phase. By the "size" of a bubble, we mean from now on the number of the bubble's 4-simplices plus the number of its minimal necks 7 (i.e. the volume the bubble would have after replacing all minimal necks by ordinary 4-simplces). In contrast: the "volume" of a bubble still refers to just the number of 4-simplices of that bubble. Figure 11 shows, as a function of κ 2 , the average fractional volume of a triangulation contained in size 6 bubbles. The complementary fractional volume would therefore be the one contained in bubbles of size larger than 6. It can be seen that at small values of κ 2 , deep in the crumpled phase, the volume is dominated by contributions from large bubbles (in fact one very large bubble containing almost all the 4-simplices) whereas at large values of κ 2 , deep in the elongated phase, the dominant contribution to the total volume comes from the size 6 bubbles. The swap in the dominance occurs not exactly at the pseudo-critical point but at a somewhat larger value of κ 2 = κ ds 2 which, however, seems to coincide with the expected infinite volume critical value of the coupling κ cr 2 ≈ 1.33. As the curves in Fig. 11 show no volume dependency in the elongated phase and may get close to  Figure 11: The left-hand figure shows the ensemble averaged fractional volume contained in bubbles consisting of five or less 4-simplices (i.e. bubbles of size 6) as a function of κ2 at system size N4 = 32k. With increasing κ2, it can be seen that at some point κ2 = κ ds 2 (N4), the curve goes above 0.5 and therefore the two classes of bubbles, those consisting of more than five and those consisting of five or less 4-simplices, change their roles as dominant and non-dominant parts of the system. The right-hand figure shows the same quantity but for three different system sizes, N4 = 32k (dark blue), 48k (dark red), 64k (dark green), in a neighbourhood of the corresponding pseudo-critical points: κ pcr 2 ({32k, 48k, 64k}) = {1.258, 1.271, 1.280}. Due to the finite volumes, κ pcr 2 (N4) is smaller than κ ds 2 (N4) which seems to be volume independent and is very close to the expected infinite volume critical value of the coupling κ cr 2 ≈ 1.33.
unity only in the limit κ 2 → ∞, it should be clear that the characterisation of this phase as consisting of just size 6 bubbles is not adequate. The reason is, that if a system would consist of size 6 bubbles only, Pachner 3-moves could only be applied to 3-simplices which are part of minimal necks, such that these moves would necessarily destroy those necks and thereby produce bubbles of size larger than 6. As can be seen in Fig. 12, which shows as a function of κ 2 the quantity ∆necks (n), i.e. the average change of the number of necks under a n-move 8 , indeed, more and more 3-moves change the number of necks as κ 2 increases. But nevertheless, the size 6 bubbles seem to be the dominant building blocks of the elongated phase. Figure 12: The figure shows for a system of size N4 = 32k, as a function of κ2, the average change of the number of necks caused by the next n-move. The 0-and 4-moves always remove or add a "volume 5" bubble and a corresponding neck, while the 2-move never changes the number of necks. For the 1-and 3-moves ∆necks (n) changes as function of κ2 as the fraction of 3-simplices, which allow for a 3-move and are also part of a minimal neck changes. For κ2 > 1. 6, it seems that the average 3-simplex which allows for a move, is rather part of more than one neck than of no neck, which is why ∆necks (3) drops below -1. The behaviour of ∆necks (1) follows from the fact that the 1-moves is the inverse of the 3-move. Thus for κ2 > 1.6, the triangle that will be created by applying a 1-move to one of the 1-simplices that allow for such a move, will be rather part of more than one neck than of no neck. At the pseudo-critical point κ pcr In the pseudo-critical region, in order to change from a rather crumpled to a rather elongated state, the system has to go through the process of producing and growing new baby-universe branches on top of the large "mother universe", until almost the whole volume fits into them while a (distinguishable) "mother universe" disappears. This is illustrated in Fig. 13 where, as a function of κ 2 , it is shown how the total volume of a system with N 4 = 32k is distributed, on average, over bubbles of different sizes.
Regarding this process of creating new baby-universe branches, which necessarily starts with the creation of a new size 6 bubble, it is interesting to note that applying a 4-move to a 4-simplex contained in a bubble that consists of only five 4-simplices results in a triangulation which is slightly more restrictive with respect to further application of 3-and 0-moves, as compared to the case where a 4-move is applied outside of such a "volume 5"-bubble. This is explained in more detail in Fig. 14 for the two-dimensional case, and in figure 15, we show that the effect of this mechanism is indeed observable: the latter figure shows, for different system sizes, the average numbers of possible Pachner n-moves at the pseudo-critical point as a function of N 2 /N 4 . Comparing the graphs in Fig. 15 with Fig. 7 in order to identify which N 2 /N 4 -interval corresponds to which phase, we see that the numbers of possible moves undergo an abrupt change precisely in the region where the valleys of the corresponding graphs in Fig.7 are located. These abrupt changes are just as one would expect from the phenomenon described in Fig. 14 which occurs a soon as a 3-move is applied to a 3-simplex that is part of a minimal neck between a "volume 4"-and a "volume 5"-bubble (which leads to a "volume 11"-bubble i.e. a "size 12"-bubble with just one neck): the number of possible 3-and 0-moves is lower and the number of possible 2-moves higher than in a configuration of the same size but without this particular "volume 11"-bubble.
At the beginning of this section, we mentioned also using the branching factor, at least on an intuitive level, as a criterion to decide if a piece of triangulation is in the elongated or crumpled phase. The branching factor itself turns out not to be a good criterion to distinguish between the two phases as its average value drops again with increasing κ 2 for κ 2 > κ pcr 2 (see Fig. 16). On the other hand, the related average neck density or "branching factor per size" is monotonic and seems to yield a meaningful criterion to distinguish between the two phases (see Fig. 17). In figure 18 we show how the total volume distributes over bubbles with different neck densities but it remains to be understood what the exact value of the critical density, which seems to be around 1/9, should be and where it comes from.  Figure, by applying a 2-move to the general triangle (abc) (which is not part of a "volume 3"-baby-universe), we have replaced the triangle by a baby-universe "A" (horizontally shaded) consisting of the three triangles (abg), (bcg) and (cag). We could continue by applying 1-moves to all the 1-simplices of the neck (abc) which would replace the link (ac) by (gf), (cb) by (gd) and (ba) by (ge). In the right hand figure, we have instead produced another "volume 3"-baby-universe "B" (vertically shaded) with neck (agc) inside the already existing baby-universe "A". In this case, we can only apply a 1-move either to (ac) and (ag) or to (ac) and (gc), as applying a 1-move to both (ag) and (gc) would result in a double link (hb). It is also clear that the insertion of the vertex "h" into one of the triangles of "A" makes it impossible to remove the vertex "g" by a 0-move. If "h" had instead been inserted into e.g. the triangle (acf), not just "h", but also "g" could still be removed. Similarly in the four dimensional case, a piece of triangulation produced by applying a 4-move inside a "volume 5"-baby-universe "A" to produce a new "volume 5"-baby "B" of "A", has the following effects: -the number of possible 0-moves is not increased, as the 4-move that created "B" has also destroyed an already existing location where a 0-move could have been applied: the vertex that was in the center of "A". This does not happen, if a 4-move is applied to a 4-simplex that is not part of a "volume-5"-bubble, -all the 3-simplices that are only part of the neck between "A" and "B" (but not of the neck between "A" and the rest of the triangulation), would allow for a 3-move, but by performing one of these 3-moves, the remaining ones become impossible. If "B" were not the baby of a "size 5,4 or 3"-bubble, 3-moves could in general be applied to all of its neck's 3-simplices, -and finally, in four dimensions, there is also an effect on the 2-move: the application of a 3-move to one of the 3-simplices that are part of only the neck between "A" and "B", removes another location where a 0-move had been possible (as it destroys "B"), but at the same time, generates three new locations where a 2-move could be applied.  Comparison with Fig. 7 shows, that the strong changes (.e.g. for N2/N4 = [2.398, 2.402] in the N4 = 64k case) happen at the location of the valley of the corresponding graph in Fig. 7 Figure 16: The figure shows, as a function of κ2, the average branching factor, i.e. the average number of necks, of bubbles which are neither the largest bubble in the system nor volume 5 bubbles, which are just the terminating leaves of a baby-universe branch. As the branching factor decreases again after the phase transition at κ2 ≈ 1.258, a large branching factor alone is not a good indicator for a bubble to be in the elongated phase.  Figure 17: The figure shows, as a function of κ2, the average density of necks of the largest bubble in the system, i.e. "number of necks of bubble"/"size of bubble" (where the size is again given be the sum of the numbers of necks and 4-simplices). As the largest bubble can be assumed to correspond to the crumpled phase for κ2 < κ pcr 2 (N4) ≈ 1.258, but for κ2 >> κ pcr 2 (N4) is just the slightly largest of many almost equally sized bubbles, which are all part of the elongated phase, this shows, that it is the neck-density of a bubble rather than its total neck number which distinguishes between the two phases. Note also that for size 6 bubbles, the neck density is always ≥ 1/6 ≈ 0.167 which according to this figure is clearly elongated, as it should be. Figure 18: The figure shows for a system of size N4 = 32k, how the total volume is distributed over bubbles with different neck densities, i.e. different ratios "number of necks of the bubble"/"size of the bubble" and how this distribution changes with κ2.

Balls in Boxes Model
It is well known that the balls in boxes model [6,7] describes nicely the qualitative features observed in EDT simulations. In the canonical formulation the model describes the statistical ensemble of a fixed total number N of balls distributed in a varying number M of boxes: where p (q) is the probability for a single box to contain q balls and κ is a coupling introduced to control the number of boxes. As shown in [6], the qualitative behavior of the model depends only on the sub-exponential factors of the single-box-occupation probability p (q), since a redefinition For power like sub-exponential weight factors it was shown in [7] that the free energy has a singularity at κ cr = − log (ζ (β)) for β ∈ (1, ∞), where ζ (x) is the Riemann Zeta-function. The phase is called fluid for κ > κ cr and condensed for κ < κ cr . For β > 2 the phase transition is 1 st order whereas for β ∈ n+1 n , n n−1 the transition is n th order. The order parameter for the transition is the first derivative of the free energy (3.8) with respect to κ, which yields the average number of boxes divided by the number of balls which vanishes in the condensed phase and equals 1 in the fluid phase.
The relation to the 4-dimensional EDT model is normally established by identifying the triangles (or nodes) in the triangulation with boxes and the number of balls in a box with the number of 4simplices that share the corresponding triangle (node). In this way, the coupling κ 2 of the EDT model nicely takes over the role of the κ in the balls in boxes model, the average Regge curvature becomes in the thermodynamic limit the analogue of the order parameter (3.9), and the β in (3.7) is related [7] to the β used in EDT models with a modified measure term [2]. Alternatively, one can identify the bubbles or baby-universes with the boxes and the number of necks of each bubble with the number of balls in the corresponding box [8]. This yields an effective theory for EDT in the form of a branched polymer model, in which the bubbles are the vertices and the necks correspond to links between the nodes (as in the figures 2, 8, but ignoring the different sizes of the nodes). While the latter yields just an effective theory, the problem with the former correspondence is, that due to geometric constraints, the interplay between the number of 4-simplices per triangle (or per node) and the number of triangles (nodes) itself is much more involved than the interplay between the balls and boxes in the balls in boxes model. We would therefore like to propose a different correspondence in which the numbers of "balls" and "boxes" are less constrained. To this end, let us focus on the largest bubble of a triangulation which we will from now on also call base-manifold: this largest bubble is made up of elementary building-blocks consisting of 4simplices and minimal necks 9 . Now consider these elementary building-blocks of the base-manifold as boxes and the number of 4-simplices contained in them as balls. An elementary building-block that is an ordinary 4-simplex corresponds to a box containing a single ball, whereas an elementary building-block consisting of a minimal neck corresponds to a box that contains as many balls as there are 4-simplices in the baby-universe branch behind that neck. The 4 dimensional EDT model therefore corresponds to a balls in boxes model with N 4 balls, where each box contains at least one ball and the number M of boxes can vary from 6 (minimal size for a (combinatorial) simplicial 4sphere) to N 4 (no necks in the triangulation). The canonical EDT partition function could therefore be interpreted as the κ = 0 case of the more general partition function  In terms of (3.10) the average size of the largest bubble and the corresponding susceptibility shown in Fig. 19 could for example be expressed as (3.12) A histogram for the M/N 4 -distribution at the pseudo-critical point is shown in Fig. 20 and should be compared with Fig. 2 of Ref. [6]. As can be seen, Fig. 20 looks much more like Fig. 2 of Ref. [6] than the N 2 /N 4 -distribution shown in Fig. 7, which would be the corresponding quantity according to the old identification: triangles → boxes, 4-simplices → balls. In particular, M/N 4 is a nice order parameter as it tends to zero in the elongated phase, while N 2 /N 4 remains finite. In what follows, we will show that the Z (κ 2 , N 4 , M ) appearing in (3.10) can be written as where the first factor Z 0 (κ 2 , M ) corresponds to the average number of ways a base-manifold consisting of M elementary building blocks (i.e. minimal necks or ordinary 4-simplices) is realized at κ 2 , and the second factor, the sum, is the corresponding probability for N 4 4-simplices to fit into the M elementary building blocks, with p (κ 2 , n) being the probability for a single elementary building block, to have volume n.
To write Z 0 (κ 2 , M ) and p (κ 2 , n) more explicitly, we need the micro-canonical partition function Z 1 (N 2 , N 4 ) that counts the number of possible triangulations with N 2 triangles, N 4 4-simplices and which have a boundary of the form of a minimal neck. For N 4 ≥ 5 each such triangulation can be obtained by removing a 4-simplex from a corresponding triangulation without boundary, that has the same number N 2 of triangles but consists of (N 4 + 1) 4-simplices. Thus Z 1 (N 2 , N 4 ) can be expressed in terms of the ordinary micro-canonical partition function Z (N 2 , N 4 ) as 10 : where 1/5! is the symmetry factor of a 4-simplex and (N 4 + 1) is the number of possibilities to remove one 4-simplex from a triangulation of size (N 4 + 1).
The corresponding canonical partition function is then: where T 1 (N 4 ) is the set of triangulations with a minimal boundary that consist of N 4 4-simplices.
We can now express Z 0 (κ 2 , M ) in terms of (3.15) and (1.6) by noting that the number of ways in which M elementary building blocks can be glued together to form a base-manifold, is the same as the number of ways to form a triangulation, consisting of M 4-simplices, that does not have any neck. We can therefore write: (3.16) where the first term within the brackets on the second line corresponds to the number of triangulations consisting of M 4-simplices and N 2 triangles, while the second term, which is a sum over all possibilities to form a triangulation of size M by glueing two triangulations, each with a minimal boundary, along their boundaries (in [3], this second term was used to measure the entropy exponent by baby-universe counting), subtracts the subset of these triangulations that in addition possess at least one minimal neck 11 , such that the whole bracket yields the number of triangulations with M 4-simplices, N 2 triangles and no necks.
The probability distribution p (κ 2 , n) required for the second factor in (3.13) is given by is the canonical partition function for the interior of triangulations at κ 2 that consist of N 4 4simplices and possess the boundary of a 4-simplex.
The reason for subtracting the whole boundary from the action in (3.17) is, that these terms are already taken into account in Z 0 (κ 2 , M ) and we want to avoid over-counting 12 .
After having written the EDT partition function in the generalized form (3.10), i.e. in terms of a base manifold and its elementary volume elements, which can be excited to form "baby universes", some comments are in order: 1. The terminology "base-manifold" and "elementary building blocks" already suggests that we would like to look at the triangulations, observed in EDT simulations, in a slightly nonstandard way. The main reasons for such a re-interpretation are the following: • it seems that the base-manifold, with all elementary volume elements in the "ground state" (such that they are just ordinary 4-simplices), can be mapped on a corresponding Lorentzian or causal triangulation, • although the boundaries of the elementary volume elements are always minimal, their volume can now change in a discrete manner. This makes EDT to fit a little better into the quantum gravity picture provided by spin-foam models.
2. The altered physical interpretation suggests, that the thermodynamic limit should be taken by sending M , the number of elementary building blocks of the base manifold, to infinity instead of (just) N 4 13 .
3. According to [16], the phase transition is associated with a change of sign in the effective curvature. In the crumpled phase, the base-manifold (or "mother universe") has negative curvature: there are two singular vertices, which could be seen as the centres of two hyperbolic 4-balls (each of them formed by many 4-simplices that are all incident to the same central 11 Again, corrections due to changing symmetry factors (when gluing two triangulations along their minimal boundary) might be necessary in (3.16). 12 Including the boundary terms as in (3.15) into the action for the elementary building-blocks does not work as now at least three boundaries (not just two) meet at each boundary-triangle. The boundary action of an elementary building-block T i is therefore given by S (∂T i ) = κ 2 ∆∈∂T i 1 n(∆) , where ∆ runs through all triangles in ∂T i and n (∆) is the number of boundaries which contain ∆. The p (κ 2 , N 4 ) would then depend (through the n (∆)) on the connectivity of the base manifold, which is highly undesirable. 13 As each elementary building block of the base-manifold contains at least one 4-simplex, M → ∞ implies N 4 → ∞, but the converse is not true. vertex) that are glued along their boundary to form a topological 4-sphere. Without a term in the action that prevents the base-manifold from shrinking, it seems to be favourable for a triangulation to collapse into a baby-universe tree as soon as the singular vertices disappear. Running simulations at quasi fixed M instead of quasi fixed N 4 (but with κ 4 sufficiently large), or at a non-zero value of the new coupling κ, would prevent the triangulation from such a collapse. The base-manifold should then survive the disappearance of the singular vertices and develop a positive effective curvature itself (instead of generating the positive effective curvature by producing many small bubbles), which would give rise to a new phase and a new phase transition that could be of higher than 1 st order. 4. For finite systems, the role played by the new coupling κ in (3.10) is related to that of the anisotropy factor in CDT as κ affects the ratio of the average diameter 14 (∼ average time needed to pass through) and volume of the elementary building-blocks of the base manifold.
In a follow-up paper we will try to verify the above assumptions and study the properties of (3.10) in more detail. An interesting question is of course whether for some values of κ, (3.10) yields a 2 nd or higher order transition in κ 2 (or the fixed ρ = N 4 /M , or fixed M version of (3.10)) and if, when integrating out the volume fluctuations of the elementary building blocks of the base manifold, κ and κ 4 can be combined to yield a kind of effective cosmological constant, such that one recovers the form of the original Euclidean Einstein-Regge action. Alternatively one could interpret the additional weight in (3.10) as a measure term.

Conclusion
Our study confirms the qualitative findings of [5,9]: for κ 2 ≈ κ pcr 2 (N 4 ) we find for N 4 ≥ 32k a clear double peak structure in the N 2 distribution, which becomes more pronounced with increasing system size (and there is no sign that the two peaks will eventually merge again in the thermodynamic limit). This is characteristic of a weak 1 st order transition. A finite size scaling analysis of the 4 th order Binder cumulant of the N 2 distribution confirms this further. As the phase transition is 1 st order, finite systems should allow for coexisting phases in a neighbourhood of the pseudo-critical point. It turned out to be difficult to give a precise criterion to distinguish "locally" between the two phases but a candidate could be that bubbles with a neck density ρ necks > ρ cr necks , can be considered as corresponding to the elongated phase, where ρ cr necks is not yet known exactly but seems to be around 1/9. Bubbles with ρ necks < ρ cr necks would then correspond to the crumpled phase. Finally we proposed a new correspondence between the EDT and "balls in boxes" models which leads to a generalization of the EDT partition function (with an additional parameter κ) and a modified interpretation of triangulations contributing to the EDT partition sum in terms of a largest bubble or "mother universe" and its elementary building blocks, which can undergo volume excitations such that their interior could also be interpreted as a baby-universe branch. The additional coupling κ enriches the phase structure of the model which could now possibly contain a 2 nd order phase transition line. In the Appendix, we propose and motivate a change in the EDT path-integral measure which introduces tunable parameters (r n ). For appropriate choices of the r n , the order of the phase transition of the ordinary EDT model might also change to second order.

Acknowledgement
We thank J. Smit, A. Goerlich and J. Jurkewicz for helpful discussions.

A Appendix: Geometric Probabilities and Path-Integral Measure
It has recently been suggested [2] that the 1 st order transition of the EDT model could perhaps be changed into a 2 nd order transition by a change of the measure in the partition sum (1.5). However, none of these attempts has proved successful. Here we motivate and derive a new proposal for a measure that could have the desired properties.
The measure that is normally used is the trivial one for which we have To motivate a particular form of the measure z (T ), assume we are currently in a triangulation T which possesses f n (T ) locations where a n-move could be applied. This means that T has 4 n=0 f n (T ) "neighboring" triangulations. Now think of each location in T where a move can be applied as something similar to a site in an Ising spin system where a spin-flip can occur. But in our case, the "spin-flip" consists of the application of a Pachner-move, e.g. a Pachner n-move that flips a piece of triangulation, spanned by (5 − n) 4-simplices into one that is spanned by (n + 1) 4simplices but has the same boundary. As the regions where different moves are possible can overlap, the flip of one region will in general destroy some of the other regions where flips were possible and instead create new ones. We therefore have a fluctuating number of degrees of freedom (see Fig. 21) and the system is much more involved than an Ising system.  Nevertheless, one could argue that, as long as the couplings κ 2 and κ 4 are turned off, all possible moves in a triangulation T should be considered as equally likely, just as in the Ising case. More generally: one could assign different probabilities r n to different move types n ∈ {0, . . . , 4}, as long as r n = r 4−n . For example, one could choose r n proportional to the local volume of a n-simplex 15 that allows for a n-move (the local volume is the same for n and (4 − n)-simplices which allow for a move).