The higher-order phase transition in toroidal CDT

We investigate the transition between the phases $B$ and $C_b$ observed in four-dimensional Causal Dynamical Triangulations (CDT). We find that the critical properties of CDT with toroidal spatial topology are the same as earlier observed in spherical spatial topology where the $B-C_b$ transition was found to be higher-order. This may have important consequences for the existence of the continuum limit of CDT, describing the perspective UV limit of quantum gravity, which potentially can be investigated in the toroidal model.


Introduction
Numerical Monte Carlo simulations applied to lattice field theories became an important tool of contemporary physics. The famous example is Lattice Quantum Chromodynamics (QCD) which has grown up from its childhood and now goes hand-by-hand with experiments and beyond, e.g. by investigating very non-trivial QCD phase diagram in the regime of coupling constants non-tractable by perturbative calculus. Despite many open questions, QCD has a well defined ultraviolet limit, where it becomes non-interacting asymptotically free theory and thus the high energy behaviour can be investigated perturbatively. The opposite thing happens when one tries to formulate a quantum theory of gravity (QG) by applying standard quantum field theory techniques to Einstein's General Relativity (GR). In that case the perturbative expansion around any fixed classical metric field fails at high energies due to the perturbative non-renormalizability of such a formulation [1,2]. However, as conjectured by Steven Weinberg in his seminal paper [3], QG can be asymptotically safe, i.e. it can admit a well behaved non-perturbative high energy limit defined in the vicinity of a non-trivial fixed point of the renormalization group flow, where quantum gravity becomes scale-invariant and thus can be extrapolated to arbitrarily large energy scale. If the asymptotic safety scenario is valid 1 then (in the ultraviolet regime) QG must be formulated in a background-independent non-perturbative way making lattice approaches well suited to tackle this problem. In such formulations one discretizes geometric degrees of freedom on the lattice with (4-dimensional) lattice 'volume' N 4 and with a minimal (cut-off) spacing a, and in the ultraviolet regime one would like to get rid of the discretization by taking a continuum limit of a → 0 and N 4 → ∞ such that N 1/4 4 · a is related to some physical length. In order to obtain non-trivial physical observables in the continuum limit, where a → 0 and N 4 → ∞, one would also like to have appropriately divergent correlation lengths c ∼ N 1/4 4 . Thus in a lattice approach the continuum limit should be associated with a higher order (continuous) phase transition. Therefore studies of the phase structure and orders of phase transitions are important steps towards defining an ultraviolet limit in a lattice formulation and thus testing the asymptotic safety scenario for gravity.
One of the most successful attempts of the lattice formulation of quantum gravity is that of Causal Dynamical Triangulations (CDT) (for reviews see [10,11]), in the sense that it has a rich phase structure, where some of the transitions are higher order, which potentially can be used to define continuum limit and that it additionally has a well behaved low energy limit consistent with GR. CDT is based on the path integral formalism and makes only a few assumptions on the geometry of quantum space-time, namely it requires that the geometry can be globally foliated into space-like hypersurfaces, each with the same fixed topology Σ. The model is using the discretization of space-time following the method proposed by Regge [12]. The three-dimensional spacial states are constructed by gluing together in all possible ways regular tetrahedra with a common link length a s to form a triangulation of a three-dimensional space with a (closed) topology Σ. The topology of states is fixed during the evolution of geometry in time, being the origin of the name causality in the model. To join states at different times t we need two types of 4-dimensional simplices. Tetrahedra become bases of 4-dimensional simplices {4, 1} (and {1, 4}) with four vertices at a time layer t and one at t+1 (resp. t−1). In our notation the simplex {i, j} has i vertices at a time t and j vertices at a time t + 1. The time links are assumed to have a common link length a t which may be different than a s . To complete the manifold structure two additional simplex structures are necessary. These are {3, 2} and {2, 3} simplices. Pairs of simplices share a common three-dimensional face (tetrahedron). The construction works both for systems with Lorentzian signature and, after Wick rotation, for systems with Euclidean signature. Each space-time configuration can be interpreted as Lorentzian or Euclidean. The possibility of performing Wick rotation is crucial if we want to use numerical methods to analyze the properties of the model. In the following, we assume the Euclidean formulation is used. The discretization described above means that the four-dimensional volume of all {i, j} simplices depends only on the type of a simplex. Similarly other geometric properties, like the angles, are universal for all simplices of a particular type.
The studied object is the Feynman amplitude Z, which is expressed as a weighted sum over manifolds T joining the initial and final geometric states separated by time T . The weight is assumed to be expressed as a discretized version of the Hilbert-Einstein action S EH (T ) where C(T ) is the symmetry factor of a graph representing the manifold. In practice the choice of the initial and final states is replaced by assuming the system to be periodic with the period T . The discretized version of the Hilbert-Einstein action takes the form [13] where N i,j denotes the number of 4-dimensional simplicial building blocks with i vertices on hypersurface t and j vertices on hypersurface t ± 1, and N 0 is the number of vertices in the triangulation. κ 0 , ∆ and κ 4 are bare coupling constants. κ 0 and κ 4 are related to Newton's constant and the cosmological constant, respectively, and ∆ depends on the ratio of the length of space-like and time-like links in the lattice. In the Monte Carlo simulations of CDT the parameter κ 4 , which is proportional to the cosmological constant, is tuned such that one can take infinite-volume limit. As will be explained later, in numerical simulations we perform a series of measurements for systems with increasing (fixed) volume N 4,1 and try to determine the limiting behaviour for N 4,1 → ∞. In the consequence the phase diagram presented in Figure 1 depends only on two bare couplings κ 0 and ∆. It is remarkable that such a simple model has a rich phase structure with four phases with very different physical properties. The analysis of the phase structure and, in particular, the order of phase transitions is fundamental to relate the model to a possible theory of quantum gravity.

