The phase diagram of the multi-matrix model with ABAB-interaction from functional renormalization

At criticality, discrete quantum-gravity models are expected to give rise to continuum spacetime. Recent progress has established the functional renormalization group method in the context of such models as a practical tool to study their critical properties and to chart their phase diagrams. Here, we apply these techniques to the multi-matrix model with $ABAB$-interaction potentially relevant for Lorentzian quantum gravity in 3 dimensions. We characterize the fixed-point structure and phase diagram of this model, paving the way for functional RG studies of more general multi-matrix or tensor models encoding causality and subjecting the technique to another strong test of its performance in discrete quantum gravity by comparing to known results.


Introduction
Understanding the quantum properties of spacetime is a fascinating challenge. A variety of approaches is currently being explored concertedly; more recently with an increased interest in understanding relations between different perspectives. The matrix/tensor model approach [1][2][3][4][5] is located at a particular vantage point within this "landscape" of theories, with potential links to a number of different ones: Firstly, its origin in two-dimensional gravity is closely linked to string theory. Secondly, its generalization to higher dimensions is closely connected to a model that is being explored in the context of the AdS/CFT conjecture [6,7]. Thirdly, in matrix models a tentative connection to an asymptotically safe fixed point in the vicinity of two dimensions has been found [8] and conjectured in higher dimensions [9]. Fourth, this class of models provides a combinatorial approach to dynamical triangulations, complementing computer simulations of the latter [10][11][12]. This rich set of connections motivates a concerted study of matrix and tensor models and in particular the search for a universal continuum limit in these. Functional Renormalization Group techniques have been developed for these models [8,13,14] and are being applied in [9,[15][16][17] and related tensorial group field theories [18][19][20][21][22][23][24][25][26][27][28][29][30], with the potential to discover a universal continuum limit beyond perturbation theory, see also [31,32] for related developments with the Polchinski equation. To set the stage for our studies, we first provide a more in-depth overview of the relevant quantum-gravity approaches.
The above approaches center on the path integral where the integration runs over all field histories, given by spacetime metrics g of the d-dimensional manifold M up to diffeomorphisms thereof at fixed spacetime topology. Working with the Einstein-Hilbert action, one can make sense of the weak-field limit in an effective-field theory framework [33], but encounters perturbative non-renormalizability, resulting in a loss of predictivity, beyond [34,35]. A possible way out of this issue is to replace the Einstein-Hilbert action by one which allows for a unitary and perturbatively renormalizable QFT. However, this might dispense with micro-causality as in the case of higher-derivative gravity [36][37][38] or with Lorentz invariance as in Hořava-Lifshitz gravity [39], see [40] for a recent review. An alternative pathway to quantize gravity within the continuum formulation of the path integral is explored by the asymptotic-safety program. To sidestep the problems of the perturbative quantization, this approach is based on an interacting fixed point in the Renormalization Group (RG) flow for gravity in the UV [41][42][43]. If it exists, such a fixed point provides a well-defined continuum limit to the path integral. At the same time, it generalizes perturbative renormalizability by ensuring that the low-energy limit is parameterized by only a finite set of free parameters, namely the relevant directions of the fixed point. There are several techniques suitable to explore asymptotic safety in gravity, falling into the two broad categories of lattice approaches and continuum approaches. A much used method pioneered for gravity by Reuter [43] is provided by the functional Renormalization Group (FRG) [44,45], see [46] for a review. At its core lies the implementation of the Wilsonian idea of a coarse graining operation which progressively eliminates short scale fluctuations. Indeed, all explicit calculations within truncated RG flows find evidence for the existence of such a fixed point, defining the Reuter universality class, providing compelling indications for asymptotic safety in Euclidean gravity, see, e.g., [47][48][49][50][51][52][53] for recent reviews and introductions. Rephrased in a lattice-like language, such a universality class enables one to take a universal continuum limit. Open questions in this approach have been discussed in [52,54] and include the fate of background independence, given the assumption of an auxiliary background metric therein [43,55,56]. Moreover, since the signature of the setup is Euclidean and one can in general not expect the Wick rotation to exist in a quantum gravitational context [57][58][59][60], it is open how to relate these results and in particular the feature of asymptotic safety to Lorentzian quantum gravity; for a first step in this direction see [61]. Similarly to the characterization of other interacting fixed points, a concerted use of several different techniques could in the future provide a qualitatively and quantitatively robust grasp of the fixed point and its properties. In the case of quantum gravity, background-independent, Lorentzian approaches to quantum gravity are of particular interest to explore as techniques that can complement the FRG results.
One such promising approach to evaluate the path integral over geometries, possibly extended by a sum over topologies, and to discover a universal continuum limit, a.k.a. asymptotic safety, is by means of a sum over discrete triangulations, together with exchanging the continuum action with its discretized reformulation. In spite of the physical discreteness that some of these settings exhibit, such as, e.g., Loop Quantum Gravity [62][63][64], taking a universal continuum limit is a central goal also in this setting. For instance, in group field theory [65][66][67][68], promising results regarding asymptotic freedom [69,70] and asymptotic safety [22,71] have been obtained, see [67] for a review. Similarly, RG techniques are being developed and applied in the search for a critical point in spin foams, see [72,73]. In this setting, studies of the phase diagram of quantum gravity have only recently started [74][75][76][77]. In contrast, in the Euclidean and Causal Dynamical Triangulation approaches (EDT, CDT) [11,12,78,79], much is already known about the phase diagram. Going beyond two dimensions, early numerical studies using Monte-Carlo methods [80][81][82][83] have only recovered unphysical geometries [84][85][86] in the Euclidean setting, though the inclusion of additional terms in the action or measure of the path integral, associated with additional tunable parameters, might change the situation [87,88]. It is argued that these pathological configurations are the result of topology change leading to spaces called branched polymers which are built from one-dimensional branched-out filaments [82,84,89,90]. At the classical level, as long as Morse geometries are excluded [91][92][93], it is long known that topology change leads to a degenerate local light cone structure and thus to a violation of micro-causality [94]. Thus, in the CDT approach discrete spacetime configurations with spatial topology change are excluded [95]. This leads to a much better behaved theory with the potential to produce a phase with physically relevant, i.e., extended macroscopic geometries [11,12,79], bordered by a second-order phase transition [96][97][98][99] that enables a continuum limit. Yet, a key challenge in the numerical approach to dynamical triangulations remains to follow RG trajectories towards the continuum limit [100] and calculate the scaling spectra to determine and characterize the universality class.
Matrix [1] and tensor models [2][3][4][5] sit at a confluence of several of these approaches. They encode random discrete geometries in dimensions d ≥ 2. Their general idea is to represent d − 1 simplices corresponding to building blocks of geometry as rank d-tensors. The tensor action encodes how to glue these building blocks together to construct d-dimensional discrete geometries. These correspond to the Feynman diagrams in the perturbative expansion of the tensor path integral. This establishes a duality between tensor models and the discrete gravitational path integral. Indeed, in their simplest form matrix and tensor models can be understood as generators of EDTs. However, apart from the case in d = 2 [1], their continuum limit so far only leads to geometrically degenerate configurations [2,3], the same way as in EDT, or planar ones [101]. Recent results indicate the possibility for non-trivial universality classes [9,15]; however, the geometric properties of the corresponding phases have not been investigated yet. If a universal continuum limit within a phase with desirable geometric properties can be taken, asymptotic safety can be confirmed in a background-independent fashion and with straightforward access to the scaling exponents. Thus, the critical role of universality in these models has been emphasized in [16,102]. Yet, in studies beyond two dimensions, the inclusion of causality has remained an open question -similarly as in the continuum approach to asymptotic safety. The success of CDT may be taken as a hint to consider additional structure enforcing micro-causality to recover higher-dimensional physical continuum geometries from such models. This idea has already been implemented in the context of matrix models for 2d quantum gravity, giving rise to a description equivalent to CDT in (1 + 1) dimensions [103]. In order to explore the impact of Lorentzian signature in higher dimensions, one could naturally try to impose such causality conditions on a model for tensors of rank d ≥ 3. Yet, there already exists a proposed correspondence between a multi-matrix model with CDT in 2+1 dimensions [104]. More precisely, it corresponds to a Hermitian two-matrix model with ABAB-interaction which also has intimate connections to vertex-models of statistical physics [105][106][107][108]. It leads to a variant of CDT defined on an enlarged ensemble of configurations which also allows for specific degeneracies of the local geometry, dubbed "wormholes" in the literature [104]. The purpose of this article is to chart the phase structure of this matrix model and in particular to study its continuum limits by means of the FRG methodology. The application of the FRG in the discrete quantum gravity context facilitates a backgroundindependent form of coarse graining where the number of degrees of freedom serves as a scale for a Renormalization Group flow. This program was started by analyzing matrix models for 2d Euclidean quantum gravity in [8,13] and has by now been extended to tensor models for Euclidean quantum gravity in 3 and 4 dimensions [9,14,15,109], see [16] for a review. Related developments for non-commutative geometry [110,111] and tensorial group field theories [18][19][20][21][22][23][24][25][26][27][28][29][30] exist. Recently, a first FRG analysis of the above-mentioned causal matrix model for CDT in 1 + 1 dimensions has been completed [112].
The article is organized as follows: In Sec. 2 we discuss causality in the matrixmodel context and review the relation of the ABAB-matrix model to CDTs, following Refs. [104,113,114]. In Sec. 3, we briefly review functional RG techniques and apply them to the ABAB matrix model. We then present results for the phase diagram and fixed-point structure. In Sec. 4 we review what is known about the phase diagram in the literature [104-106, 113, 114] and compare our results to it. Finally, in Sec. 5 we discuss implications of our results and future directions.

Causality and matrix models
Spacetime is rather distinct from space, both at the conceptual as well as mathematical level. Therefore, it is crucial to take this difference into account in quantum gravity, with the distinct phase diagrams of CDT and EDT constituting a clear example of the impact of causality. In the matrix and tensor model approach, the additional structure imposed on discrete configurations by causality can be implemented through additional degrees of freedom: In Refs. [103,112], this is done by an external matrix, whereas the ABAB model uses two dynamical matrices to generate configurations which carry imprints of causality. This strongly motivates the further development of FRG techniques for multi-matrix/tensor settings, such that similar developments can be made possible in higher dimensions. In this section, we will review the relation of the ABAB matrix model to 2+1 dimensional discrete spacetime configurations. In particular, we will follow Refs. [104,113,114] to also review the connection to CDT. The ABAB matrix model is defined by the following partition function  where A and B are Hermitian N × N matrices and J A and J B are the respective external (N × N -matrix) sources. Its Feynman diagrams are ribbon diagrams with two distinct types of lines; such that the duals of the three types of vertices correspond to three distinct squares, cf. Fig. 1. This already highlights that the model reduces to the standard twodimensional gravity case, whenβ → 0. In the presence ofβ, 2+1 dimensional structure can be encoded [104].
In CDTs, one only considers geometries with Lorentzian signature which admit a global foliation in proper time t and disallow topology change such that micro-causality is rigidly maintained. The basic building blocks are pyramids and tetrahedra with distinct timelike and spacelike edges, of which three distinct types make up the discrete configuration, cf. Fig. 2, see also Refs. [104,[115][116][117]. The correspondence with the ribbon graphs of the ABAB matrix model arises when considering the spatial hypersurface at t + a/2, also indicated in Fig. 2, where the three distinct types of squares dual to the vertices of the matrix model, cf. Fig. 1, appear. In the t + a/2 spatial planes, quadrangulations are formed as, e.g., in Fig. 3, the duals of which are the ribbon graphs of the matrix model. This already suggests that any CDT configuration can be encoded in terms of a Feynman diagram of the matrix model, which would motivate setting the CDT partition function equal to the free energy of the matrix model, as usual in the correspondence between triangulations and matrix models. To see this in more detail and further discuss whether or not there is an exact correspondence, let us follow the discussion in [104].
Indeed, the entire information on the CDT partition function can be encoded in the onestep propagator which provides the transition amplitude between the spatial hypersurfaces at t and t + a. It is this property that enables a connection to the ABAB matrix model. Indeed, it has been shown in [104] that the Euclideanized transition amplitude between the triangulation ∆ of the spatial hypersurface Σ t and ∆(Σ t+a ) is given by (cf. Eq. 6 in [104]) whereT is the transfer matrix. In the above expression, N 41 is the number of squares at t and there are N 14 squares in t + a. The sum is over all intermediate quadrangulations, and Figure 2: The three fundamental building blocks of CDT in 2 + 1 dimensions which interpolate between two consecutive spatial hypersurfaces at integer times t and t + a, cf. Fig. 1 in Ref. [104]. The numbers at the pyramids and the tetrahedron refer to the number of vertices in the quadrangulations at constant integer time t and t + a. In between, their intersections with the t + a 2 -plane are shown, giving rise to an equilateral quadrangulation thereof in terms of blue, red and bi-colored squares. The colorization emphasizes that 2 + 1 dimensional information is encoded in a 2d setting with colors. N CDT (N 41 , N 14 , N 22 ) denotes the total number of quadrangulations at fixed N 22 . We note that this expression holds for spaces of spherical topology triangulated by a large number of squares. Further, the above expression assumes the discretized Einstein-Hilbert action, such that λ and κ are related to the bare cosmological constant Λ and bare Newton coupling G N as follows and a is the lattice spacing. The free energy of the matrix model N 2 F = − ln Z[0, 0] is given by With the identifications this makes the close relation to the CDT-partition function obvious. However, the number of configurations generated by the CDT model is smaller than that generated by the matrix model, i.e., N MM > N CDT : The matrix model generically generates disconnected subgraphs via so-called touching interactions (see below). This is most easily seen by switching off the interactions in one sector,ᾱ 2 → 0 and integrating out B. The resulting single-matrix model contains multi-trace terms. In the dual triangulation picture, these lead to branched trees of spherical bubbles, see Fig. 4. Clearly, in the generalized situation whereᾱ 2 is reinstated, touching interactions will also be present and yield such undesirable quadrangulations. Such pathologies of the local geometry are disallowed by construction in CDT [11]. Despite these differences, the similarities between the ABAB matrix model and CDTs reinforce the more general notion that causality can be imposed on matrix and tensor models by enlarging the field content of the model and introducing a second matrix/tensor that ultimately enables a distinction of spacelike and timelike edges in the dual triangulation. This motivates us to perform a functional RG analysis of the ABAB matrix model. On the one hand, a comparison to existing results in the literature will enable us to conduct a novel, powerful test of the performance of this technique. On the other hand, this will pave the way for future studies of multi-field models that encode causality in the interaction structure.

FRG method for matrix models
The FRG is a powerful and versatile tool to implement the Wilsonian renormalization program. In a standard, local field-theoretic setting, given the Euclidean path integral, one introduces a regulator function which suppresses the functional integration of modes below a given momentum cutoff k which correspond to the slow modes in the Wilsonian perspective. By progressively lowering the values of k, one carries out the complete integration over all modes. Hence, instead of performing the path integral at once, it is computed in a momentum-shell-wise fashion. The central object in the FRG is the so-called flowing action Γ k which interpolates between the classical action S (when k → ∞) and the full effective action Γ (when k → 0). It satisfies a flow equation [44,45] which has a simple one-loop structure, making it very efficient for practical calculations, where approximations have to be employed. As the full propagator enters, nonperturbative physics is captured, despite the one-loop structure. At the formal level, the equation is exact and no approximation enters its derivation, rendering it formally equivalent to the path integral. For a recent review of the FRG in a broad range of contexts from condensed matter to quantum gravity, see [46]. The Wilsonian picture as described above relies on a background structure which provides a notion of momentum scales. Such a momentum scale is used as a coarse-graining parameter. In the case of matrix models for quantum gravity, there is no such background. Indeed, these models can be thought of as pre-geometric, with a smooth spacetime and notions of distance emerging in the continuum limit. Yet, as introduced in Ref. [118], a notion of renormalization group can naturally be defined if the number of entries of the matrices is taken as the coarse-graining parameter. Hence, integrating out "fast modes" corresponds to integrating out the outermost rows and columns. In Ref. [13], this idea was translated into an exact flow equation for matrix models, paving the way for similar developments in group field theories [18] and tensor models [14]. More recent developments can be found in Refs. [9,15,17,111], see also [16]  We demand that the regulator ∆S N has the structure of with R N being independent of the random matrix M . It is required to satisfy the following three properties: which suppress the matrix entries in the block a, b = 1, ..., N and facilitate that the "UV" matrix entries with a, b > N are integrated out [13]. It is then a formal manipulation to show that the flowing action 1 Γ N [M ] defined by Indices were suppressed for simplicity. The Tr represents a sum over the indices and ∂ t ≡ N ∂ N . Such a derivative should actually be a finite difference at finite N . However, we are interested in the large-N limit, justifying the explicit use of the derivative. The flowing action Γ N contains all terms compatible with the symmetries of the model. Thus, it can be expanded as with O[M ] denoting all the operators satisfying a given symmetry andḡ I standing for "dimensionful" couplings. By expanding the right-hand side of Eq. (3.8) in terms of the same basis O I , it is possible to project out the flow of each couplingḡ I , i.e., to compute the beta functions of the theory. In standard local field-theoretic settings, one is interested in the dimensionless version of the couplings and their flow as these contain the information on (quantum) scale invariant points. In the present case, there is no local notion of scale; nevertheless "dimensionless" couplings can be defined that absorb factors of N . Specificially, the couplingḡ I is related to its dimensionless counterpart g I by the relationḡ I = N [d I ] g I , where [d I ] denotes the canonical scaling dimension ofḡ I . Since the renormalization group parameter N is dimensionless, the assignment of canonical dimensions to couplings does not follow from simple dimensional analysis. Actually, as discussed in [8,14,15], the scaling with N is fixed by demanding that, at large-N , the system of beta functions is non-trivial and autonomous, i.e., the flow does not depend explicitly on N . One can follow similar arguments in the framework of Dyson-Schwinger equations, see, e.g., [119] as well as the Polchinski equation [31]. The beta functions of the dimensionless couplings are thus The function ∂ tḡI is read off from the flow equation (3.8) by a suitable projection onto the basis defined by Eq. (3.9).
In this work, we are interested in multi-matrix models. As in standard field theory with multiple fields, the flow equation can easily be derived in this setting. For concreteness, we consider a model with two interacting Hermitian matrices A and B (but the argument extends to a generic set of matrices) with the following partition function, with J A and J B being the sources associated, respectively, to A and B and S[A, B] is the classical action of the model. The regulator ∆S N is introduced as where Φ I ≡ {A, B} and R ab,cd N,IJ (a, b) satisfies the properties (3.4), (3.5) and (3.6) transferred to the multi-matrix case. Thence, the flow equation is derived in full analogy to the singlematrix model, leading to where Tr represents a sum over matrix indices. Equation (3.13) easily extends to a generic number of matrices. In the next subsections, we will explicitly apply it to the case of the two-matrix model related to CDT in 2 + 1 dimensions, the details of which we specified above.
As it is crucial for the developments that follow, we emphasize an important property of the flow equation: While the path integral Eq. (3.11) requires the specification of a classical action S, the flow equation is independent of the classical action. Instead, it provides a local vector field in the space of all couplings, indicating how the dynamics changes under a finite RG step. A classical action can be specified as an initial condition to integrate this flow and obtain the effective action. On the other hand, the flow equation can also be used to search for special points in the space of couplings, which correspond to fixed points of the RG flow. Such fixed points then give rise to a particular proposal for a classical (or microscopic) action as the starting point of the flow.

General setup of the flow equations
To obtain the set of beta functions of the couplings, one has to project the flowing action onto the couplings of the corresponding operators. In practice, this is achieved by means of a series expansion of the rhs of Eq. (3.13) in terms of powers of A and B, the P −1 F expansion. To this end, one rewrites the denominator of the flow equation as with the fluctuation matrix F N and inverse propagator defined by The right-hand side of the flow equation is then expanded as wherein STr denotes a sum over the available matrices and the indices thereof. In the next step, we split the Hermitian matrices A and B into symmetric and anti-symmetric parts With this parameterization, we compute the variations for the Hessian. After doing that, the next step typically is to choose a particular field configuration which facilitates the projection of the right-hand side of the flow equation onto the corresponding beta function of interest.
For the purposes of this work, it suffices to take A 2ab = B 2ab = 0. The Hessian then has the following structure, In the very last step, after having computed the beta functions which also map combinatorial differences between operators of the same power in A and B, we project the remaining field configurations onto A 1ab = a 1 δ ab and B 1ab = b 1 δ ab , see Sec. 3.3. The regulator is chosen in a diagonal form 20) and Therein, Z A , Z B are wave-function renormalization factors and the details of the cut-off procedure are captured by the function so that the regulator is closely modelled after Litim's cutoff [120]. With these expressions, the computation of the inverse propagator P N , Eq. (3.15), is straightforward. A relevant detail repeatedly used in the P −1 F expansion for concrete calculations is then given by the approximation valid at leading order in 1/N . In this expression, we introduce the anomalous dimensions as η A/B = −∂ t ln Z A/B . Finally, as discussed in Refs. [13,16] and Sec. 3.1, the couplings in matrix and tensor models have an inherent dimensionality in spite of having no natural notion of momentumscale. It dictates their behavior with respect to rescalings in N . In our model, one has the following rescalings for the dimensionful couplings where α 1/2 and β are "dimensionless" couplings. This assignment of scaling dimensionality is chosen in such a way as to facilitate a 1/N -expansion of the beta functions where the leading coefficient is O(N 0 ). Because of this choice a sensible continuum limit exists. Equipped with this we are ready to calculate the beta functions from the flow equation. (3.26) Compared to the analysis of the single-matrix model [8,13], we stick to the sign conventions of Refs. [105,106] in which the interaction terms carry a negative sign.

Beta functions and fixed points for the ABAB matrix model
From Eq. (3.26), one can extract the fluctuation matrix F N which is the remaining ingredient we need in order to compute the P −1 F expansion. Its non-vanishing entries are given by: In the following, we present the computations relevant to obtain the set of beta functions from the flow equation. First, we briefly summarize our line of action: Given the simple truncation, we only need to compute the P −1 F expansion up to order 2. At 0th order one obtains a field-independent constant which can be absorbed in the normalization of the path integral and is thus of no interest hereafter. From the 1st-order-contributions, we compute the expressions for the anomalous dimensions η A and η B . Then, at 2nd order we yield the beta functions for the couplings.
Starting off with the 1st order of the expansion, we have (3.33) In the next step we focus only on single-trace contributions and thus discard any occurring double-trace terms (which have the structure Tr(. . .)Tr(. . .)). Then, as discussed in Sec. 3.2, the second part of the projection is applied, namely inserting A 1ab = a 1 δ ab and B 1ab = b 1 δ ab and carrying out the summations. At large N , this leads to In the last step we used the rescaling of the dimensionful couplingsᾱ 1/2 , as given in Eq. (3.24). By comparison with the left-hand side, we extract the anomalous dimensions which results from the solution of an algebraic equation.
At 2nd order of the expansion, we obtain where S 1 to S 4 are given by N (g, b)) , N (g, h)) , N (a, h)) , . (3.37) Again focusing on the single-trace terms only and setting A 1ab = a 1 δ ab and B 1ab = b 1 δ ab , at large N one finds In the last step we used the rescaling of the dimensionful couplingsᾱ 1 ,ᾱ 2 andβ, see Eqs. (3.24) and (3.25).
We therefore obtain the following beta functions, All three beta functions reflect the exchange symmetry A ↔ B. The same symmetry determines the fixed-point structure: A fixed point must either be symmetric under the exchange α 1 ↔ α 2 , or come with a counterpart such that the pair of fixed points can be mapped into each other under this exchange. Furthermore, for α 1 = α 2 = β the system of beta functions is governed by a single equation, reflecting the fact that the model then exhibits a hidden U(1) symmetry which is otherwise broken [121][122][123][124][125][126].
Comparing with previous work, we find agreement with the beta function for the singlematrix quartic coupling reported in [13], as expected (note the difference in sign of the quartic term in the truncation). As a side remark, we mention that when extending the single-trace truncation used here by the next higher-order term, no contribution of type Tr (AB) 3 is included in the effective action since it would violate the Z 2 -symmetry of the model. Then only the beta functions for α 1 and α 2 receive further contributions in agreement with the single-matrix model limit [13].

