Characteristics of the new phase in CDT

The approach of Causal Dynamical Triangulations (CDT), a candidate theory of nonperturbative quantum gravity in 4D, turns out to have a rich phase structure. We investigate the recently discovered bifurcation phase \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{b}$$\end{document}Cb and relate some of its characteristics to the presence of singular vertices of very high order. The transition lines separating this phase from the “time-collapsed” B-phase and the de Sitter phase \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{dS}$$\end{document}CdS are of great interest when searching for physical scaling limits. The work presented here sheds light on the mechanisms behind these transitions. First, we study how the B–\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{b}$$\end{document}Cb transition signal depends on the volume fixing implemented in the simulations, and find results compatible with the previously determined second-order character of the transition. The transition persists in a transfer matrix formulation, where the system’s time extension is taken to be minimal. Second, we relate the new \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{b}$$\end{document}Cb–\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{dS}$$\end{document}CdS transition to the appearance of singular vertices, which leads to a direct physical interpretation in terms of a breaking of the homogeneity and isotropy observed in the de Sitter phase when crossing from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{dS}$$\end{document}CdS to the bifurcation phase \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{b}$$\end{document}Cb.


Introduction
The asymptotic safety program is an attempt to describe quantum gravity as an ordinary quantum field theory. To overcome the well-known nonrenormalizability of the perturbative quantization, the program needs to assume the existence of a nonperturbative fixed point in the ultraviolet (UV). Cona e-mail: ambjorn@nbi.dk b e-mail: jakub.gizbert-studnicki@uj.edu.pl c e-mail: andrzej.goerlich@uj.edu.pl d e-mail: jerzy.jurkiewicz@uj.edu.pl e e-mail: n.klitgaard@science.ru.nl f e-mail: r.loll@science.ru.nl crete continuum calculations using the so-called functional renormalization group equations lend support to this assumption [1][2][3][4][5][6][7], but they necessarily involve truncations. Since the reliability of these truncations is ultimately difficult to quantify, it is important to obtain independent evidence for the existence of a UV fixed point from alternative, nonperturbative methods.
Defining a quantum theory by using a lattice regularization is a well-tested method for obtaining nonperturbative results. The arguably most spectacular results of this kind have been obtained in lattice QCD, where the underlying theory is renormalizable, but many observables cannot be calculated by perturbative methods. Lattice field theories are also well suited to finding nonperturbative UV fixed points, which typically are associated with second-order phase transitions. This means that the first step in a fixed point search consists in localizing phase transition points or lines in the space of bare coupling constants.
In nongravitational lattice field theories the lattice approximates a piece of fixed, flat background space-time and the lattice spacing a acts as a UV cutoff. Given that in General Relativity space-time itself becomes dynamical, it is natural that in a corresponding lattice field theory the lattices themselves should become dynamical entities also. This is precisely what happens in the approach of Dynamical Triangulations (DT) [8][9][10][11][12][13][14][15][16][17][18][19][20][21] and its Lorentzian counterpart, Causal Dynamical Triangulations (CDT) [22][23][24][25][26][27][28][29][30]. Curved space-times, which are summed over in the gravitational path integral, are represented in the lattice regularization by ddimensional "lattices" constructed from elementary building blocks, d-dimensional simplices of lattice link length a, which are glued together in all possible ways compatible with topological and other constraints one may impose. Note that the simplices are not "empty", but are pieces of flat spacetime, such that by assembling them one obtains continuous, piecewise flat manifolds, the said triangulations. The working hypothesis is that in the limit as a → 0 this set of piecewise linear geometries becomes dense in the set of all continuous geometries, assuming that a suitable metric can be defined on the latter. We focus on the CDT rather than the DT approach to nonperturbative quantum gravity, because only in the CDT case one has observed a second-order phase transition which potentially can be used to obtain a UV scaling limit of the lattice theory. 1 Moreover, considering its conceptual simplicity and simple action (see Eq. (3) below), CDT turns out to have a remarkably rich phase diagram, as a function of the bare inverse gravitational coupling κ 0 and the asymmetry parameter . The existence of three distinct phases with corresponding transition lines between them is one of the classic CDT results [26,27,36]. There are two phases A and B in which no meaningful (from the point of view of General Relativity) semiclassical limit seems to exist, a conclusion one arrives at by monitoring the dynamics of the total spatial volume of the universe in time. By contrast, the phase C does display physically interesting behaviour, in that the dynamics generates a quantum universe whose large-scale properties match those of a four-dimensional de Sitter space. While the A-C phase transition was subsequently shown to be first order, the B-C transition turns out to be a second-order transition [37,38], opening the exciting possibility of finding a UV fixed point and an associated continuum theory.
Recently, this picture has been further refined with the discovery of a new transition line cutting diagonally through phase C and dividing it into two regions [39,40]; see Fig. 1. 1 A phase transition observed in DT was originally thought to be second order [19][20][21], but subsequently shown to be first order [31,32]. Recent attempts to enlarge the coupling constant space of DT in search of second-order transition points have so far not been successful [33][34][35].
A first investigation of the order of the new phase transition has not yielded a conclusive answer on whether it is of first or higher order [41]. Since it has now become clear that there are two phases instead of the single phase C, it is a good time to settle on a definite name and notation for them. To ensure continuity with the previous situation and at the same time be descriptive we suggest "de Sitter phase" (C d S ) for the phase above the new phase transition ("above" in the usual κ 0phase diagram), and "bifurcation phase" (C b ) for the phase below the transition. The transition formerly known as the A-C transition then becomes the A-C d S transition, and the former B-C transition becomes the B-C b transition. New is the de Sitter-bifurcation transition The properties of the de Sitter phase C d S coincide with those previously associated with phase C, including the de Sitter-like scaling of the spatial volume. A de Sitter-like scaling is also observed in the bifurcation phase C b , but is modulated there by other dynamical effects, as became apparent when studying the behaviour of the spatial volume in the context of the so-called effective transfer matrix introduced in [42]. In this setting one studies the CDT system with a minimal total number of time steps t tot , typically t tot = 2, compared to the usual t tot = 80. While in the latter simulations inside phase C the entire (de Sitter) universe is visible, in the transfer matrix setting one only has access to a thin "slice" of the universe. Of course, one has to investigate carefully to what extent the two systems describe the same physics (including phase structure and phase transitions), and to isolate finite-size and finite-time effects. Several of the results presented below contribute to this issue.
A major new result found in the transfer matrix approach is the new phase transition C d S -C b , between a phase where the three-volume of adjacent constant-time slices tends to align (C d S ) and a phase where the volume profile is modulated such that the volumes of alternating slices align (C b ). The latter results in a two-peak structure when one plots the volume-volume correlator of neighbouring slices as a function of their (oriented) volume difference [39]. This motivated the term "bifurcation phase", since the corresponding plot in the de Sitter phase C d S has only a single peak. Below, we will uncover a dynamical mechanism behind the bifurcation transition C d S -C b and give it a more direct interpretation in terms of symmetry breaking. At the same time, this will shed some light on the geometric nature of the bifurcation phase, which at this stage is only incompletely understood.
The reason why such an understanding is not straightforward has to do with the nonperturbative character of the dynamics, which is determined by the interplay between the action and the entropy, that is, the number of configurations (triangulated space-times) for given values of the action. An example of this is the behaviour of CDT near the secondorder B-C b transition. The original investigation [37,38] exhibited unusual features, some of them more reminiscent of a first-order transition. Interestingly, as we will see, these first-order aspects disappear when one employs a different prescription for fixing the overall space-time volume. By performing a quantitative analysis of the entropy factor near the transition, we will give a common explanation for both of these phenomena below.
All results presented in this work contribute to the understanding of the dynamical mechanisms determining the behaviour and phase structure of nonperturbative systems of higher-dimensional (in this case four-dimensional) geometry, about which relatively little is known, compared to the wellstudied case of two-dimensional gravity of either signature. To the extent that these properties are driven by "entropic effects", one would expect them to be largely independent of the details of the CDT set-up, and therefore not necessarily confined to this particular approach to nonperturbative quantum gravity.
The remainder of this paper is organized as follows. After a short summary of some vital ingredients of the CDT approach in Sect. 2, we concentrate in Sect. 3 on the second-order B-C b phase transition. We explain a curious dependence of the transition signal on the choice of volume fixing found in previous work by carefully analyzing the entropy factor underlying this behaviour. In the appendix we show that a simple ansatz for this factor can reproduce the characteristic shapes of the transition signals. Section 4 is dedicated to a closer examination of the new bifurcation phase C b . It is performed by simulating an ensemble of CDT configurations with minimal time extension t tot = 2, which is found to display the same phase characteristics and phase transitions as the more customary large-time ensemble. We obtain a quantitative understanding of the properties of the bifurcation phase in terms of a vertex of very high order that appears on one of the two spatial slices of the system. This enables us to give a direct interpretation of the C d S -C b phase transition in terms of symmetry breaking, in this case, the breaking of the homogeneity and isotropy of the average geometry observed in the neighbouring de Sitter phase C d S . A summary and conclusions are presented in Sect. 5.

