The phase structure of Causal Dynamical Triangulations with toroidal spatial topology

We investigate the impact of topology on the phase structure of four-dimensional Causal Dynamical Triangulations (CDT). Using numerical Monte Carlo simulations we study CDT with toroidal spatial topology. We confirm existence of all four distinct phases of quantum geometry earlier observed in CDT with spherical spatial topology. We plot the toroidal CDT phase diagram and find that it looks very similar to the case of the spherical spatial topology.


Introduction
Asymptotic safety provides an exciting possibility of formulating a nonperturbative and background independent theory of Quantum Gravity. The idea was put forward in a seminal paper of S. Weinberg [1] who proposed to extend the notion of renormalizability into the nonperturbative regime defined at non-Gaussian ultraviolet fixed point(s) (UVFP) of renormalization group trajectories. If the renormalization group flow originating in the UVFP lies (in an abstract space of coupling constants) on a hypersurface of finite dimension then only a finite number of (important) couplings are attracted to the UVFP. Such couplings could in principle be determined by a finite number of experiments making the whole formulation finite and predictive for arbitrarily large energy scale. As a result it might be possible to overcome the well-known problem of (perturbative) nonrenormalizability of gravity treated as a conventional quantum field theory (QFT) expanded around any fixed background geometry [2,3]. There are known examples of perturbatively nonrenormalizable but asymptotically safe QFTs [4,5,6], and for gravity the Weinberg's conjecture is strongly supported by functional renormalization group studies [7,8,9,10,11,12] which provide growing evidence that gravitational UVFPs really exist. A key difficulty remains: in order to investigate an asymptotically safe theory in full-glory one is forced to apply nonperturbative tools. Such tools are provided in the research program of Causal Dynamical Triangulations (CDT) [13] which is an attempt to quantize gravity based on a latice regularization of the nonperturbative path integral (1) over spacetime geometries [g], i.e. equivalence classes of metrics g with respect to diffeomorphisms, and S HE is the classical Einstein-Hilbert action. To give precise meaning to expression (1) CDT introduces a lattice regularization by constructing (continuous) spacetime geometries from four-dimensional simplicial building blocks. It is assumed that spacetime inside each such block is flat and the nontrivial geometry depends on how these blocks are glued together (e.g. a local deficit angle in D dimensions is encoded in the number of simplices sharing a given D-2 dimensional subsimplex). By gluing simplices together one obtains piecewise linear simplicial manifolds, also called triangulations, and the formal path integral (1) is defined as where the sum is over triangulations T , and S R is the Einstein-Hilbert action for a triangulation obtained following Regge's method for describing piecewise linear simplicial geometries [14].
A theory of quantum gravity should describe spacetime at the Planck scale, where one may expect large fluctuations of the geometry. Such fluctuations could in principle lead to changes in spatial topology. Therefore there is a longstanding discussion if in the path integral (2) one should allow for topology changes and, if so, which topologies should be taken into account. Including topological fluctuations was considered in the easiest case of two-dimensional toy models in the Euclidean formulation [15]. It was shown that in this case a naive understanding of the topology of the Universe breaks down and if we want to stay in a formalism similar to QFT we must suppress such fluctuations. The problem is that there are many more geometries of complicated topology than there are of simple topology and any sum over geometries is dominated by these complicated topologies and plainly divergent. 1 The problem is somehow eased in the Lorentzian formulation where one is able to define two-dimensional models with spatial topology fluctuations but at a cost of restricting the class a geometries taken into account to those causally 'well-behaved' [16,17,18,19]. In higher dimensions the situation becomes much worse, and in particular in four dimensions the problem is not even well posed as four-dimensional topologies are not classifiable. Therefore in CDT one adopts a pragmatic point of view by considering only spacetimes admitting a global proper time T foliation into spatial hypersurfaces Σ of fixed topology. This requirement is compatible with imposing an additional causal structure of global hypebolicity on admissible geometries (triangulations) which enter the path integral (2), such that each triangulation T is topologically M = Σ × T . The idea of an imposed global time-foliation also appears in Hořava-Lifshitz gravity [20], which is an attempt to define the UV completion of general relativity by introducing anisotropic scaling of space-time coordinates in the high-energy regime. It has been actually shown that the continuum limit of two-dimensional CDT is compatible with the two-dimensional projectable Hořava-Lifshitz gravity [21,22], and that both approaches share many features in three [23,24,25,26] and four [27,28,29] space-time dimensions.
In four-dimensional CDT, each spatial layer of integer (lattice) time t is constructed from equilateral tetrahedra with space-like links a s . The neighbouring spatial layers at t and t ± 1 are linked by four-(dimesional-)simplices, with additional time-like links a t , glued together in such a way that also all intermediate spatial layers between t and t ± 1 have the requested fixed topology Σ. One can show that the four-dimensional simplicial complex obeying the CDT topological restrictions can be constructed from just two types of building blocks (see Fig. 1), the (4, 1) simplex with 4 vertices (a tetrahedron) in time t and one vertex in t ± 1, and the (3, 2) simplex with 3 vertices (a triangle) in t and 2 vertices (a link) in t ± 1. In this case the Einstein-Hilbert-Regge action takes a form [30] S R = − (κ 0 + 6∆) N 0 + κ 4 N (4,1) + N (3,2) where N (4,1) , N (3,2) and N 0 denote the total number of (4, 1) simplices, (3, 2) simplices and vertices, respectively. The action includes three bare dimensionless coupling constants κ 0 , ∆ and κ 4 . κ 0 is inversely proportional to the bare Newton's constant, ∆ is related to the ratio of the length of space-like and time-like links on the lattice (a 2 t = −α a 2 s , where α > 0 is called the asymmetry parameter) and κ 4 is proportional to the bare cosmological constant. Such a formulation is coordinate-free (all geometric degrees of freedom are expressed by topological invariants -lengths and angles) and nonperturbative (all possible triangulations are included in the path integral (2)). It is also manifestly background independent as there is no preferred triangulation put in by hand but macroscopic geometry emerges dynamically from fluctuations. It is important to note that no ad hoc discreteness of spacetime is assumed from the outset, and the discretization appears only as a regularization of the path integral (1). The finite length a of the links in a triangulation constitutes an ultraviolet cutoff which is intended to be removed in the continuum limit a → 0, which should be consistent with the perspective UVFP of the renormalization group trajectories. In a lattice formulation, as CDT, the UVFP should appear as a phase transition of second or higher order, where infinite correlation lengths enable one to go simultaneously with the lattice spacing a → 0 and the linear lattice size L → ∞, such that physical lengths aL remain constant and thus observable quantities expressed in physical units are kept fixed. Therefore analysis of the CDT phase structure and the order of the phase transitions constitute first steps in a quest for the ultraviolet limit (in CDT formulation) of quantum gravity. At the same time it is very important to be able to correctly reproduce the infrared limit compatible with classical Einstein's general relativity.
This paper is organised as follows. After reviewing some technical details regarding the numerical implementation in section 2 we summarize the state of the art of CDT in section 3. In section 4 we define the order parameters and we explain how to study the CDT phase structure. In section 5 we present the results of numerical simulations performed with the spatial topology fixed to that of a three-torus. A discussion and a summary of the results obtained in this work are presented in section 6.