Fixed-point structure in the symmetric limit
The system analyzed above lends itself to an enhancement of the symmetry to an A ↔ B exchange symmetry, which entails a smaller theory space since it requires setting α 1 = α 2 ≡ α, as well as η A = η B . In addition, when also α = β holds, the symmetry of the system is enhanced even further, since the model then displays a U(1) invariance where the matrices A and B behave as a vector under U(1) transformations. In these symmetry-enhanced cases we obtain the fixed-point structure in Tab. 1, where we use the convention that the critical exponents θ i are the eigenvalues of the stability matrix, multiplied by an additional negative sign, i.e., Herein, g denotes the vector of all couplings. For the U(1) symmetric fixed point C, the first critical exponent, θ 1 = 1.07, is also recovered within the U(1) symmetric theory space, where only a single quartic interaction exists. The second critical exponent, θ 2 = 0.43, encodes the relevance of U(1)-symmetrybreaking perturbations. Accordingly, we conclude that the enhancement of the symmetry in the IR requires tuning, as the fixed point is not IR attractive from outside the U(1) symmetric theory space.
At α = 1.5, we observe a singularity of the flow arising from the non-polynomial structure for the anomalous dimension, cf. Eq. (3.35) beyond which we have discarded any further zeros of the beta functions. These are characterized by anomalous dimensions beyond the  Table 1: Fixed-point values of the couplings and respective critical exponents in the symmetric limit where α 1 = α 2 ≡ α holds. By SP (saddle point) we denote fixed points with one relevant direction, which can act as either IR or UV fixed points of the flow. regulator bound η max = 1, which is required for a well-defined regulator, see [15,127] for more details. Further, they exhibit very large critical exponents O(10), and therefore fail the a-posteriori-check of the self-consistency of our truncation. The latter is based on the assumption of near-canonical scaling, and therefore requires deviations of the critical exponents from the canonical dimensions to be O(1). The fixed-point candidates are illustrated in Fig. 5. The critical exponents for these are compatible with our assumption of nearcanonical scaling. Further, the anomalous dimension η is relatively small. Accordingly, we only find changes of the fixed-point properties at the percent level when we compare to the perturbative approximation, where the ηs arising from loop factors are neglected.