CDT set-up in a nutshell
We will briefly review the ingredients of the CDT construction and their notation, to the extent that they are needed in the rest of the paper. A comprehensive description of the set-up can be found in [43]. The regularized CDT implementation of the path integral for pure gravity takes the form of a sum over distinct causal triangulations T . After Wick rotation, it is schematically given as the partition function where S E H (T ) is the Einstein-Hilbert action of the piecewise flat manifold T (originally due to Regge) and C T denotes the order of the automorphism group of T , a number equal to 1 in the generic case that the triangulation T does not possess any such symmetries. A triangulation can be thought of as assembled from elementary building blocks, the fourdimensional simplices, which in the standard CDT formulation come in two types, depending on their edge length assignments.
Recall that the interior, flat geometry of a d-dimensional simplex (a "d-simplex") is completely fixed by its edge lengths. CDT configurations have two types of edges, spacelike and time-like. All space-like edges have the same proper length squared a 2 , and all time-like edges the same proper length squared −αa 2 , where α > 0 and a denotes a UV cutoff that will be taken to zero as the regularization is removed. After Wick rotating, which amounts to an analytic continuation of the parameter α to the negative real half-axis in the complex α-plane [43], the triangulations still have two different edge lengths (unless α is set to unity), namely, where α > 7/12 to satisfy triangle inequalities. In addition to the Minkowskian geometry of its simplicial building blocks, the causal character of CDT quantum gravity is reflected in the gluing rules for the foursimplices, which are such that the causal (=light cone) structure of each triangulation T is well defined. In standard CDT this is achieved through the presence of a stacked structure associated with the presence of a discrete time parameter t. 2 A causal triangulation consists of a sequence of threedimensional spatial triangulations, each labelled by an integer t, with four-dimensional space-time simplices interpolating between adjacent slices of constant times t and t + 1. In the present work, the spatial slices will have the topology of the three-sphere.
The two four-simplex types mentioned above are precisely those that are compatible with this stacked or layered structure. They are the (4, 1)-simplex (together with its time-reflection, the (1, 4)-simplex) and the (3, 2)-simplex (together with the time-reflected (2, 3)-simplex). A (4, 1)simplex shares a purely space-like three-simplex (spanned by four vertices) with the three-dimensional triangulation at time t and a single vertex with the spatial triangulation at time 2 There is an alternative version of CDT, using so-called locally causal dynamical triangulations (LCDT) [44], where the causal structure is only implemented locally, without referring to a preferred global lattice time slicing. This can be achieved by introducing new types of building blocks (with edge lengths still given by Eq. (2)). In three space-time dimensions, this approach has produced results compatible with those of CDT [45,46], at the expense of considerable additional computational complexity. Fig. 2 The two types of four-simplex appearing in CDT, the (4, 1)-simplex (left) and the (3, 2)-simplex (right), interpolating between neighbouring spatial slices of constant integer time t. Space-like edges are drawn in blue, time-like ones in red t + 1, whereas a (3, 2)-simplex shares a two-dimensional space-like triangle (spanned by three vertices) with the slice at time t and a space-like edge (spanned by two vertices) with the slice at time t + 1. It follows that a (4, 1)-simplex has 6 space-like and 4 time-like links, and a (3, 2)-simplex has 4 space-like and 6 time-like links (see Fig. 2). Analogous statements hold for the (1,4)-and (2, 3)-simplices when interchanging t and t + 1.
Since there are only two geometrically distinct building blocks, the Einstein-Hilbert-Regge action (including a cosmological constant term) assumes a simple form in terms of the global "counting variables" N i (T ), i = 0, 1, . . . , 4, which for a given triangulation T count the number of idimensional simplices contained in T . Below, we will use the numbers N 0 of vertices and N 4 of four-simplices. It will be essential to keep track of the separate numbers N (4,1)  is the number of (4, 1)-and (1, 4)-simplices together. Since they occur frequently in our formulae, we will use N 41 := N (4,1) 4 and N 32 := N (3,2) 4 as a shorthand notation. Of course, we have N 41 (T )+ N 32 (T ) = N 4 (T ) for any T . In terms of these, we can finally write the gravitational action as [43] where κ 0 is the bare inverse Newton constant, κ 4 (up to a κ 0 -dependent shift) the bare cosmological constant, and is an asymmetry parameter that depends on the finite, relative scaling α between time-and space-like links introduced in (2). Details of this algebraic dependence will not concern us here, other than the fact that vanishes for equilateral simplices, that is, (α = 1) = 0. In the nonperturbative regime investigated by CDT, plays the role of a coupling constant. To emphasize various aspects of the action (3), whose motivation will become clear in subsequent sections, we can rewrite it in a number of equivalent ways, Equation (4) is a straightforward reshuffling of terms, Eq. (5) is a rewriting of (3) using the difference x := N 41 − N 32 , while (6) results after performing a linear redefinition of the coupling constants according to κ 0 := κ 0 + 6 , κ 41 := κ 4 + 2 and κ 32 := κ 4 + .
In the actual CDT computer simulations the lattice volume is kept (approximately) constant, by adding a volume-fixing term S fix to the bare action (3). This means there are de facto only two tunable bare couplings, κ 0 and , as illustrated by the phase diagram of Fig. 1. Two different quadratic volume fixings have been used in the literature, either fixing the total number of four-simplices toN 4 by setting or fixing the number of (4, 1)-simplices to some target valuē N 41 by setting where ε in both cases denotes an appropriately chosen small, positive parameter. Inside the "old" phase C and well away from the phase transitions B-C b and A-C d S one does not expect results to depend on the volume fixing used, since at a given (κ 0 , ) the two four-simplex types occur approximately in a fixed ratio [43]. However, as already mentioned above, some measurements at the B-C b transition appear to depend on the volume fixing, a phenomenon that will be explained in Sect. 3.