How to perform numerical simulations
The study of the regularized path integral (2) in four spacetime dimensions requires using numerical methods. This is possible by applying a Wick rotation of the proper time coordinate from Lorentzian (real) time t (L) to Euclidean (imaginary) time t (E) = −it (L) . Such a Wick rotation is well defined in CDT due to the assumed spacetime foliation into Cauchy hypersurfaces of constant proper time. The causal structure enables one to change time-like links into space-like links by changing the asymmetry parameter α → −α and accordingly the length of the time-like links becomes a 2 t = α a 2 s , α > 0. As a result the four-simplices become parts of Euclidean space. At the same time one should change the Lorentzian action S R by analytically continuing the bare coupling constants κ 0 , ∆ and κ 4 from eq. (3), which are all analytic functions of α in the lower half of the complex α plane. 2 Consequently the path integral of CDT (2) becomes a partition function of a statistical theory of random geometries One should keep in mind that the class of admissible (Euclidean) triangulations T (E) which enter the partition function (4) keeps track of the imposed Lorentzian structure by the proper time foliation constraint, and it is not the same as for Euclidean Dynamical Triangulations (EDT) where such a constraint is absent. EDT have not been able to correctly reproduce the suitable infrared limit [31,32,33,34,35] nor to define a continuum limit 3 [37,38], irrespective of the spacetime topology chosen [39]. The problem with EDT probably lies in the fact that the formulation is Euclidean from the outset and thus the distinction between space and time is lost and time has to be reconstructed dynamically, which leads to its pathological behaviour. These problem seems to be cured in CDT, where time is treated semi-classically. The above construction makes it possible to study the partition function (4) using numerical Monte Carlo techniques. One starts from an arbitrary simple initial triangulation (configuration of simplices), consistent with the requested fixed spacetime topology, and updates its geometry using Monte Carlo moves. The moves are local (they create, delete or reconstruct simplices in the closest neighbourhood of some randomly chosen place in a triangulation), causal (preserve fixed spacetime topology as well as the foliation structure) and ergodic (any triangulation obeying CTD topological restrictions is achievable from any other triangulation by a sequence of moves). By using the Metropolis algorithm the numerical code applies the moves in such a way that the system performs random walk in the space of admissible triangulations and, after a thermalization period, statistically independent triangulations T (E) are generated with probabilities consistent with Boltzman weights ∝ exp(−S . As a result one can use generated triangulations to estimate expectation values or correlators of observables. In a typical Monte Carlo simulation one sets the values of the bare coupling constants: κ 0 , ∆ and κ 4 . Simulations show that for fixed values of κ 0 and ∆ the number 2 The functional form of the bare Euclidean Einstein-Hilbert-Regge action S (E) R is the same as for the Lorentzian action S (L) R from eq. (3). Of course the form of κ 0 , ∆ and κ 4 as functions of α is modified, but anyway they are bare coupling constants and as such their functional dependence of the bare Newton's constant, cosmological constant and α is irrelevant. 3 Authors of Ref. [36] conjecture that a continuum limit may exist if bare coupling constants are fine-tuned in a specific way. of triangulations with a fixed lattice volume N 4 = N (4,1) + N (3,2) is, to leading order, proportional to exp(κ c 4 N 4 ), i.e. it grows approximately exponentially with N 4 . κ c 4 is a function of κ 0 and ∆. Due to this entropic effect, and thanks to the fact that in the bare CDT action (3) the κ 4 is conjugate to N 4 , the leading behaviour of the partition function (4) ) and the partition function is divergent for κ 4 < κ c 4 . In numerical simulations it is more practical to fix the total lattice volume N 4 and in such a case one should also choose κ 4 ≈ κ c 4 (N 4 , κ 0 , ∆). 4 The volume is usually controlled by introducing an additional volume fixing term to the bare action. In this work we use a quadratic volume fixing 5 which makes the total number of (4, 1) simplices oscillate aroundN (4,1) , and the impact of such volume fixing can easily be removed from the numerical data. Finally, before starting the numerical simulations, one should choose the fixed topology of spatial slices Σ and the length / boundary conditions for the proper time axis T . For practical reasons the CDT simulations are usually done for a periodic time axis of fixed length with t tot spatial slices, resulting in a global spacetime topology M = Σ×S 1 . Most of the previous results of CDT were obtained for the spatial topology of a threesphere Σ = S 3 , and in this article we focus on the spatial topology of a three-torus Summing up, in the numerical simulations described below we set: M = T 3 × S 1 , t tot = 40 or 4, = 0.00002,N (4,1) = 80000 or 160000 and we scan the parameter space spanned by κ 0 and ∆. For each data point we adjust κ 4 to the critical value κ c 4 (N (4,1) , κ 0 , ∆), see Fig. 2, so we are effectively left with a two-dimensional parameter space: (κ 0 , ∆). 4 By performing simulations this way one in fact investigates the properties of Z(κ 0 , ∆, N 4 ) which is related to Z(κ 0 , ∆, κ 4 ) by a Laplace transformation. One can investigate the infinite volume limit as well as reconstruct Z(κ 0 , ∆, κ 4 ) from extrapolations of measurements done for different N 4 .
5 Note that technically we fix N (4,1) instead of N 4 , but for fixed values of κ 0 and ∆ the ratio N 4 /N (4,1) is approximately constant and volume fixing method does not have much impact on the results.

State of the art
The key assumption of CDT is a choice of fixed spatial topology used in computer simulations. Most of the numerical studies performed in the past were for the specific choice of a three-sphere and time-periodic boundary conditions resulting in the global spacetime topology M = S 3 ×S 1 . This particular choice led to many interesting results, including discovery of four distinct phases of spacetime geometry, historically called A, B, C dS [40,28] and C b [41,42], see Fig. 3.  In phase A spacetime disintegrates into many causally uncorrelated baby universes with very short proper time extension, see Fig. 4 where the typical spatial volume profile n t ≡ N (4,1) (t) (i.e. the number of simplices with 4 vertices in lattice time t) is plotted. 6 The Hausdorff and spectral dimensions of triangulations inside this phase are approximately equal two. The phase is the CDT analogue of the branched polymer phase observed earlier in EDT. In phase B the spacetime geometry collapses into a single spatial slice, see Fig. 4, with (probably) infinite Hausdorff and spectral dimensions. The phase is the CDT analogue of the crumpled phase of EDT and it does not have a physical interpretation. A nontrivial result is the existence of the phase C dS , also called the de Sitter phase, where one observes the dynamical emergence of a large scale four-dimensional geometry [40,43] consistent with a semiclassical (Euclidean) de Sitter universe [44,45], see Fig. 4. At the same time the spectral dimension shows a non-trivial scale dependence and ranges from ≈ 2 for short scales to 4 at large scales (diffusion times) [46,47,48]. In the phase C dS the quantum fluctuations of the spatial volume are well described by the minisuperspace reduction of the Hilbert-Einstein action [45,49], and the contribution of such fluctuations vanishes in the infinite volume limit. This phase can be interpreted as the infrared limit of CDT, consistent with Einstein's GR. The recently discovered phase C b , also called the bifurcation phase [41], has a very nontrivial spacetime geometry. The volume profile n t resembles the one observed in the C dS phase, see Fig. 4, but spatial homogeneity is strongly broken by the appearance of compact spatial volume clusters concentrated around vertices with macroscopically large coordination numbers present in the every second spatial layer [42,50]. This phase is still being studied carefully and conclusive physical interpretation of its geometry has not yet been found.
The phases are separated by first order (A − C dS ) and second (or higher) order (B−C b ) [51,52] phase transition lines. The recently discovered C dS −C b phase transition is also second (or higher) order [53,54] and all the phases might meet in a common point. 7 This may in principle allow one to establish the perspective continuum limit by approaching a second order transition line or a multiple point from the physically interesting phase C dS , thereby defining a smooth interpolation between the low and high energy regimes of CDT [55,56].
All results described above were obtained for a spherical spatial topology Σ = S 3 , but recently we have been interested also in imposing toroidal spatial topology. One reason this is of interest is the background independence of CDT. No background geometry is imposed. However, in the case of spherical spatial topology we saw semiclassical four-dimensional (Euclidean) de Sitter-like configurations emerge, around which there were well defined quantum fluctuations. By changing the spatial topology to a toroidal topology Σ = T 3 one would expect different semiclassical solutions to dominate, so if the emergence of semiclassical geometry is a universal aspect of CDT one would expect to observe completely different geometries at least for choices of bare coupling constants where one obtained semiclassical configurations for spherical spatial topology. Further, we were inspired by the functional renormalization group research adapted to the ADM-formalism, where the authors of [57,58] started to investigate time foliated (Euclidean) spacetimes with topology T d × S 1 . Recently a similar study was performed for the topology S d × S 1 [59] leading to conclusions that the renormalization group flow linking IR and UV fixed points is essentially independent of the spatial topology chosen. Our results suggest that in CDT with a fixed spatial topology of a three-torus (Σ = T 3 ) there exists a semiclassical phase C, similar to phase C dS earlier observed for the topology of a three-sphere (Σ = S 3 ) [60,61]. However the dynamically generated (flat) background geometry of the new semiclassical phase C observed in toroidal topological conditions is completely different than the (four-sphere) geometry of phase C dS observed for spherical topological conditions. As mentioned above this is a quite non-trivial result and it gives strong support to the idea that there is a phase of CDT where semiclassical geometries emerge. We have in addition shown that the spatial volume fluctuations are still well described by a suitable minisuperspace reduction of the (toroidal) Einstein-Hilbert action in this phase, and one is even able to measure quantum corrections with much higher precision in the toroidal case.
In the current work we want to investigate the existence of other phases in CDT with toroidal topological conditions Σ = T 3 and to check if the phase diagram is similar to the case of spherical topology. As phase transition studies are very resource consuming, in this paper we focus on the phase structure itself, and the order of the phase transitions will be investigated in forthcoming articles.

Order parameters for the phase transitions
Before investigating the phase diagram of CDT one should define appropriate order parameters (OP s) that capture differences between generic triangulations observed in various phases of quantum gravity. The differences are usually caused by breaking of some global or local symmetries of the triangulations making the order parameter jump (first order transition) or rapidly change (second or higher order transition) from one phase to another. At the same time there is some freedom in choosing the order parameters and one should do it carefully to get clear signals of the transitions.
In this work we focus on four order parameters, selected in such a way that they capture the most important differences between the CDT phases. We use the experience gained in previous studies of CDT with spherical spatial topology, where various phases of geometry were observed (see Section 3) and similar order parameters were used in the study of the phase transitions. The order parameters discussed in this paper can be divided into two groups. The first group comprises order parameters which capture very global properties of CDT triangulations, such as which were previously used in the analysis of the A−C dS and B−C b transitions observed in the case of spatial topology S 3 [51,52]. The A − C dS transition was related to the time-translation symmetry breaking from (a symmetric) phase A to a (less symmetric) phase C dS where one could see a spacetime blob structure, i.e. some macroscopic extension of the universe in time direction. This time-extended universe persisted in phase C b and the symmetry was further broken in (the least symmetric) phase B, where the universe collapsed to a single spatial layer, see Fig. 4. Such symmetry differences resulted in the order parameter OP 1 being large in phase A, medium in phases C dS and C b , and small in phase B, see Fig. 5. At the same time the breaking of causal connections of neighbouring spatial layers caused the order parameter OP 2 to be small in (time uncorrelated) phases A and B and large in (time correlated) phases C dS and C b , see Fig. 5.
The second group of order parameters focuses on microscopic properties of the underlying CDT triangulations, namely the shape of the spatial volume profiles n t : and the existence of vertices of very large order: where O(v) is the vertex coordination number, i.e. the number of simplices sharing a given vertex v, thus OP 4 is just the coordination number of the highest order vertex in a triangulation. When one looks at the volume profiles n t , see Fig. 4, one observes that the profile is quite smooth in phase C dS where there are no big differences in n t and n t+1 . In phase C b the volume profile narrows in time direction causing slightly bigger differences between n t and n t+1 . The profile is much less smooth in phase A, where a kind of "zig-zag" shape is observed and very non-smooth is phase B where there is a sudden jump in the slice where all spatial volume is concentrated. As a result, the OP 3 is small is phase C dS , medium in phases A and C b , and large in phase B, see Fig. 5. Last but not least, the existence of high order vertices is related to the formation of spatial volume clusters which are characteristic for the bifurcation phase C b [42,50] and phase B (one huge volume cluster in the collapsed slice), resulting in large OP 4 inside these phases. In phases C dS and A some spatial volume concentrations may also form from quantum fluctuations but they do not form any distinguished large-scale structures, making the OP 3 small, see Fig. 5. The order parameters are summarised in Table 1. The right plot shows the OPs as a function of κ 0 for fixed ∆ = 0.6, which corresponds to a horizontal line in the phase diagram in Fig. 3 starting in phase C b (κ 0 <∼ 1.0), going through phase C dS (∼ 1.0 < κ 0 < 4.7), to phase A (κ 0 > 4.7). The OP 1 , OP 2 and OP 4 were rescaled to fit into a single plot, and due to a very large range of OP 3 (a few orders of magnitude) we plot the rescaled ln OP 3 instead of OP 3 , the rescaling in the left plot being identical to that in the right plot. We expect that in the case of toroidal spatial topology Σ = T 3 there exist phases similar to those observed for spherical spatial topology and that we should see a similar behaviour of the order parameters defined above. Thus one can scan the (κ 0 , ∆) parameter space and measure the averages OP 1 , ..., OP 4 in order to identify the various phases, if they exist. In order to establish a precise position of the phase transition in the parameter space one can also look at the susceptibilty of an order parameter which should peak at the phase transition point. For each case we will choose the set of the order parameters which gives the clearest signal / noise ratio. We will also analyze or Note that there is always some freedom in choosing the order parameter, and instead of OP one can use any monotonic function f (OP ) of OP . For small fluctuations of OP around the mean value OP one has Up to a numerical constant, the χ * OP and χ * * OP are obtained for f (OP ) = √ OP and f (OP ) = ln OP , respectively, and such a choice is useful when an order parameter changes by a few orders of magnitude at a phase transition. The approximation (13) seems to work very well, see Fig. 6 where the results obtained by a redefinition OP → f (OP ) in the raw data and by using approximation (13) cannot be optically distinguished.

Phase structure in toroidal spatial topology
We begin with a rough scan of the CDT parameter space for fixed spatial topology of a three-torus and time-periodic boundary conditions, i.e. spacetime topology M = T 3 × S 1 . We choose the lattice volumeN (4,1) = 80000 and the length of the time period t tot = 40. For such a choice of simulation parameters one could observe all four CDT phases in the the case of spherical topology. We split the parameter space 8 (κ 0 , ∆) into a grid of equally separated points (see Fig. 7) and run simulations in which we measure spatial volume profiles n t (see Fig. 8) and the order parameters OP 1 , ..., OP 4 described in Section 4 (see Fig. 9). 8 The κ 4 coupling constant is adjusted to the critical value κ c 4 (N (4,1) , κ 0 , ∆), see Section 2 for details.   Fig. 7, and the right plot shows the OPs as a function of κ 0 for fixed ∆ = 0.6, which corresponds to a horizontal line in the phase diagram in Fig. 7. The OP 1 , OP 2 and OP 4 were rescaled to fit into a single plot, and due to a very large range of OP 3 (a few orders of magnitude) we plot the rescaled ln OP 3 instead of OP 3 , the rescaling in the left and right plots being the same. The qualitative behaviour of the OPs is the same as for phases A, B and C dS in CDT with spherical spatial topology (see Fig. 5 and Table 1). Approximate positions of the B − C and the C − A phase transitions are marked by dashed lines on the left and right plots, respectively.
By looking at these data one can easily notice three distinct phases of spacetime geometry, denoted A, B and C. The position of phases A and B on the phase diagram in Fig. 7, as well as the spatial volume profiles (series of uncorrelated spatial slices in phase A, and time-collapsed volume structure in phase B) and values of the order parameters (see Fig. 9 and Table 1) suggest that the phases are in one-to-one correspondence to phases A and B of CDT with spatial topology Σ = S 3 (see Fig. 3).
Existence of phase C in spatial topology Σ = T 3 was already reported in [60,61], where the authors noticed that the spatial volume profile n t is highly correlated and can be well described by a toroidal minisuperspace-like action with small quantum corrections. It can be argued that this is the toroidal analogue of the semi-classical phase C dS , earlier observed for the spatial topology Σ = S 3 (see Fig. 3). This is further confirmed by behaviour of the order parameters OP 1 , ..., OP 4 (see Fig. 9 and Table 1).
The fact that we don't see the toroidal analogue of the bifurcation phase C b in the above data is not very surprising as the (expected) constant volume profile of spatial slices causes the volume of a single slice to be too small (n t ∼N (4,1) /t tot = 2000) to allow for a creation of large volume clusters and consequently for the emergence of highorder vertices. The same effect appeared in the spherical topology for small systems where high order vertices could be observed only in slices with spatial volume higher than the, so-called bifurcation point volume. In CDT with spherical spatial topology the volume profile n t had a characteristic blob structure (see Fig. 4) with the central part volume much higher than the average volumeN (4,1) /t tot and thus the choice of N (4,1) = 80000 and t tot = 40 was good enough to let the central part volume be higher than the bifurcation point volume. As a result the high-order vertices could form inside the blob part of the spherical CDT triangulations, which seems not to be the case in the toroidal CDT with flat volume profiles.
To circumvent this obstacles, i.e. to let the average spatial volume of the toroidal CDT triangulations exceed the bifurcation point volume, we decided to pursue a detailed study of the toroidal CDT phase structure for much bigger average volumē N (4,1) /t tot = 40000, by settingN (4,1) = 160000 and t tot = 4, respectively. To use our computer resources effectively, we have also carefully fine-tuned the grid of the (κ 0 , ∆) points in which we performed numerical simulations, such that we have much higher precision near expected phase transition points, see Fig. 10 where the phase diagram is shown. These results reconfirm the findings of the initial phase transition and t tot = 4. Points denote the actual numerical simulations. The measured phase transitions are marked by solid lines with error bars denoted by shaded areas of the same color. The error bars for the C − C b transition are due to the observed hysteresis of the measured order parameters (see Section 6 for discussion), the error bars for the other phase transitions are purely due to the resolution of the measurement grid. Dashed black lines are extrapolations of the measured phase transition lines. All four distinct phases of geometry can be identified. Phase A is observed for large values of κ 0 (orange points), phase B for small values of ∆ (black points), phase C b for small κ 0 and medium ∆ (blue points), and phase C for small κ 0 and large ∆ (green points). study from Fig. 7, and also confirm the existence of the fourth phase, denoted C b (blue points in Fig. 10). The phase is again a toroidal analogue of the bifurcation phase C b observed in the spherical topology, with spatial homogeneity broken by a formation of volume clusters around high order vertices emerging in the every second spatial layer, see Fig. 11.The behaviour of all four order parameters as functions of κ 0 and ∆ is shown in Fig. 12 and their susceptibilities in Fig. 13. These results can be used to draw the phase transition lines with high precision, as presented in Fig. 10. The phase structure looks very similar to the one observed for a spherical spatial topology, see Fig. 3.  Fig. 10, and the right plot shows the mean OPs as a function of κ 0 for fixed ∆ = 0.4, which corresponds to a horizontal line in the phase diagram in Fig. 10. The OP 1 , OP 2 and OP 4 were rescaled to fit into a single plot, and due to a very large range of OP 3 (a few orders of magnitude) we plot the rescaled ln OP 3 instead of OP 3 , the rescaling in the left plot being identical to that in the right plot. The qualitative behaviour of the OPs is the same as for phases A, B, C dS and C b in CDT with spherical spatial topology (see Fig. 5 and Table 1 We studied the phase diagram of four-dimensional CDT with toroidal spatial topology (Σ = T 3 ) and time periodic boundary conditions, see Fig. 14. Our results confirm the existence of four distinct phases of quantum geometry which are direct analogues of the phases previously observed for the spherical spatial topology (Σ = S 3 ). The position of the critical lines in the (κ 0 , ∆) parameter space is also very similar in both topologies. The A − C and B − C b phase transition lines measured in the toroidal case are slightly shifted compared to the lines measured in the spherical case. Using our present data one cannot completely exclude that the shifts are real effects which might be attributed to the topology change. However it is much more likely that the shifts are due to final size effects as the positions of (pseudo-)critical points in the parameter space depend on the lattice volume and similar parallel shifts were observed in the spherical case when the lattice volume was increased [51,52,53,54]. On the one hand, the finite size effects should be much stronger in the toroidal CDT, where the minimal triangulation is much larger than in the case of the spherical CDT [60]. On the other hand, the length of the (periodic) time axis used in the toroidal CDT simulations presented herein was much shorter than in the spherical simulations, and the resulting effective volume per slice was much larger in the toroidal case.
The critical region near the point where the phase transition lines meet is the most interesting place to concentrate on, since, following the asymptotic safety arguments, it is a natural candidate for the physical UV limit of CDT [55,56]. This region could not be studied with a sufficiently high precision for the spherical topology. The reason was purely technical: the local Monte Carlo algorithm we use to update triangulations was very inefficient in this critical region. As a result we couldn't make precise measurements there and our conjecture about the existence of a common "quadruple point" where all four phases meet was based on an extrapolation of the measured phase transition lines. The situation is different for the toroidal case, where the Monte Carlo algorithm works fine in the critical "corner" region of the parameter space, as can be seen from a plot shown in Fig. 14. In this case the common "quadruple point" seems less likely. The C b − C phase transition line is now shifted slightly to the left and tilted in the (κ 0 , ∆) plane compared to the C b − C ds phase transition line observed in the case of spherical spatial topology. Consequently, in the toroidal case, it seems more likely that we have the "old" triple point around (κ 0 = 4.50, ∆ = −0.05) where A, B and C phases meet and a "new" triple point around (κ 0 = 3.75, ∆ = −0.02) where the bifurcation phase C b meets the phases B and C. As a result there exists a region in the parameter space, where one can observe a direct transition line between the semiclassical phase C and the collapsed phase B, see Fig. 14. We will concentrate our future precise measurements in the "critical corner" region to determine accurately the phase diagram. Although the grid of measured points presented in the article seems rather dense, we cannot exclude that the exact shape of transition lines is more complicated in the sense that they may bend and meet in one "quadruple point" (again see Fig. 14). In this precise study we will also measure a sequence of volumes which should enable us to analyze finite size effects affecting the infinite volume position of the phase transition lines more accurately. The results will be published in forthcoming publications.
In the data presented above we did not measure the order of the phase transitions. Such a study would require massive numerical simulations to be performed for a suitable choice of various lattice sizes enabling one to extrapolate the results to the infinite volume limit and to measure critical exponents. Nevertheless we have made some initial observations. First, in CDT with toroidal spatial topology one can observe a clear C −A transition signal (peak in susceptibilities) both for the OP 1 and OP 2 parameters, as it was the case for the spherical CDT. Nevertheless the behaviour of the parameters at the phase transition is quite different from the behavior originally observed for CDT with spherical spatial topology where one could observe that the order parameters jump between two clearly separated metastable states, one for phase A and one for phase C dS [52]. This suggested it was the first-order transition and this was confirmed by a detailed finite size analysis. In the toroidal case the order parameters change smoothly between the two phases and one does not observe any separation of states. While this could be an indication of a higher order transition, it is more likely that it reflects that one is using a constraint of the four-volume (namely (5) where N (4,1) is kept fixed) different than that used in [52] (where N 4 was kept fixed). A similar phenomena were observed for the B − C b transition [50]. Clearly the only way to settle the issue is to perform a carefull finite size analysis to determine the order of the transition.
Secondly, the C b − C dS phase transition in case of CDT with spherical spatial topology was found to be a second (or higher) order transition [53,54]. Now, in the toroidal case, one can observe a clear hysteresis of all measured order parameters when moving from phase C to phase C b or the opposite. There is a large region in the parameter space, denoted by the orange shaded area between phases C and C b in Fig. 14, where the value of an order parameter depends on its initial value, or more precisely on the geometric configuration (triangulation) used to initiate Monte Carlo simulations. If one starts with an initial triangulation from phase C and makes simulations inside the hysteresis region the generated triangulations persist in phase C, if one starts instead with a triangulation from phase C b the system persists in phase C b , even for very long simulation runs (a few months of CPU time, or a few ×10 12 attempted Monte Carlo moves). Of course if one goes deep enough into phase C or alternatively into phase C b (outside the shaded region in Fig. 10) the system finally thermalizes to phase C or phase C b , respectively, independent on the starting configuration. This is illustrated in Fig. 15, where we show the OPs as a function of ∆ for fixed κ 0 = 2.0. For each ∆ we plot the measurements done in many independent Monte Carlo runs with different initial configurations. The hysteresis region between phases C and C b is clearly visible for ∆ ≥ 0.38. Again, while this can indicate a first order transition, and thus a change of transition order with topology, it could also be an algorithmic issue with the Monte Carlo simulations since our configurations are relatively small and the toroidal topology clearly is more constraining than the spherical topology. Again the only way to settle the issue is to perform a proper finite size analysis, which is of course quite computer demanding. Summing up, the general phase structure of CDT with toroidal spatial topology is very similar to phase structure of the CDT with spherical spatial topology and even the positions of phase transition lines are almost the same. We have observed some qualitative difference in phase transitions but it is impossible to say presently if this implies that the order of some of the phase transitions should change with topology. It would indeed be somewhat surprising if the transitions can be related to continuum physics, in particular UV physics which one would imagine was related to short distance phenomena. Short distance phenomena should be insensitive to topology. However, the precise nature of the transitions requires further studies which will be presented in forthcoming publications.