Phase transitions in MC simulations of lattice field theories
According to Ehrenfest's classification, the order of a phase transition depends on the behaviour of the thermodynamic free energy. If all first n − 1 order derivatives of the free energy are continuous functions of some thermodynamic variable, e.g. the coupling constant of the lattice theory, and the n-th order derivative exhibits a discontinuity at the transition point then the transition is the n-th order phase transition. Here we are especially interested to distinguish between the first-and the higher-order phase transitions, as the continuous limit of the lattice field theory should be associated with the latter type. The derivatives of free energy are related to order parameters, which capture differences of thermodynamic properties of the system in two different phases separated by the transition point. For a first-order transition one should observe a discontinuity of the order parameter at the transition point and for the higher-order transition the order parameter should be continuous but its derivatives, e.g. its susceptibility, should diverge. Unfortunately measuring the (dis-)continuity of the (derivatives) of an order parameter in numerical simulations is a tedious task. Actually, in numerical Monte Carlo simulations, which are always performed for a finite lattice size N 4 , one does not even observe phase transitions per se. The finite lattice size and the finite lattice spacing make all thermodynamic functions and their derivatives finite, even though they can become arbitrarily large for large lattice sizes. One should therefore carefully analyze finite (lattice) size effects and, if possible, take the infinite (lattice) volume limit N 4 → ∞.
As phase transitions are usually related to breaking some symmetries of the studied lattice field theory, one can define order parameter(s) OP which capture these symmetry differences between various phases of the theory in question. One then usually performs numerical Monte Carlo (MC) simulations for some fixed lattice volume N 4 in many points of the theory parameter space (see e.g. the CDT phase diagram in Figure 1) to find regions where the order parameter rapidly changes, see e.g. Figure 4 where we show the mean value OP of the four order parameters (for their definitions see equation (9)) used in CDT phase transition studies measured in the B − C b transition region. The precise position of the phase transition is signaled by a peak of the susceptibility of an order parameter related to its first-order derivative with respect to some thermodynamic variable, see e.g. Figure 5. For a finite lattice volume N 4 one can only determine a position of the (volume dependent) pseudocritical point. Positions of such points may in general depend on the order parameter or the method used. Only in N 4 → ∞ limit they must coincide. Let ∆ c (N 4 ) be the pseudo-critical value of the thermodynamic variable ∆, e.g. the coupling constant, measured for a given phase transition for the lattice volume N 4 . The typical (large) volume dependence is where the critical exponent ν is one for a first-order transition and larger than one for a higherorder transition. Thus by making a series of measurements of ∆ c (N 4 ) for different lattice volumes N 4 one can establish a value of the critical exponent ν and in effect determine the order of the phase transition. Another way of distinguishing between the first-and the higher-order phase transitions in numerical Monte Carlo studies is to analyze the behaviour of the order parameter(s) measured precisely at (or in practice as close as possible to) the transition point. For a first-order transition the discontinuity of an order parameter can appear in its MC history as jumps between two different states. In such a case, the histogram of the order parameter measured at the pseudo-critical point should show two separate peaks centered around the values generic for the two different phases.
Here one should also carefully analyze finite size effects related to the finite lattice volume N 4 fixed in the numerical studies. The separation of the peaks in the MC history histogram can either increase or decrease with the lattice volume which can imply the first-or the higher-order transition, respectively. If the separation of the states, generic for the first-order transition, is large enough one typically observes a hysteresis at the transition region. In order to check that, one can run two separate series of Monte Carlo simulations, one initiated with configurations generic for one phase and the other one initiated with configurations generic for the other phase. If hysteresis is present then one can observe a (statistically) different behaviour of the two series in the transition region, e.g. the pseudo-critical points measured in the two different series could be shifted versus each other. If hysteresis is absent the results of the two series should (statistically) agree. Running two independent series initiated with different staring configurations is also a good way of checking thermalization of the Monte Carlo data, i.e. checking if the MC simulation has run for long enough to reach the proper statistical equilibrium and thus if measurement data can be collected, see e.g. Figure 2.  (9)) measured in two independent MC simulations of CDT with toroidal spatial topology with exactly the same parameters, i.e. N 4,1 = 300k, T = 4, κ 0 = 2.2, ∆ = 0.048. One simulation was initiated with a configuration from phase B (blue line) and the other one started from a configuration from phase C b (orange line). Both data series statistically agree from ca 40000 sweeps (1 sweep = 10 7 attempted MC moves). Data from earlier MC time history, called the thermalization period, are excluded from final measurements.