A second look at the B-C b transition
We begin by examining the transition between phase B and the bifurcation phase C b . It has been known for some time to be a second-order transition, and thus potentially interesting for continuum physics. The original investigation of what was then called the B-C transition was performed at fixed N 4 , implemented by a volume fixing of the form (7), for volumes of up to N 4 = 160k [37,38]. The order parameter chosen to study the transition was conj( ) := N 41 −6N 0 , which is the expression conjugate to at fixed N 4 , as can be read off from (4). The analysis required some care, because the probability distribution of conj( ) measured at the transition exhibited a double-peak structure. This is unusual, because a double peak is typically associated with a first-order transition, where it is brought about by a jumping of the order parameter between two metastable states on either side of the transition. However, in the case at hand a careful analysis of finite-size effects in terms of observables like the Binder cumulant, particularly suited to distinguishing between firstand higher-order transitions, all pointed towards a secondorder transition.
We have found it convenient to work with another order parameter, the quantity x = N 41 − N 32 introduced earlier.
Looking at the action (5), one observes that x would be conjugate to for fixed N 4 if we also held N 0 fixed (which we do not). Using x instead of conj( ) as an order parameter corresponds to approaching the transition line along a slightly different phase space trajectory, and it leads to an equivalent result for its probability distributionP(x). 3 The results forP(x), measured at fixed N 4 = 40k and for time extension t tot = 80, are shown in Fig. 3 and display the same kind of double peak as in the original work [37,38]. Note that the relative height of the two peaks in the distributionP(x) depends on the coupling . We define the critical value c as the value where the peaks have the same height. 4 Following a space-time configuration and measuring its x-value as a function of Monte Carlo time, one finds that x is located close to one of the peaks for some time and then makes a very rapid change to the other peak where it again stays for some time. (Examples of Monte Carlo time histories of order parameters, albeit in a slightly different context, are depicted in Fig. 9 below.) This is precisely the behaviour expected at a first-order transition, for sufficiently small volumes. However, for a genuine first-order transition such a cross-over between different phases will be suppressed as the system size goes to infinity. The absence of such a behaviour for increasing volume led to the more detailed investigation of [37,38], with the outcome that the B-C b transition in CDT appears to be of higher order.
Somewhat surprisingly, when repeating the same measurements with N 41 rather than N 4 kept fixed, we found no