Fixed-point structure in the asymmetric case
The more general theory space for the ABAB-matrix model does not require the exchange symmetry between A and B to hold at each point during the flow. Accordingly, there can even be fixed points outside the symmetry-enhanced theory space in which α 1 = α 2 . Such fixed points appear in pairs which are mapped into each other under the mapping A ↔ B. Further, the symmetry-enhanced fixed points from Tab. 1 are also necessarily fixed points of the more general flow. Within the larger theory space, they can be IR attractive or repulsive, i.e., symmetry enhancement can occur as a natural consequence of the flow. In Tab. 2, we list fixed-point candidates and their properties. We highlight the symmetry-enhanced fixed points from Tab. 1 in italics; the additional critical exponent which indicates whether the symmetry-enhancement is an automatic consequence of the flow (which requires θ < 0) is indicated in bold. As expected, fixed points with α 1 = α 2 come in pairs which map into each other under α 1 ↔ α 2 , cf. third and fifth as well as fourth and sixth line of Tab. 2. It is interesting to observe that at these fixed points, one of the two sectors features no self-interactions, and can therefore be integrated out exactly in the path integral.
An inspection of Tab. 2 highlights that symmetry-enhancement requires fine-tuning, as the additional critical exponent that characterizes the flow of the symmetry-enhanced fixed points in the larger theory space is positive for both cases. Accordingly, the surface α 1 = α 2 is likely to be an IR repulsive surface, since both fixed points inside it are. Given that these two fixed points differ by one in the number of their relevant directions, there is a separatrix linking the two, providing a complete trajectory within the symmetry-enhanced fixed point that links a nontrivial universality class in the UV to a nontrivial universality class in the IR. In addition, for α 1 = α 2 = β the symmetry of the system is further enhanced. The corresponding surface is strongly IR repulsive since the fixed point inside it has three relevant directions. Thus any quartic perturbation away from this symmetry enhancement grows under the RG flow. The fixed points and separatrices which serve to divide the space of microscopic couplings into distinct phases, are shown in Fig. 6.
As expected, we recover θ = 1.07 for the fixed point with α 1 * = 0 = β * (and α 2 * = 0 = β * , respectively). Both of those correspond to the interacting fixed point in the single Hermitian matrix model. In the single-trace approximation, the critical exponent converges to θ = 1 from above. As discovered in [16], an alternative prescription for the critical exponents, in which η = const is implemented, improves this result to θ → 0.91, whereas the exact result, corresponding to the double-scaling limit, is θ = 0.8. In the present paper, we work with a simple truncation and do not aim at accuracy in the scaling spectrum, therefore we only provide the critical exponents evaluated according to Eq. (3.42).