Another quantity of interest is the Binder cumulant
which is always non-positive because (OP 2 ) 2 − OP 2 2 ≥ 0, and it reaches a minimum at the pseudo-critical point ∆ c (N 4 ), because there fluctuations are maximal. In the numerical MC simulations one can measure the (volume dependent) value of the Binder cumulant minimum for different (fixed) lattice sizes N 4 and then analyze its behaviour in the large volume limit N 4 → ∞. In the case of a higher-order phase transition the probability distribution of the order parameter OP approaches a Dirac delta around OP in the infinite volume limit. And then 2 Note that here we use a definition of the Binder cumulant which is shifted (by a −2/3 constant) versus the original Binder's formulation [14,15,16]: x 4 x 2 2 . The definition (5) was used in previous CDT phase transition studies [17,18,19] and thus we keep it in order to ease comparison with these results. The virtue of using our definition is that, as explained in the text, the deviation of (critical) B OP from zero with rising lattice volume may signal a first order transition, while the convergence to zero is characteristic of a higher order transition. One could as well use the original Binder's definition and look at the deviation from 2/3. B min OP (∞) should equal 0. In the case of the first-order transition the distribution of the parameter OP is a sum of two distributions centered at expectation values characteristic for the two different phases. In the infinite volume limit, when these distributions approach Dirac delta functions, the minimum of the Binder cumulant becomes: where OP B and OP C b are expectation values of the observable OP at two different phases, say "B" and "C b ", and the relative strength of Dirac delta functions is assumed to be 2 , respectively. In Table 1 we summarize methods used in numerical MC simulations of lattice field theories to distinguish between the first-and the higher-order phase transitions. We will then apply these methods in Section 4 to analyze the B − C b transition in CDT with the toroidal topology of spacial slices.