Fig. 3
Probability distributionP(x) of the order parameter x, measured at three different couplings close to the critical point c ≈ 0.0220, for total volume N 4 = 40k and κ 0 = 2.2 trace of a double-peak structure for any of the order parameters considered. The distribution of x (which for constant N 41 coincides with the distribution of N 32 ) is shown in Fig. 4. As explained in more detail in Sect. 3.2 below, we have determined the (pseudo-)critical value c from a peak in the susceptibility χ(x) = x 2 − x 2 under variation of , where the distributionP(x) has maximal width. Thus it appears that for fixed N 41 the situation is consistent with that of a typical second-order transition.
In the following, we will demonstrate that the observed dependence of the distribution of x on the volume fixing has its origin in the function that counts the number of configurations (including their symmetry factors 1/C T ) for given values of the counting variables N 0 , N 41 and N 32 , the entropy (factor) where T (N 0 , N 41 , N 32 ) denotes the set of triangulations with fixed N 0 , N 41 and N 32 . Using the action in the form (6), the partition function can now be written as We will apply Monte Carlo techniques to extract the entropy N (N 0 , N 41 , N 32 ). In order to measure this function over a whole range of values in the (N 41 , N 32 )-plane, as we would like to do, an efficient method is to modify the action in a controlled way such that one probes smaller regions. By adding quadratic terms, x Generic Cross section to the action (6), one ensures that the Monte Carlo simulations probe a well-defined, not too large region in the vicinity of a prescribed point (N 41 ,N 32 ). More specifically, a given set of numbers N 0 , N 41 and N 32 will occur with probability We have covered the region of interest by eight patches corresponding to different valuesN 41 ,N 32 , such that they overlap mutually. This allows us to adjust the relative probability distributions measured in the different patches to a common probability distribution, which is determined up to a common normalization factor. We could in principle have chosen different values for the three couplings κ 0 , κ 41 and κ 32 in the various patches, but we keep them constant across all patches and equal to the reference valuesκ 0 ,κ 41 andκ 32 .
To simplify the comparison between fixing N 4 and N 41 , we integrate out the number N 0 of vertices weighted by eκ 0 N 0 to obtain the "reference" probability distribution where the normalization factor C ensures that the probabilities add up to one. The distribution (13) can be extracted from the measured quantities PN It is understood that during the matching process for the overlap regions the various PN 41 ,N 32 have been normalized rela-tive to each other such that after multiplication with exp(S fix ) and summing over N 0 only a single common normalization factorC is needed, as already mentioned above. The righthand side of Eq. (14) therefore describes a single, joint probability distribution, which by construction no longer depends onN 41 andN 32 .
Rather than working directly with P(N 41 , N 32 ), we have found it convenient to work with its logarithm which can be interpreted as (minus) the free energy of the system. The density plot of the measured free energy (15)  In connection with Fig. 3 we already reported on the double peak in direct Monte Carlo simulations of the probability distributionP(x) observed for fixed N 4 . Remarkably, the same double peak can be reproduced by taking a cross section where again x = N 41 − N 32 . This is illustrated in Fig. 6 (blue dotted curve). The fact that we can reconstruct the double peak in this way shows that the saddle-shaped geometry of the free energy F(N 41 , N 32 ) is responsible for this structure. In other words, in the volume range considered, the occurrence of such a double peak is caused by "entropy", in the sense of the distribution of configurations contributing to the path integral, and is not an indication of the presence of a firstorder transition.
3.2 Single-peak structure and transition for fixed N 41 By contrast, for fixed N 41 , implemented by adding the volume-fixing term (8) to the action, the distributionP(x) is well approximated by a concave function with a single "Gaussian-like" bump as illustrated by Fig. 4. The violet curve shows the results of standard Monte Carlo simulations forP(x), while the blue line represents The corresponding cross section through the (N 41 , N 32 )plane is given by the vertical grey line in Fig. 5. The two methods for determining this distribution are in perfect agreement. Note also that the maximum ofP(x) of Fig. 4 and the minimum ofP(x) of Fig. 3 occur approximately at the same point, namely, N 41 = 33k, N 32 = 8k. The free energy F(N 41 ,N 41 −x), together with a quadratic fit, is shown in Fig. 7. As mentioned earlier, by looking at where the standard deviation σ (x) of the distributionP(x) for N 41 =N 41 peaks as a function of the coupling , we can extract the critical value of . To obtain the standard deviation ofP(x) one can proceed in two different ways. One option is to simply perform Monte Carlo simulations at fixed N 41 for a number of selected values of (yellow dots in Fig. 8). The other procedure (whose results are represented by the blue dots) is more indirect and involves a reconstruction from measurement data taken at fixed .
More specifically, we have taken as a starting point the distributionP(x) displayed in Fig. 4, which was measured for fixedκ 0 = 2.3320,κ 41 = 0.9856 andκ 32 = 0.9636, and therefore corresponds to the single, fixed value¯ := κ 41 −κ 32 = 0.0220. Since N 41 is kept fixed, the relevant coupling constants areκ 0 andκ 32 . Due to the simple form of the action (6), there is an easy relation which allows us to Since we are keeping κ 41 fixed, a change in κ 32 is equivalent to a change in , in the sense that =¯ +κ 32 − κ 32 , which is exactly what we are interested in when determining the standard deviation σ (x) ofP(x). The only limitation to be taken into account when constructing σ (x) from numerical data in this way is that κ 32 should not differ too much fromκ 32 . One typically has accurate measurements ofPκ 32 (x) only for some limited range in x, which means that for |κ 32 − κ 32 | too large the centre ofP κ 32 (x) will be shifted to an x-interval wherePκ 32 (x) is poorly determined, and thus will lead to a large uncertainty in the derived distributionP κ 32 (x). As can be seen in Fig. 8, in the case at hand the two very different ways of determining the standard deviation agree remarkably well, especially with regard to the location of their peaks. This has allowed us to extract the critical value of with good accuracy as c ≈ 0.026. The fact that this differs slightly from the mea-surement at fixed N 4 is not particularly surprising, since at finite volume the two volume fixings lead to systems with different behaviour.
In the appendix, we make a simple ansatz for the free energy F(N 41 , N 32 ) in terms of several free functions at most quadratic in N 41 and N 32 , which we determine uniquely from fitting them to our data. This ansatz reproduces the features described in this section: a cross section N 4 = const results in a double-peak structure and a cross section N 41 = const in a single-peak structure for the probability distribution of x = N 41 − N 32 . At the same time, the ansatz is too simple to reproduce the observed higher-order critical behaviour at the transition. This demonstrates explicitly that the unusual double-peak structure near the B-C b transition is not necessarily related to any critical behaviour and the question whether the observed transition is of first or second order.