Review of the phase diagram of the ABAB matrix model
The two-matrix model with ABAB interaction has been exactly solved in the large N limit [105,106] using the character expansion method [128][129][130][131][132]. Its phase diagram has been analyzed there and in particular in relation to CDT in 2 + 1 dimensions in [104,113,114]. Indeed, our results look strikingly similar to theirs, constituting a strong test of the functional RG technique. In particular, we emphasize that we obtain our phase diagrams in the simplest possible truncation of the flowing action. Thus the qualitative agreement is a very promising sign for the performance of the method. Regarding quantitative results, we can compare, e.g., the critical exponents at the α 1 = 0 = β fixed point to that of the single-matrix model. As in studies of the single-matrix model in [13], the value for the relevant critical exponent deviates from the exact result corresponding to the string susceptibility γ str = − 1 2 by about 34 %. At the α 1 = 0 = α 2 fixed point the relevant critical exponent differs by about 17 % from the exact result γ str = 1 3 [106]. For fixed points C and E in Tab. 2 the relevant critical exponent is off by 21 % compared to γ str = − 1 3 [105,133]. Finally, at the tri-critical point corresponding to the fixed point H in Tab. 2, our result departs by 7 % from the exact value γ str = 0 [106,107,134]. Therefore, an extended truncation is called for to obtain quantitative results for the various universality classes. Additionally, it is an important open question what the impact of higher-order operators, in particular multi-trace ones, is on the full structure of the phase diagram.
Regarding the interpretation of the various phases, we follow the literature, in particular [104,106]. The different phases can be distinguished by exploring the behavior of order parameters. For instance, the scaling behavior of the loop average serves as an order parameter which allows to relate the length L of a loop of A/B-type of links to the enclosed area A. To make similar arguments solely using FRG technology, the composite-operator renormalization construction [135,136] developed in the context of the Asymptotic Safety program in [137][138][139], would have to be carried over to the application of the FRG in the discrete quantum-gravity context. For the symmetric model (α 1 = α 2 ) with ABAB-interaction term, it is possible to identify two distinct continuum phases. In [104,106], the first phase is identified with the line which connects the points D and C in Fig. 5, which is a separatrix in our RG flow. Indeed, initial conditions for the RG flow along that lines are the only ones that lead to the fixed point D as the IR fixed point. At this fixed point, β = 0, which means that one faces two decoupled single-matrix models. A fixed point at finite quartic coupling with one relevant direction is well-known to encode the double-scaling limit [13,118]. Thus, its geometric properties match those of 2d continuum Liouville quantum gravity [1,140,141] which is a model describing a c = 0 conformal field theory coupled to 2d Euclidean quantum gravity, with Hausdorff dimension d H = 4 and the scaling relation A ∼ L 2 [106], where A is the area of a closed loop and L a characteristic length scale of the loop.
The second phase in the symmetric model (α 1 = α 2 ), is identified in [104,106] with the line which connects points B and C. In our RG flow, this line again corresponds to a separatrix and is the unique RG trajectory that ends in fixed point B. For this phase, the geometric properties are not the same as of the other phase in the symmetric model. In spite of also having c = 0, the critical behavior is different and one finds the anomalous scaling A ∼ L 4 3 , which is that of a branched polymer [142,143]. These results indicate that touching interactions are abundant in the second phase and lead to a fractal structure of the emergent geometry.
Further, we point out that both phases are separated by a phase transition which is marked by the point C in Fig. 5. As noted in [104], it is not clear if and how this point can be given an interpretation in terms of three-dimensional geometry and CDT in particular. However, as stated in [105][106][107][108], it can be related to a conformal field theory with central charge c = 1 (i.e., a free and massless boson) coupled to 2d Euclidean gravity.
Finally, fixed point A is the free fixed point and therefore only exhibits trivial scaling behavior.
One obtains a somewhat richer phase structure for the general case α 1 = α 2 , see Sec. 3.3.2. Its geometric properties are discussed in [105,114]. As far as the relation to CDT in 2+1 dimensions is concerned, however, only the critical line for which α 1 = α 2 holds (corresponding in Fig. 6 to the line which connects the points H and G) is relevant [114]. Based on this, it has been suggested [114] to identify this phase of the matrix model with that of CDT in dimensions 2 + 1 [10,117] which, as shown by numerical evidence [115], corresponds to extended three-dimensional spacetimes with Hausdorff dimension d H = 4 for the spatial slices, in line with the value found for 2d Euclidean gravity. Notice, however, that without a more detailed grasp of the transfer matrix and thus of the Hamiltonian operator, it is not possible to more precisely characterize the interactions between the adjacent spatial slices and to understand how it analytically encodes the time evolution and thus the generation of extended three-dimensional geometries, so far only numerically observed by means of Monte-Carlo simulations in CDT [11,115].