OBSERVABLE
First-order transition Higher-order transition Critical exponent ν in ν ν scaling of ∆ c (N 4 ), eq. (4) = 1 > 1 OP histograms measured at double peaks single peak or pseudo-critical points ∆ c ( 3 The properties of the bifurcation phase C b The existence of the bifurcation phase in the CDT model with a spherical spatial topology was discovered relatively late [20,21,22]. The reason why in the early studies only three phases were discussed was that the basic observable used in these approaches was the (average) spatial volume profile of configurations. A typical setup for numerical experiments was to use systems periodic in time, with a period T usually in the range 40 -80. Using the spatial volume observable, the three phases, A, B and C, were characterized by completely different qualitative behavior. The phase A was characterized by large fluctuations of the spatial volume in the neighboring time slices. The observed average volume distribution in time corresponded to the unbroken symmetry of the time translations. In the phase B almost all spatial volume (except for the stalk, necessary to satisfy the periodic boundary conditions) was concentrated at a single time slice. This meant that for typical states in this phase the symmetry of the time translations was fully broken. The physically most interesting was the phase C, where the volume profile contained the blob and the stalk, again meaning that for a typical configuration the symmetry of the time translations was broken. Average volume distribution in the blob and its fluctuations could be very accurately explained using the effective mini-superspace model for the isotropic four-dimensional Euclidean space-time [23,24,25]. Most results were obtained for a particular point in the coupling constant space with κ 0 = 2.2 and ∆ = 0.6, where it was shown that volume distribution scaled with the total N 4,1 lattice volume in a way consistent with the Hausdorff dimension d H = 4. Similar measurements performed for decreasing values of ∆ showed that, although qualitatively the volume profile still contained a blob and the stalk, the scaling properties did not follow those determined in the de Sitter phase C. It was observed that the scaling was consistent with that predicted for systems with the Hausdorff dimension d H = ∞. The name bifurcation phase C b appeared to describe the additional property observed in the volume profile: a different behavior in the even and odd time slices when the time period T was sufficiently small [20]. It was soon realized that the reason for the observed behavior came from the breaking of the isotropy of the spatial volume distribution in the new phase. For the time slices separated by two units in time, vertices with very high coordination numbers appeared, leading to a formation of highly nontrivial geometric objects, forming a chain in the time direction. A physical interpretation of these objects was conjectured to be a result of a local signature change from Euclidean to Lorentzian [21], producing objects with some qualitative similarity to a black hole or rather a series of black points. A detailed description of the microscopic mechanism producing such effects will be the subject of a separate paper.
As can be seen in Figure 1, for decreasing values of ∆ and a fixed value of κ 0 , one observes a phase transition between the C b and B phases. The properties of this phase transition were very accurately measured in the case of a spherical spatial topology [17,18,22], although originally the phase C b was interpreted as being a part of the de Sitter phase C. Results indicated that the phase transition was higher order, a very important property from a theoretical point of view, as explained earlier. The purpose of the present analysis is to check if the position and properties of the phase transition remain the same for systems with the spatial topology Σ of a sphere S 3 and of a three-torus T 3 .
The first question to be asked is: are the qualitative properties in the C b phase similar or different when we consider systems with a different spatial topology. Again we may look at the simplest object, a volume profile for systems with the periodicity T of the same order as the one used in the spherical case. This is the observable which was found to behave differently in the C phase. The observed volume profile, in this case, was found to be flat rather than containing a blob [26,27]. The reason of such a behavior could be explained using a mini-superspace spatially isotropic model for a system with the spatial topology of a three-torus. The averaged volume profile is flat since in the toroidal case the time translation symmetry remains unbroken [26,27].
Investigations show that this is not the case in the bifurcation phase C b . The volume profile observed for the point in the coupling constant space, typical for the bifurcation phase (κ 0 = 2.0 and ∆ = 0.2) shows the appearance of a blob and the stalk, see Figure 3, the same way as it was observed in the spherical case. Also the scaling of the volume profile with the total N 4,1 lattice volume is consistent with the Hausdorff dimension d H = ∞, the same as in the spherical CDT. The analysis of the geometric properties of configurations in the bifurcation phase C b shows that also from a microscopic point of view the toroidal and spherical cases are very similar. In both topologies, we observe the high-order vertices, separated in time by two steps. The shape of the blob observed for periodicity T large enough (T ≥ 20) again scales consistently with the infinite Hausdorff dimension. The difference is observed in the stalk, which has a much larger volume for a torus than that for a sphere. This is well understood and results from the fact that a minimal 3D spatial configuration depends strongly on the topology (see [26]).
As a conclusion, one may expect the critical properties of the phase transition between the C b and B phases to be very similar in both topologically different realizations of the model. Below we show that this is indeed the case. The measurement of the critical behavior on the boundary between C b and C phases may, on the other hand, be different, or at least difficult to be determined numerically.

The B − C b phase transition in the toroidal CDT
Below, we present the results of the B −C b phase transition study in CDT with the toroidal spatial topology. The B − C b transition was earlier studied in the spherical spatial topology [17,18,22] where it was classified to be the higher order transition. As explained in Section 2 in order to investigate the phase transition one has to make a series of Monte Carlo simulations for various points in the CDT (κ 0 , ∆) parameter space 3 , around the phase transition point. In this study all measurements were taken for one fixed value of κ 0 = 2.2 and for a sequence of ∆ values. 4 In each simulation the N 4,1 lattice volume of the system (i.e. the total number of {4, 1} and {1, 4} simplices) is fixed or, more precisely, it fluctuates around the target valueN 4,1 . The lattice volume is controlled by a volume-fixing potential δV = (N 4,1 −N 4,1 ) 2 (8) added to the bare Einstein-Hilbert-Regge action of CDT (2) such that the volume is sharply peaked around a chosen value ofN 4,1 , with a well-defined amplitude of fluctuations controlled by the numerical constant. In the CDT Monte Carlo simulations one also has to set the length of the (periodic) time axis, i.e. the number of (integer) time slices T . In our case the number of time slices was equal T = 4 and the magnitude of volume fluctuations was fixed at = 0.00002. 5 In our analysis we will focus on the behaviour of four order parameters which have previously been successfully used in phase transition studies both in the spherical [18,28,29] and the toroidal [30,19,31] CDT, 6 where V 3 (t) is the spatial volume 7 in the time slice t and O(v) is the vertex coordination number, i.e. the number of simplices sharing a given vertex v. The behaviour of the order parameters in all CDT phases has been summarized in Table 2. Specifically when changing from the phase B to the phase C b the OP 1 , OP 2 and OP 4 increase in value while the OP 3 decreases, see Figure 4. The MC simulations were performed for nine different (fixed) lattice volumes, i.e. forN 4,1 = 40k, 60k, 80k, 100k, 120k, 140k, 160k, 300k, 400k. For each lattice volumeN 4,1 the approximate location of the B − C b phase transition point was found and then a series of precise measurement was performed for ∆ in the range around the expected critical value ∆ c with a resolution of 0.001. Each measurement series was performed twice, each time for a different initial triangulation: one from phase B and one from phase C b , and the two data series were compared in order to check thermalization and possible hysteresis, see e.g. Figure 2. For each lattice volumeN 4,1 and each of the two measurement series (s = B, C b ) and each of the four order parameters OP i (i = 1, 2, 3, 4) the precise position of the (volume dependent) pseudo-critical point ∆ c i,s (N 4,1 ) was established based on the peak of the OP i,s susceptibility χ OP i,s , see Figure 5 where we present the results of measurements for the lattice volumeN 4,1 = 100k. The values of ∆ c i,s (N 4,1 ) measured for different OP i and in the two data series in general coincide up to the used ∆ resolution. If the results for various OP i or for various data series are different, usually shifted not more than by the ∆ difference of 0.001, we simply take the arithmetic mean and assign a correspondingly larger measurement error, e.g. for the lattice volumeN 4,1 = 100k one has ∆ c (N 4,1 = 100k) = 0.0376 ± 0.0016.
advocated in [19] the order of CDT phase transitions does not depend on that. 6 Here we use a slightly different definition of OP 1 than in previous CDT phase transition studies, where it was: OP 1 ≡ N 0 /N 4 . Current definition is more natural when N 4,1 volume is fixed (see equation (8)) which was the case in all MC simulations described herein. 7 To ensure consistency with our earlier publications we define V 3 (t) as twice the number of spatial tetrahedra with the integer time coordinate t.  Then we fit the finite size scaling relation (4) to the measured ∆ c (N 4,1 ) values. The best fit of the true (infinite volume) critical point is ∆ c (∞) = 0.073 ± 0.004, and the best fit of the critical scaling exponent is ν = 2.7 ± 0.4 which supports the higher-order nature of the B − C b phase transition, see also Figure 6 where we plot the measured data together with the best fit of the scaling relation (4) and compare it to the fit with a forced value of ν = 1 (typical for a first-order transition) showing that the quality of the latter fit is much worse. The measured values of the true critical point and the critical exponent also agree with ∆ c (∞) = 0.077 ± 0.004 and ν = 2.51 ± 0.03 measured in CDT with the spherical spatial topology [18], giving strong evidence that the results are independent of the topology chosen (at least for the toroidal and the spherical one).  The left plot is for data series started from configuration in phase B and ∆ = 0.037 (i.e. the peak of susceptibility χ OP 1 measured for this data series, see Figure 5) while the right plot is for data series initiated in phase C b and ∆ = 0.038 (peak of χ OP 1 for this data series).
In order to corroborate this result, we have performed the detailed Monte Carlo history analysis of all order parameters at (and in the vicinity) of the measured pseudo-critical points, see Figure 7 where we plot the MC history histograms of the OP 1 measured for the example N 4,1 = 100k volume and for ∆ = 0.037 (peak of χ OP 1,B ) and ∆ = 0.038 (peak of χ OP 1,C b ). In none of the cases have we observed the double peaks in the measured histograms nor the hysteresis of the measured data series. These results support the higher-order B − C b transition.
Finally, we have analyzed the behaviour of the Binder cumulants (5) in search of minima, see Figure 8 where we plot data measured for N 4,1 = 100k. The value of pseudo-critical∆ c i,s (N 4,1 ) defined by the minimum of the Binder cumulants B OP i,s in general coincides with the ∆ c i,s (N 4,1 ) value defined by the maximum of susceptibility χ OP i,s , the possible shift is usually up to ∆ difference of 0.001. In Figure 9 we plot the measured values of B min OP i,s (N 4,1 ) ≡ B OP i,s (∆ c i,s (N 4,1 )) as the function of the lattice volume N 4,1 . 8 All Binder cumulants measured for OP 1 ,..., OP 4 visibly grow towards zero when N 4,1 is increased, which again favours the higher-order nature of the B − C b transition.

Summary and conclusions
Applying phase transition analysis methods described in Section 2 to the B −C B transition in CDT with the toroidal spatial topology we have shown that the transition is most likely the higher-order phase transition. This result is supported both by the finite size scaling analysis of equation (4) showing the best fit scaling exponent ν = 2.7 > 1, by the large volume behaviour of the Binder cumulant minima (6): B min OP i (N 4,1 → ∞) → 0 and by the lack of hysteresis / two-state jumping of the order parameters measured at the (pseudo) critical points.
The above result and also numerical values of the critical scaling exponent ν = 2.7±0.4 and the true critical point ∆ c (∞) = 0.073±0.004 are also consistent with the B −C b transition measured in CDT with the spherical spatial topology for the same fixed value of the κ 0 = 2.2 parameter, where ∆ c (∞) = 0.077±0.004 and ν = 2.51±0.03, respectively [18]. Thus the B −C b transition properties are the same in both spatial topologies. This is also the case for the A − C transition which was found to be the first-order phase transition in both topologies -the detailed analysis of the A − C transition in the spherical and the toroidal CDT for various Monte Carlo simulations' parameters (lattice volume fixing methods and lengths of the (integer) time period T ) can be found in [19]. One can therefore formulate a conjecture that CDT results including the phase structure and the order of phase transitions are independent of the spatial topology choice, which is a parameter put in "by hand". Each data point denotes OP 1 measured from last 100k sweeps (1 sweep = 10 7 attempted MC moves), data from initial thermalization period were skipped. Shaded regions between the dashed lines denote the range of the measured data. Hysteresis is clearly visible for ∆ ≥ 0.38, especially for the OP 4 parameter which is the most sensitive to the C − C b transition. This is not the case for the higher-order B − C b transition (described herein) observed around ∆ ≈ 0.05.
The question mark remains for the C − C b transition which was found to be the higher-order phase transition in the spherical CDT [28,29] 9 and has not been yet investigated in detail in the toroidal CDT. The reason is that in the toroidal CDT case one observes a very strong hysteresis in the C − C b transition region (see Figure 10) and therefore one is not able to perform precise MC measurements which would enable one to make finite size scaling analysis as it was explained in Section 2. The very strong hysteresis would suggest that the C − C b transition is most likely the first-order transition in the toroidal CDT, i.e. the order of the transition would change due to the different spatial topology. But this can be as well an algorithmic issue of the MC code used in the CDT simulations and more advanced methods should be used in order to resolve this problem. 10 In the toroidal CDT one was also able to make MC simulations in the most interesting region of the CDT parameter space, namely in the vicinity of the two "triple" points where the A − B − C and the B − C − C b phases meet (see the CDT phase diagram in Figure 1), which was not possible in the spherical CDT where MC simulations got effectively "frozen" in this region of the phase diagram. As a result in the toroidal CDT one observes the direct B − C transition which was classified to be the first-order transition, albeit with some atypical properties suggesting a possible higher-order transition [31]. Summing up, we have shown that the B − C b transition is the higher order transition which most likely makes the B − C − C b "triple" point the higher order transition point even though the B − C and the C − C b transitions are possibly the first-order transitions. The above "triple" point is thus a natural candidate for an UV fixed point for QG [34,35].