The bifurcation phase
Having exhibited one aspect of the nonperturbative dynamics of CDT near the B-C b transition, we now turn to a closer analysis of the bifurcation phase C b , including the associated, new C d S -C b transition. The results we will discuss are obtained in the framework of the so-called effective transfer matrix [42], which was instrumental in the discovery of the bifurcation phase in the first place [39]. This formulation involves the reduced transfer matrix M, whose matrix elements describe the transition amplitudes between a spatial configuration of three-volume m at time t and a neighbouring spatial configuration of three-volume n at time t + 1. They are obtained by measuring the probabilities for a system with a total time extension t tot = 2 [39] and extracting the matrix elements according to The term reduced or effective transfer matrix refers to the fact that of all the geometric degrees of freedom that characterize the three-dimensional spatial slices of constant integer time, one only keeps track of the total three-volume N 3 (t) of the slices at constant t. It is a nontrivial finding that one can reconstruct the well-known effective, "minisuperspace" action and the global dynamics of the three-volume [26][27][28][29][30] from measurements of the reduced transfer matrix alone [39,42]. It was a closer examination of the "unphysical" phases A and, more specifically, B in terms of the effective transfer matrix  [39]. We will study this new phase by concentrating on the correlations between neighbouring spatial slices. To facilitate the investigation and allow for large spatial slices we will consider the situation t tot = 2 with just two spatial slices and periodic boundary conditions. 5 Furthermore, we will keep N 4 fixed by including a term (7) in the action, and set κ 0 = 2.2 throughout.

Equivalence with large-time simulation
Note that imposing periodic boundary conditions in time can be viewed formally as studying the system at a finite temperature that is inversely proportional to the time period t tot . Certain phase transitions may disappear when the temperature increases and the time period therefore decreases. However, in previous computer simulations for t tot = 4, 6 we found no indications that the presence of the B-C b transition depends on t tot [42]. Also for the minimal time extension t tot = 2 used here we still see a clear transition signal. By way of illustration, Fig. 9 shows the measurements of two different order parameters at the B-C b transition, for N 4 = 10k kept fixed. One of them is the order of the highest-order vertex in the triangulation T , where "order" is defined here as the number of one-dimensional edges sharing the vertex, normalized to lie between 0 and 1. 6 The other one is a normalized version of 5 More precisely, we work with t tot = 4, where the space-time geometry between t = 2 and t = 4 is an identical copy of the geometry between t = 0 and t = 2. This is done to maintain a regular triangulation, where by definition any (sub-)simplex is uniquely identified by its vertices, and happens purely for convenience of our computer code. 6 Many different definitions of "vertex order" and normalization are possible, leading to qualitatively similar results. The normalization the quantity conj( ) := N 41 − 6N 0 introduced at the beginning of Sect. 3. As also discussed in Sect. 3, at fixed N 4 and large t tot one finds a double-peak structure in the probability distribution of the order parameter x = N 41 − N 32 , superficially resembling the behaviour encountered at a first-order transition. Our observations for small t tot are entirely compatible with this picture, in the sense that the order parameters depicted in Fig. 9 also display a typical first-order behaviour, jumping back and forth between two different states on either side of the transition.
The B-C b transition appears when we keep κ 0 fixed (and not too large) and, coming from inside C b , decrease the coupling . Its pseudo-critical value c (N 4 ) is a function of the system size N 4 . By studying its behaviour as a function of N 4 we have found a dependence which can be fitted well to the functional form with some non-vanishing constant c and an exponent γ ≈ 2.4 that within measuring accuracy agrees with the corresponding exponent γ = 2.51(3) determined originally for a system with large time period [37,38]. In a similar vein, one can compare the behaviour of order parameters away from the B-C b transition, into phase C b and beyond, by increasing for fixed κ 0 . As an example, Fig. 10 shows the behaviour of the order parameter O P 1 , defined as the absolute value of the difference of the average Footnote 6 continued chosen here is a division by the maximal number of edges that could meet at a vertex in a triangulation that has the same numbers of vertices in the two spatial slices as the given triangulation T . This theoretical maximum would entail that the vertex is connected by an edge to every other vertex in the same spatial slice and to every vertex in the neighbouring spatial slice. Fig. 10 The order parameter O P 1 as a function of the coupling , measured at t tot = 2, N 4 = 10k and κ 0 = 2.2, indicating the presence of a phase transition between the de Sitter and bifurcation phases spatial curvatures of two adjacent spatial slices, where N 0 (t) and N 3 (t) denote the numbers of vertices and spatial tetrahedra contained in the spatial triangulation at time t. This quantity is one of several order parameters first introduced in [40] to study the newly discovered C d S -C b phase transition. The data points shown in Fig. 10, measured at t tot = 2, are qualitatively very similar to measurements of the same quantity for large t tot [40,41]. 7 This holds for the entire range of ∈ [0, 0.6] considered here, with the C d S -C b phase transition presumably located around = 0.2. For the volume N 4 = 10k used presently, the B-C b transition lies at = −0.042 (2) and therefore well outside the measurement range of Fig. 10. Note that the -values in the two-slice system with t tot = 2 are systematically lower than those of the system with full time extension t tot = 80 of [37,38], including for the extrapolated critical value c (∞) of the B-C b transition. Comparing with the results of [40,41], where the order of the C d S -C b transition is analyzed in more detail, the same seems to be true for this transition also. This is not surprising, since the systems are genuinely different and the location of a critical point is not a universal quantity. We conclude that our simulations with t tot = 2 reproduce the same characteristics of the bifurcation phase C b and the adjacent 7 Another difference is that in previous work [40,41]  phase transitions as were already seen for the large-time system with t tot = 80. The two-slice system therefore seems well suited for a further investigation of this phase.