Discussion and outlook
In this paper, we have pursued the goal to further develop functional Renormalization Group techniques for matrix and tensor models, with a particular focus on multi-field models that could constitute the key to incorporate causality into the setting. Our leading-order FRG analysis has resulted in a phase diagram in good agreement with the results presented in [106], providing a strong test of the functional RG approach to matrix/tensor models and highlighting that already at very low truncation order, key physics aspects are well captured. To quantitatively characterize the universality classes requires extended truncations, as is obvious, e.g., from a comparison of the 2d Euclidean gravity scaling exponent that we recover in a limiting case which differs from the exact result by 34 %. In the future, a more accurate determination of scaling exponents will allow to bridge the gap to continuum techniques to understand whether the Reuter universality, in particular in continuum studies with Euclidean signature [144][145][146][147][148][149][150][151] or with foliation structure [152][153][154][155][156], can be recovered from the matrix/tensor model setting. At the same time, this would also allow to understand whether another gravitational universality class, that related to the tentative asymptotically free fixed point in Hořava-Lifshitz gravity [39,[157][158][159][160][161][162][163], can be recovered from the present setting, see [164][165][166][167][168][169][170][171][172] for studies on the relation between CDT and Hořava-Lifshitz gravity. Indeed one would expect that matrix/tensor models which encode causality should be rich enough to encode various continuum universality classes within their phase diagram. A comparison of the scaling spectrum between continuum and discrete settings can provide a check of this expectation. Such an improvement in accuracy is expected to arise both from an improvement in the truncation. The critical exponents at the fixed points we uncovered here are not incompatible with a truncation principle that follows canonical scaling: As in other matrix and tensor models, an increasing order in powers of A and B as well as number of traces of an interaction decreases the scaling dimension of the associated coupling. For fixed-point candidates for which the deviation of the critical exponent from the canonical scaling dimensions is O(1), as it is in our case, robust truncations can be constructed by neglecting those interactions which are canonically highly irrelevant ones. Further, setting up a coarse graining in matrix size N leads to a violation of the U(N ) symmetry associated with each of the matrices (where N is a UV cutoff). The resulting Ward-identities have first been solved for matrix models in [8], where it was used that they become trivial in the tadpole-approximation. The latter is well-suited to characterize the fixed point in matrix models for 2d quantum gravity and might therefore also be viable in the present case of the ABAB matrix model. Beyond the tadpole approximation, the solution of the Ward-identities has been explored in [173][174][175].
Beyond the characterization of the universality class in terms of the critical exponents, an understanding of the emergent geometries is desirable. Within the FRG, the compositeoperator formalism [135][136][137][138][139] appears well-suited to provide access to, e.g., the scaling of geometric quantities. At the purely formal level, the evaluation of the corresponding flow equations for matrix/tensor models is not expected to be more challenging than for the flow equation for the effective action itself. The main challenge therefore lies in the identification of suitable operators, where one would expect geometric information to be encoded in higher-order operators in the tensor/matrix formalism.
The relation of the ABAB matrix model to a causal model of dynamical triangulations is an example to highlight a promising route to impose causality in matrix/tensor models: To take into account the presence of both timelike and spacelike edges in a triangulation, a multi-field approach seems indicated, see for instance [176]. The present paper as well as the recent application of the FRG to a model of non-commutative spacetime with two fields [111] highlights that the FRG is well-suited to study such models. The present paper lays the ground for extensions to higher-rank models.
Beyond the encoding of causality, multi-field models could also account for the interplay of quantum gravity with matter [27,177,178], providing yet another motivation to extend the present study to tensor models in the future.