Singular vertices
A key feature of the bifurcation phase, already reported in [40], is the appearance of a single "singular" vertex 8 of very high coordination number (this is the number n c (v) of foursimplices sharing a vertex v) on every second spatial slice. Coming from the de Sitter phase and moving into the bifurcation phase by lowering , one finds that a gap opens between the coordination number of the vertex with largest n c and that of the vertex with the second-largest n c . Well inside phase C b , the maximal n c (v) in a spatial slice containing such a singular vertex is typically orders of magnitude bigger than the average coordination number in the slice. At the same time, such a vertex is also singular from a purely three-dimensional point of view, in the sense that it is also shared by an exceptionally large number of spatial tetrahedra inside the spatial slice where it is located. Another observation, made in [40], is that in simulations with large t tot and therefore many spatial slices, the singular vertices on alternating slices are associated with a four-dimensional substructure of the triangulation, which takes the form of a chain of "diamond-shaped" regions in the time direction. This substructure is embedded in the rest of the triangulation and contains a large, finite fraction of the triangulation's total four-volume.
As already remarked in [40], the presence in the bifurcation phase C b of singular vertices and the structures associated with them breaks the homogeneity and isotropy (on average) of geometry which is present in the de Sitter phase C d S . Given the way Causal Dynamical Triangulations are implemented, there is nothing in principle that prevents homogeneity and isotropy of the average universe modelled by CDT triangulations, in the limit as the lattice spacing is taken to zero. This is indeed what is observed in phase C d S , where a number of properties of the dynamically generated "quantum universe" are very well described by a minisuperspace model with built-in spatial homogeneity and isotropy [26,27,50]. Moreover, in C d S the average shape of the universe can be fitted to a de Sitter space, a maximally symmetric space-time solving the classical Einstein equations. The appearance of isolated vertices of very high coordination number in phase C b is clearly incompatible with these symmetries. Given that phase transitions in physical systems are often related to the breaking of a symmetry, it is natural to associate the C d S -C b phase transition with a symmetry breaking also, namely, of homogeneity and isotropy.

Singular vertices cause bifurcation split
In the following, we will provide further evidence that phase C b is associated with the appearance of singular vertices and that they can be viewed as the decisive characteristic of the bifurcation phase. More specifically, we will establish a quantitative relation between the "bifurcation split", the observed typical volume difference between neighbouring spatial slices [40,41], and the order of the singular vertex present. We will set = 0, which for the volumes considered places us in the bifurcation phase, and at a safe distance from either of the adjoining phase transitions.
To analyze the geometry of the triangulations with t tot = 2 in greater detail, we will use a variant of the notion of vertex order, which for a given vertex v counts the number of (4, 1)-simplices between the two slices that share the vertex v and have a spatial three-simplex in common with the spatial slice not containing v. Using this definition, 9 we will call O max the maximal vertex order occurring in a given twoslice configuration. When a singular vertex is present, O max will coincide with the order of this vertex. Like in our earlier discussion of the matrix elements (19) of the reduced transfer matrix, we will use the letters m and n to denote the threevolumes of the two adjacent spatial slices. In addition, by definition, m will denote the volume of the slice that contains the vertex of maximal order, and n the volume of the slice that does not. Note that if a singular vertex v s is present in the spatial slice of volume m, O max ≤ n is the three-volume of the intersection of the (half-)diamond with tip v s and the spatial slice of volume n.
In Fig. 11 we show the distribution of the highest vertex order O max versus the volume difference n − m of the two spatial slices. One can roughly distinguish two regions. Below O max ≈ 300, the configurations contain no singular vertex in the sense that there is no significant gap between O max and the orders of the other vertices. A closer analysis reveals that for fixed O max in this region, the distribution of the volume differences is approximately Gaussian around n−m = 0. In other words, neighbouring slices preferentially have equal volumes. From previous investigations [39] we recognize the latter property as characteristic for configurations inside the de Sitter phase C d S . These configurations by no means dominate the dynamics of the bifurcation phase studied here, but the system makes occasional excursions to them, at least for the space-time volume we are considering. This will be further corroborated by data presented below.
The vast majority of configurations lie in the region where O max 400. Around O max = 400 a gap opens between O max and the distribution of the orders of the remaining vertices that becomes larger as the value of O max increases, signalling the appearance of a singular vertex. At the same time, at fixed O max , the configurations are now peaked around a nonvanishing volume difference. 10 This is typical for the bifurcation phase C b , where the effective transfer matrix n|M|m has a double-peak structure as a function of the volume difference n−m (and at fixed m+n), unlike the single peak found Fig. 12 Expectation value of the highest vertex order O max as a function of the difference n −m of the spatial-slice volumes (same specifications as in Fig. 11) in C d S . It entails that the two-slice volumes preferentially differ by a finite amount |n − m| = 0.
The interesting new finding from our data is that the expectation value O max depends linearly on this "bifurcation split" n − m between the two spatial volumes, where again the slice with the lower volume m is the one containing the highest-order vertex. This linear relation is illustrated in Fig. 12. Extrapolating n −m down to zero one obtains a vertex order of around 400, in agreement with Fig. 11. We conclude that the bifurcation phenomenon, observed in previous studies of the effective transfer matrix [40], seems to be a function of the appearance of singular vertices.
The particular choice of coupling constants for which the above results have been obtained is associated with specific expectation values for both the highest vertex order O max and the bifurcation split n −m. Not surprisingly, these variables have Gaussian-like distributions around their mean values. For example, O max has an approximate Gaussian distribution peaked at 675 with standard deviation around 50. Although it is not very visible on the scatter plot of Fig. 11, there are therefore many fewer configurations with vertex order 500 or 900, say, than there are with vertex order 700. Furthermore, for each given value of O max the width of the (Gaussian) distribution of n−m is approximately the same and coincides with the one determined by the effective action associated with the effective transfer matrix. This implies that the width is not a function of the vertex order for fixed values of the coupling constants.
The much rarer configurations with O max 400 have a special status, a fact that becomes clear when studying the maximal vertex order as a function of Monte Carlo time. As shown in Fig. 13, O max fluctuates around 675. Since there is a gap in the vertex order distribution below the maximal value, and since vertex orders can only change by relatively small amounts in each Monte Carlo update, the highest-order vertex usually remains located firmly in one of the two spatial slices. However, occasionally O max takes a very fast dip to a value below 500, which means that the distinguished, singular vertex disappears. After such a dip, a new singular, highest-order vertex appears randomly on either one of the spatial slices. We do not yet understand in detail how this process works, but the excursions occur seldom and their durations are much too short in Monte Carlo time to be explained as random processes associated with the Gaussian distribution of O max . The configurations with O max ≤ 500 in Fig. 11 constitute less than 0.1% of the total number of configurations.
Finally, we would like to understand whether there is just one singular vertex in a given spatial slice or whether further vertices with exceptionally high order can appear in the same slice when the system size goes to infinity. Studying O max as a function of the total space-time volume N 4 , we found that the relative order of the singular vertex, that is, O max divided by N 41 , grows with N 4 , as shown in Fig. 14. Recall that O max coincides with the four-dimensional volume of the (half-)diamond whose tip is the singular vertex, which in turn is bounded by N 41 . Since the measured ratio O max /N 41 is a finite fraction of 1, there can be at most a finite number of similarly "singular" vertices in the limit N 4 → ∞. Presumably it is just the single singular vertex on every second spatial slice we see at lower volumes. However, the detailed interpretation of this infinite-volume limit requires some care. The point is that just taking N 4 → ∞ at fixed coupling constant corresponds to changing the real, effective coupling constant. The presence of such a volume dependence is apparent from Eq. (22), which describes how the pseudo-critical c (N 4 ) increases with increasing N 4 . In the case at hand, for sufficiently large volume N 4 (larger than what we have considered here), our present choice = 0 will therefore no longer lie in phase C b , but in phase B. With this caveat in mind, our data indicate that in the infinite-volume limit, the CDT ensemble with t tot = 2 contains just one singular vertex.

Summary and conclusion
In this paper we have investigated the bifurcation phase C b recently discovered in CDT quantum gravity. We first reexamined the B-C b phase transition (formerly called the B-C transition). The order parameter used previously to determine the order of this transition exhibited an unexpected dependence on how the total space-time volume was fixed in the simulations: keeping the total number N 4 of foursimplices fixed resulted in a double-peak distribution for the order parameter, whereas keeping the number N 41 of foursimplices of type (4,1) fixed yielded only a single peak. A careful examination of the entropy factor N (N 41 , N 32 , N 0 ) revealed that in the volume range considered it has a rather complicated form as a function of N 41 and N 32 , which completely explains the observed behaviours of the order parameter for the two different volume fixings. These findings reconfirm that the double-peak structure seen for N 4 = const in no way contradicts the earlier conclusion that the B-C b phase transition is of second order [37,38].
The fact that the new C d S -C b transition was discovered in simulations with a short total time extension t tot , to determine the so-called effective transfer matrix, raised the question of whether the choice of t tot (as a long or short compactified time direction) has an influence on the observed phase structure. In the measurements presented above we have not found any indication that this is the case. The B-C b transition is still present for the system with t tot = 2, with a signal compatible with that observed for t tot = 80. Earlier work [40,41] had already shown that the new C d S -C b transition between the de Sitter and the bifurcation phase is also present for large t tot , and clearly visible for appropriate choices of order parameters. There are preliminary indications that this transition could be of higher order too [41], but more extensive simulations are needed to obtain more conclusive results.
The equivalence between long and short t tot motivated our further study of the properties of the bifurcation phase by considering the somewhat simpler two-slice system. We showed that the behaviour of the highest-order, "singular" vertex that appears in this phase is directly related to the previously observed tendency of the neighbouring spatial slices to develop a non-vanishing mean volume difference or "bifurcation split". More specifically, the maximal vertex order O max scales linearly with this volume difference. This gives us a more detailed, geometrical understanding of the mechanism behind the bifurcation split: a finite fraction of the (4, 1)-simplices between the two spatial slices clusters into a half-diamond whose tip is the singular vertex. This halfdiamond forms a substructure, which is embedded in the rest of the triangulation and leads to a corresponding "excess" of three-volume of the slice not containing the singular vertex.
At the same time, the appearance of a singular vertex 11 when crossing into phase C b from the de Sitter phase signals a breaking of the homogeneity and isotropy of geometry present in the de Sitter phase on scales above the cutoff scale. It suggests that the bifurcation-de Sitter phase transition can be associated with the breaking of a symmetry, a situation common in non-geometric statistical systems.
From this point of view the C d S -C b phase transition resembles the phase transition between the branchedpolymer and the crumpled phase in (Euclidean) Dynamical Triangulations. The DT configurations in the branchedpolymer phase appear to be homogeneous and isotropic (although not in any sense that is associated with a fourdimensional space-time), while configurations in the crumpled phase are characterized by the appearance of two distinguished, singular vertices of very high order and a singular link in between them [47][48][49]. Unfortunately, in this purely Euclidean quantum gravity model the phase transition between the two phases is only a first-order transition, even in extended DT models with an additional coupling constant, as already mentioned in the Introduction. In CDT we may be in the more exciting situation that the analogous C d S -C b phase transition is of second order, like the B-C b transition, and therefore it may be used to define a continuum theory of quantum gravity.
From a more general perspective, our investigation has given us additional insights into the type of mechanisms that can drive the nonperturbative dynamics of systems of (a priori) higher-dimensional geometry and the appearance of phase transitions, our understanding of which is rather limited. A conclusion we can already draw at this stage is that the phase structure of Causal Dynamical Triangulations in four dimensions, despite the presence of only two tuneable bare parameters, is amazingly rich and presents us with further opportunities to uncover viable continuum theories of quantum gravity.
tions g i , they can be determined uniquely from the data, without assuming any specific functional form for them. Also the constants c i get determined uniquely. Substituting the extracted functions and constants into (24) yields a function F(N 41 , N 32 ) that agrees with the directly measured F(N 41 , N 32 ) up to noise in the data.
The additional constraints we impose are where for each of the two variables z = N 41 , N 32 the average f (z) z of a function f (z) is defined as Taking into account the relations (25), three of the unknown quantities can be constructed directly from the measured function F(N 41 , N 32 ), namely, where the first equation involves a double average. We can find the remaining functions g 3 (N 41 ) and g 4 (N 32 ) by solving an eigenproblem. Let us define the two matrices One can show that in order to extremize (30), g 3 (N 41 ) must be an eigenvector of the matrix T , and g 4 (N 32 ) an eigenvector of the matrix T . To minimize E, one must choose the eigenvectors corresponding to the largest eigenvalues, a condition that fixes g 3 (N 41 ) and g 4 (N 32 ). The largest eigenvalue is positive and has the same value c 2 2 for both matrices, which also fixes the (positive) constant c 2 .
Having determined all functions g i and constants c i , the resulting model function F (N 41 , N 32 ) differs from the original empirical function F (N 41 , N 32 ) only by what looks like noise, which suggests that in the selected region the functional form assumed in (24) is very accurate.
The extracted functions g 1 (N 41 ) and g 2 (N 32 ) are shown in Fig. 15. Although g 1 (N 41 ) is approximated very well by  (24), extracted from the measured free energy F(N 41 , N 32 ) as described in the text. Again, the ranges of the parameters N 41 and N 32 on the horizontal axis have been rescaled and shifted. The red curves are best fit quadratic functions a quadratic function in a neighbourhood around its minimum, that is, in the range of N 41 we have been considering, this cannot possibly be true for its entire range. The reason is that a quadratic dependence would imply an entropy growth proportional to e +const·N 2 41 , which would contradict the fact, proven in [51], that the number of triangulations can grow at most exponentially with N 41 . The quadratic fit displayed in Fig. 15 must therefore be a local approximation arising as an expansion of a slower growing function of N 41 .
The extracted functions g 3 (N 41 ) and g 4 (N 32 ) are shown in Fig. 16. Because both g 2 (N 32 ) and g 4 (N 32 ) are well approximated by quadratic polynomials, the corresponding distri-butionP(x) for fixed N 41 (cf. Eq. (17)) is almost Gaussian. The function g 3 (N 41 ) is nearly linear so g 3 (N 41 ) · g 4 (N 32 ) results in a decreasing width of the distribution P(N 41 , N 32 ) as N 41 grows (the larger N 41 , the more negative is the coefficient in front of N 2 32 in F (N 41 , N 32 )). Going back to our earlier Fig. 6 depicting the distributionP(x) of x for fixed N 4 , the green curve is based on the ansatz (24) with quadratic functions g i . 12 It fits the data based directly on the measured free energy F(N 41 , N 32 ) (indicated by the blue dots) quite well.
To summarize, we have demonstrated that a simple ansatz like (24), with quadratic functions g i (x) can reproduce the observed features of the probability distributionsP(x) and P(x), without at the same time reproducing any signal of critical behaviour. This further supports earlier assertions that the appearance of a double peak inP(x) is not necessarily related to any specific form of critical behaviour.