Towards anisotropic cosmology in group field theory

In cosmological group field theory (GFT) models for quantum gravity coupled to a massless scalar field the total volume, seen as a function of the scalar field, follows the classical Friedmann dynamics of a flat Friedmann-Lema\^itre-Robertson-Walker (FLRW) Universe at low energies while resolving the Big Bang singularity at high energies. An open question is how to generalise these results to other homogeneous cosmologies. Here we take the first steps towards studying anisotropic Bianchi models in GFT, based on the introduction of a new anisotropy observable analogous to the $\beta$ variables in Misner's parametrisation. In a classical Bianchi I spacetime, $\beta$ behaves as a massless scalar field and can be used as a (gravitational) relational clock. We construct a GFT model for which in an expanding Universe $\beta$ initially behaves like its classical analogue before"decaying"showing a previously studied isotropisation. We support numerical results in GFT by analytical approximations in a toy model. One possible outcome of our work is a definition of relational dynamics in GFT that does not require matter.

In cosmological group field theory (GFT) models for quantum gravity coupled to a massless scalar field the total volume, seen as a function of the scalar field, follows the classical Friedmann dynamics of a flat Friedmann-Lemaître-Robertson-Walker (FLRW) Universe at low energies while resolving the Big Bang singularity at high energies. An open question is how to generalise these results to other homogeneous cosmologies. Here we take the first steps towards studying anisotropic Bianchi models in GFT, based on the introduction of a new anisotropy observable analogous to the β variables in Misner's parametrisation. In a classical Bianchi I spacetime, β behaves as a massless scalar field and can be used as a (gravitational) relational clock. We construct a GFT model for which in an expanding Universe β initially behaves like its classical analogue before "decaying" showing a previously studied isotropisation. We support numerical results in GFT by analytical approximations in a toy model. One possible outcome of our work is a definition of relational dynamics in GFT that does not require matter. A major challenge for discrete approaches to quantum gravity is the derivation of an effective (emergent) continuum description which can be compared with classical general relativity or more general gravitational theories. The challenge arises on many levels, for instance in recovering the usual notions of a spacetime manifold from combinatorial structures [1]; recovering an effective description in terms of coordinates and restoring the continuum notion of diffeomorphisms or coordinate changes [2]; and understanding the intricate interplay between a continuum and semiclassical limit. Deriving such a description, however, is crucial for understanding the phenomenology of such quantum gravity theories and ensuring their compatibility with observation given that, e.g., new fifth-force degrees of freedom at low energies would have to be compatible with tight experimental bounds [3]. A common approach in this situation is to restrict to situations of high symmetry, in particular spatially homogeneous cosmology or spherically symmetric black holes. While by assumption they no longer include all degrees of freedom, symmetry-restricted models would be expected to capture at least some phenomena of the underlying theory (as they do in classical general relativity), while also connecting directly to phenomenology given the obvious relevance of cosmological and black hole spacetimes. A prominent example is loop quantum gravity (LQG), whose cosmological sector has been studied in loop quantum cosmology [4] while there are also a number of effective black hole models including LQG discretisation effects [5].

CONTENTS
Effective cosmological models have recently been constructed in the GFT approach to quantum gravity [6], itself closely related to the spin foam definition of LQG dynamics [7] and to matrix and tensor models [8]. Just as LQG, GFT models are fundamentally defined in terms of discrete, combinatorial structures, interpreted loosely as "quanta of spacetime". One can associate geometric notions such as areas, volumes and angles to these degrees of freedom by incorporating concepts from LQG, but there is no simple way of giving these an effective continuum interpretation, in particular given that the conceptual status of a given discretisation in GFT is similar to a Feynman graph in quantum field theory, i.e., only one term in an infinite expansion. What can be done relatively straightforwardly, however, is to derive dynamical equations for global observables such as the total (spatial) volume of a certain geometry, which can then be contrasted with globally homogeneous cosmological models. These equations are usually derived for certain GFT states whose properties make them good candidates for spatially homogeneous geometries. After some prior groundwork [9] a breakthrough in this line of research came when, in GFT models for quantum gravity coupled to a massless matter scalar field, a "relational" volume observable (corresponding to the volume of space for a given value of the scalar field) was shown to satisfy the Friedmann dynamics of general relativity at low energies while also replacing the classical Big Bang singularity by a bounce [10]. Similar results have been obtained using different methods and from different starting points [11][12][13][14], emphasising their robustness: one can understand the main properties of the resulting cosmological dynamics from classical solutions for a single field mode. While important for the phenomenology of GFT and for connecting to approaches such as loop quantum cosmology, these results have so far been restricted to the case of a flat homogeneous, isotropic Universe. 1 Some studies have included anisotropies perturbatively [16,17] showing that they decay leading to isotropisation, but there is so far no characterisation of, e.g., an anisotropic Bianchi cosmology.
Here we take the first steps towards the study of anisotropic cosmologies in group field theory, focussing on the simplest possible case of Bianchi I cosmology with local rotational symmetry so that two out of the three directional scale factors are taken to be equal. There are at least two, initially quite separate, challenges involved in this extension of past work. The first is to find a characterisation of anisotropies in group field theory, i.e., to define an observable that can distinguish isotropic and anisotropic geometries and quantify the amount of anisotropy. Here the key idea we use is the Misner parametrisation of Bianchi models (see, e.g., [18]; we will also review this below) in terms of a volume degree of freedom and two relative anisotropy variables, the Misner variables β ± . In a classical Bianchi I model β ± behave as free, massless scalar fields in a flat FLRW geometry 2 , which have already been studied in GFT. On the other hand, the discreteness of geometry in GFT means that we cannot simply take over a continuum definition of anisotropy, so the construction of an analogue β ± variable requires careful thought. The second challenge is to understand which simplifying approximations used in past work need to be relaxed in order to allow for anisotropies in the effective description. For instance, while the work of [10] only used "isotropic" states interpreted as describing simplicial building blocks for which all faces have equal area, it is known (see, e.g., [13]) that this microscopic restriction to isotropy is neither necessary nor sufficient to obtain the correct (flat FLRW) Friedmann dynamics: the more relevant assumption is to restrict to a single field mode in the Peter-Weyl expansion in representation data. Hence, in order to describe anisotropic geometries, multiple Peter-Weyl modes must be taken into account, but it is not clear how many (and which) modes are needed to capture physical anisotropies.
In this paper we show how to tackle the first challenge; we define an β ± analogue with a clear geometric interpretation, quantum ambiguities that disappear for large quanta, and correct physical properties -constant velocity and hence linear evolution -at least for a certain cosmological period of time, before the isotropisation observed in [16,17] sets in and the anisotropy disappears. The second challenge is partially addressed, given that the β ± dynamics partially match expectations from classical relativity, but the observed isotropisation does not correspond to a classically expected behaviour and, more importantly, anisotropies do not backreact on the effective Friedmann equation as expected. This suggests that while our constructions will be useful for future work, our model needs further refinement to reproduce the correct physics of a classical Bianchi Universe.
Any monotonically evolving quantity in a cosmological model can be used as a relational clock: all other dynamical variables can be written, at least in principle, as functions of this "clock". In a vacuum Bianchi I model, the anisotropy variables β ± have this property and hence, in contrast to what is often done in quantum cosmology and in particular in GFT, no coupling to matter would be needed to be able to express the dynamics in relational terms. The fact that we have defined a new quantity with monotonic evolution in GFT cosmological models hence raises the possibility of defining relational evolution in GFT without adding matter fields, which might help in understanding the "problem of time", or possible dependence of dynamics on the choice of clock, in GFT (see [20] for some work on this issue in models with multiple possible clocks).
The remaining parts of the paper are structured as follows. In section II we review basic ideas of the GFT formalism, possible definitions of a canonical Hilbert space quantisation, and application to cosmology: we show how one can derive effective Friedmann equations by restricting to a single field mode and neglecting interactions. Readers familiar with GFT cosmology may skip this review section. Similarly, section III is a review of classical FLRW and Bianchi cosmologies written in relational terms, using a scalar field clock; this is the classical theory that any effective description of GFT can be compared to. Section IV includes the main new results: we motivate the introduction of an effective β ± variable used to characterise anisotropies in GFT. We then propose models based on a few Peter-Weyl modes and study the effective dynamics of both the anisotropies and the spatial volume, comparing both with the dynamics of general relativity. We also propose a simplified "toy model" in which some of our main results, in particular the linear growth in anisotropy which matches classical expectations, can be derived analytically rather than numerically as in the main part. We conclude in section V. An appendix gives details on the classical and quantum geometry of tetrahedra as used in LQG and GFT.

II. SHORT REVIEW OF GFT COSMOLOGY
In this section we briefly summarise past work on deriving effective cosmological dynamics from group field theory. In this past work, effective Friedmann equations were obtained after truncating the full dynamics and choosing simple GFT states, following two different approaches. The two approaches, which we will call algebraic and deparametrised, will be introduced in section II B.

A. Basics of group field theory
We are interested in GFT models for simplicial gravity coupled to a (free, massless) scalar field χ. In such models one defines a group field ϕ whose arguments are d elements of a Lie group G (hence the name "group field") and a real variable corresponding to the scalar matter field χ ∈ R, where K can be R or C. When applied to four-dimensional quantum gravity, d = 4 and one usually takes G to be the local gauge group of general relativity: G is typically SO(3, 1) or SL(2, C) in the Lorentzian case, SO(4) or Spin(4) in the Euclidean case, or their rotation subgroup SU (2) which is the gauge group of loop gravity (or the Ashtekar-Barbero formulation of classical general relativity [21]). The last choice is the one we will use here. To implement a notion of discrete gauge invariance in the resulting simplicial gravity description, one usually requires invariance of the field under the right diagonal group action, The action has the general form where for a real fieldφ = ϕ. Here dg stands for an integration over d copies of the group, using the Haar measure normalised to unity. The action is therefore split into a quadratic part and an interaction part V . The kernel K is assumed to respect the symmetries associated to a minimally coupled massless scalar field on a curved background, namely shift (χ → χ + c) and sign reversal (χ → −χ) symmetries. For this reason, K cannot depend explicitly on χ, but it is in general a differential operator in χ, which does not involve odd powers [10,22]. χ usually plays the role of a relational time variable, so that other observables are defined relative to χ, as is common in many approaches to quantum cosmology, in particular in loop quantum cosmology. The kinetic term can be written as an expansion in derivatives with respect to χ, which is usually truncated after the second term: starting from the simplest K that only includes a constant "mass term" suggested by the relation to spin foam models [23], a Laplacian term is generated by radiative corrections [24]. One could stop here, given that no higher derivative terms are required, or make the weaker assumption that higher derivatives are present but suppressed, i.e., |K (2n) /K (0) | |K (2) /K (0) | n for n > 1. For concreteness, we will follow previous work (see, e.g., [9,25]) and assume that K has the minimal form where m and M are coupling constants and ∆ g I is the Laplace-Beltrami operator acting on the I-th group argument. Evidently, this corresponds to (4). To construct an interaction term in simplicial gravity models, one can think of the group field ϕ as representing a (d − 1)-simplex. The interaction term then describes the gluing of such simplices to form d-dimensional structures, here a d-simplex, similarly to what happens in tensor models [8,26]. This means that an appropriate interaction term consists of products of fields that are paired according to a pattern which encodes the combinatorics of a d-simplex. In four dimensions, we could for example choose (here for a real field) where λ is a coupling constant, (dg) 5 means integration over G 20 , and V (g 1 I , ..., g 5 I ) is a product of ten Dirac delta distributions on the group, imposing appropriate matching between group elements appearing as arguments of the fields ϕ(g j I , χ), in order to encode the pattern of gluings needed to form a four-simplex out of five tetrahedra (see right panel of figure 1). Such an interaction should allow the structure of four-dimensional spacetime to emerge from the dynamics of the theory.
We will now fix G d = SU (2) 4 . Then the field ϕ(g I , χ) can also be represented as a four-valent spin network node. The g I are associated to the links dual to the faces of a tetrahedron (see left panel of figure 1), and the matter field χ "sits" on the node dual to the tetrahedron.
It is often convenient to use the Peter-Weyl theorem to express the group field as where the D (j I ) m I ,n I (g I ) are Wigner matrices and ϕ j,ı m (χ) ≡ ϕ j 1 ,...,j 4 ,ı m 1 ,...,m 4 (χ) are complex functions 3 . The intertwiners I j,ı n ≡ I j 1 ,...,j 4 ,ı n 1 ,...,n 4 arise because of the invariance under group multiplication from 3 If the group field is real, the complex Peter-Weyl coefficients satisfy the reality condition [12,27] ϕ j,ı m (χ) = (−1) Graph of five vertices, representing a four-simplex as gluing of five tetrahedra (three-simplices). Here each field is represented by a node with four legs, and a shared link corresponds to matching group arguments.
the right (2). Recall that intertwiners are elements of the Hilbert space of the tetrahedron (or four-valent node) H 4 = Inv 4 I=1 H j I , where each H j I corresponds to the Hilbert space of an irreducible unitary representation of SU (2) (see appendix A). The sums are over representations j I ∈ N 0 /2 and magnetic indices m I , n I ∈ [−j I , j I ], while ı labels the possible intertwiners for each set of given spins. One can picture each pair (j I , m I ) as living on a link emerging from the node, while ı lives on the node itself.
The main motivation behind these definitions is the connection with spin foam models; one can understand GFTs as a more fundamental quantum field theory-like framework (or quantum gravity theory) into which spin foam models coming from LQG can be embedded. Namely, the Feynman expansion of the GFT partition function generates a sum over graphs Γ that can be seen as dual to simplicial complexes (here again for a real field), where n V (Γ) is the number of interaction vertices and A Γ is the Feynman amplitude. The Feynman integrals over momenta are discrete sums (as the space on which the field theory is defined is compact) over group representations and intertwiners, associated with faces and edges of the twocomplex dual to Γ. The GFT Feynman amplitudes are exactly what is usually called spin foam amplitudes [23,28]. In the LQG spin foam approach, one often makes a choice of a fixed twocomplex, but a sum over two-complexes would be necessary if one hopes to capture the continuum dynamics of quantum gravity. The Feynman expansion (9) provides precisely such a sum over two-complexes, with relative weights determined by λ, and hence a generating functional for the covariant quantisation of LQG. By choosing K and V in the GFT action (3) one can generate models which are related in a precise way to different spin foam models.
We should mention that just as for the kinetic term K, the choice of interaction term(s) V may also be guided by considerations beyond the correspondence with simplicial gravity. For instance, additional interaction terms could be generated by renormalisation group flow (see, e.g., [29] for an analysis of models with G = U (1)), or included in a general effective field theory description.

B. Canonical quantisation of GFT
The traditional way of thinking of quantum GFT is through the path integral as given in (9), and its expansion into Feynman graphs and amplitudes. As we discussed, this path integral directly connects to the covariant (spin foam) setting of LQG. More recently however, a lot of work has focused on the canonical quantisation of GFT, with the main goals of connecting to the canonical setting for LQG and -most importantly for us here -extracting effective cosmological dynamics which are more easily defined in a canonical setting, as they are in LQG where loop quantum cosmology has so far only been derived through canonical methods.
Two main approaches to canonical quantisation have been established in the GFT literature: one based on a kinematical Fock space of nondynamical spin network-like states on which dynamical equations are imposed in a suitable (usually mean-field) approximation [30], and one where a time variable is selected before quantisation and used to directly obtain a physical Fock space and physical (relational) Hamiltonian [12]. These two types of quantisation are to some extent analogous to two approaches in usual canonical LQG: on the one hand, in LQG one can define a kinematical Hilbert space whose states -though SU (2) invariant and satisfying the (spatial) diffeomorphism constraint -do not yet satisfy any dynamics. Physical states would only emerge after one has also implemented the Hamiltonian constraint (see, e.g., [31]). Since this can usually not be done exactly, one can try to implement dynamics approximately, for example using coherent states to study a semiclassical effective theory [32][33][34]. On the other hand, one can introduce matter (often pressureless "dust" scalar fields) to define a gauge-fixed classical theory, in which the matter is used as a preferred standard of time. This deparametrised approach leads to a physical Hilbert space since the Hamiltonian constraint can be solved for the momentum canonically conjugate to "dust time" [35,36]. While this reduced Hamiltonian approach breaks time reparametrisation invariance by requiring a preferred time coordinate before quantisation, it automatically describes physical quantum states which satisfy dynamics, and allows bypassing the unresolved problem of implementing the Hamiltonian constraint on the kinematical Hilbert space.
In GFT, there is no Hamiltonian constraint with the same status as in canonical quantum gravity, i.e., related to diffeomorphisms in time. The two approaches we will review nevertheless share the same two different starting points of either imposing dynamical equations on an abstract Hilbert space, or directly defining a physical Hilbert space by choosing a matter clock. While in general these two types of quantisation cannot be equivalent, in the application to cosmology they lead to very similar effective dynamics, essentially since one requires semiclassical states with small quantum corrections to classical GFT solutions.
vacuum |0 satisfyingφ J |0 = 0. The Fock space contains fundamental quanta (pictorially "atoms of space" [38]) created and annihilated by field operators, which are interpreted as equivalent to LQG spin-network vertices. As in the kinematical Hilbert space of LQG, there is no notion of dynamics so far: this will only be implemented "on average" later on, similarly to what happens in condensed matter physics. The Fock vacuum |0 is a state with no topological nor geometrical meaning and information. Then the one-quantum Hilbert space has states interpreted as (open) spin network vertices decorated with four Peter-Weyl labels (or group elements) and a scalar χ. At this stage one usually defines general one-body operators (see, e.g., [30]) where O J is the LQG matrix element of the desired operator evaluated on a single spin network node, e.g., the volume 4 V j,ı which plays a particularly important role in cosmology. To give explicit formulae, the number and volume operators read Notice that these operators are defined "relationally", in the sense that they are given as functions of the matter field χ, associated to the creation and annihilation operators. At the kinematical level, this χ-dependence does not represent time evolution: operators defined for different χ are independent. One then usually assumes that equations of motion are satisfied on average, where S is the action of our GFT model. Specifying suitable coherent states to compute expectation values will then give effective cosmological equations, as explained in the next section.

Deparametrised approach
If, on the other hand, one wants to follow a deparametrised approach, one can work with a real GFT field. Working in the spin representation, equation (3) then reads [12] 4 Concretely, after choosing a basis of intertwiners ı that diagonalises the LQG volume operator on a spin network node, one associates a volume V j,ı , given by the eigenvalue of this operator, to the quanta generated by ϕ † J (χ). See appendix A for more details on the spectrum and eigenstates of this LQG volume operator. These properties can be derived from a more general setting of quantisation of tetrahedra in terms of SU (2) recoupling theory [39][40][41], without working necessarily within LQG.
χ is defined as in (4). In this paper we will assume (5) so that K J = m 2 − ∂ 2 χ − M 2 I j I (j I + 1) (using the fact that the Laplace-Beltrami operator acts as a Casimir on Wigner matrices, ∆ g D (j) mn (g) = −j(j + 1)D (j) mn (g)), but we can keep K (0) J and K (2) J general for the time being. Now one can define the conjugate momentum π J (χ) to the group field ϕ J (χ), and the Legendre transform of the Lagrangian L with respect to χ gives a relational Hamiltonian which determines the dynamics of any observable O through Poisson brackets dO/dχ = {O, H}.
The field and its momentum are promoted to operators with the usual canonical equal-time commutation relation The key difference with the previous approach is that these operators already satisfy dynamical equations, implemented through the Heisenberg equations of motion. This has a cost: we needed to specify our time variable once and for all from the very beginning to define the conjugate momentum of the field and the relational Hamiltonian.
One can now define creation and annihilation operatorsâ † J andâ J as in any bosonic quantum field theory (not to be confused withφ † J andφ J ) with their own (equal-time) commutation relations, and use these to construct a physical Fock space. This space is "smaller" than the one introduced in the earlier algebraic approach since the states are already interpreted as physical states, not subject to any constraints. 5 Dynamics for any operatorÔ are given by relational Heisenberg equations The relational Hamiltonian operator for the free theory (i.e., for V [ϕ] = 0) can be expressed as a sum of single-mode HamiltoniansĤ = JĤ J , and theĤ J can themselves be written in terms of creation and annihilation operators. For a mode such that K where J . For the case of interest in the rest of this paper, where the kinetic term is (5) so that K J is as given below (15), Only modes for which the argument of the square root is positive will correspond to such a Hamiltonian. The coupling m j only depends on the spin labels, so magnetic and intertwiner indices can be dropped in many expressions. The Hamiltonian (20) is a squeezing operator which does not leave the Fock vacuum invariant, but creates pairs of excitations with opposite values for the magnetic indices. Expanding space arises from the instability of the Fock "vacuum", realising a type of "geometrogenesis" (the term was coined in [43] and then used in the GFT literature [9,44,45]). The rate of squeezing or expansion is determined by |M J |; for models in which |M J | takes a maximum for some J (such as, in our case, for the modes of lowest j I ), these modes will always dominate asymptotically [25]. This instability is to be contrasted with modes for which K (0) J and K (2) J have the same sign; these modes are stable with constant particle number and quickly become insignificant compared to the unstable modes, which is why they are usually ignored. Of course, this also means that to obtain a realistic cosmology we have to assume that at least one mode has a Hamiltonian of squeezing type (20).
In order to extract the simplest cosmological interpretation, we are only interested in the total particle number and volume operators, here given bŷ where we distinguish the notation from (13) by underlining them to keep in mind that these are different operators from the ones defined in the algebraic approach. We can however compare the two approaches at the level of expectation values, as illustrated below.
The operators we are interested in here are always a sum of single-mode operators. Such onebody operators are extensive and often simplify because one is only interested in certain modes.

C. Emergent FLRW Universe from free theory and single mode
In order to link this theory with cosmology, some simplifying assumptions are needed. First of all, the cosmological sector of GFTs often deals with regimes in which interactions can be neglected, so we will only consider the kinetic term of the action (i.e., the free theory). This approximation is often interpreted [9,12] as corresponding to homogeneity given that interactions between spinnetwork nodes (and correlations between the GFT quanta) are negligible. Given the instability of the theory, such an approximation can only hold for a finite amount of time before the number of quanta is too large [10].
Moreover, we now focus on a single Peter-Weyl mode. This restriction is motivated by computational simplicity and the fact, mentioned above, that often one mode quickly dominates dynamically so that this approximation becomes better and better with time. There is evidently a certain clash with the first approximation of negligible interactions, which gets worse over time. In any case, this second approximation means that a cosmological spacetime expands or contracts by modifications to the combinatorial structure of the spin network (i.e., by changing number of the GFT quanta), rather than by changing the spin labels on the network (i.e., transitioning between GFT quanta of different spin representations). Furthermore, one can use the insights gained from loop quantum cosmology, where all the spins are usually fixed to only one value [4], to motivate this single-mode restriction in GFT. Such assumptions behind loop quantum cosmology models are often motivated by a suggestion that cosmological expansion or contraction is indeed realised by a changing graph structure in full LQG [19].
One then also needs to specify a particular type of states, and here coherent states are used to implement a notion of semiclassicality similar to what is often done in quantum cosmology: in a macroscopic Universe quantum fluctuations over expectation values should be small. Different choices of coherent states are examined with respect to this criterion, e.g., in [13].

Algebraic approach
The next step in the algebraic approach is to define states in which (14) is to be evaluated. The simplest choice is given by field coherent or "condensate" states, whose key property iŝ This property allows replacing (after normal ordering) all field operators by the collective variable σ J (χ) (sometimes called condensate wavefunction), which explicitly depends on χ; the quantum equation of motion (14) reduces to the classical GFT equation of motion. This mean-field approximation can be interpreted as the idea that cosmology arises from the "hydrodynamics of quantum gravity" [38]. Since all quanta are characterised by a single quantum state, it has been conjectured that these states may correspond, when coarse-grained, to homogeneous cosmological spacetimes.
In this mean-field approach one can then derive effective dynamics for expectation values of the operators of interest, such asN (χ) and in particularV (χ), basically ignoring fluctuations: starting from the definitions (13) and making use of the property (23), one finds the expectation values (also assuming that only a single Peter-Weyl mode contributes) The expectation values (24) can now be evaluated for a solution of the free theory K J σ J (χ) = 0, given that we have decided to neglect interactions. In particular, we will again assume the singlemode kinetic kernel Using the notation of [25], the general solution to this equation of motion then reads where the coefficients generally depend on the mode J. Substituted into (24), one finds the volume Crucially, this form for the volume corresponds to a cosmological solution that satisfies Friedmannlike dynamics while also containing a cosmological bounce. Even though the original work [10] was based on equilateral tetrahedra which were supposed to encode isotropy, one can retain four different spins j = (j 1 , j 2 , j 3 , j 4 ) for the faces of the building blocks; the key ingredient for this result lies in the single-mode restriction alone.
To make this explicit, one can introduce the "GFT energy" , which is a conserved quantity associated to χ translations (but whose interpretation was not specified in [10,25]). One can easily show that (27) satisfies an effective Friedmann equation where Q J := 2m j (α + J α − J ) is another conserved quantity, associated with the U (1) symmetry of the complex GFT. In [10], it was proposed that this same quantity also plays the role of conjugate momentum to the scalar field, so that by defining an effective energy density ρ = Q 2 J /(2V 2 ) and ρ c = m 2 j /(2V 2 j,ı ), one can rewrite (28) in a form similar to loop quantum cosmology with a term −ρ/ρ c , plus an extra term proportional to the GFT energy. In these terms, ρ c would (for a fixed mode) represent a universal upper bound for the density ρ, which allows a bounce at a minimum nonsingular volume, just like in loop quantum cosmology. However, it is not clear why Q J would be interpreted as the conjugate momentum of the matter field, as it would make more sense to associate this role to E J (see also, e.g., [20] for more on this).
The large volume limit of (28) needs to be consistent with the classical Friedmann equation (V /V ) 2 = 12πG, which we will review in the next section. In this case, large volumes (or low energy densities) are obtained when |χ| is big enough so that we are away from the bounce. Hence agreement with the classical theory requires the identification between the GFT coupling and Newton's constant as m 2 j = 3πG. This identification should be seen as the "emergence" of Newton's constant from more fundamental GFT parameters.
Finally, we note that an effective Friedmann equation very similar to (28), with some additional corrections, emerges from the "effective relational" approach of [14] which uses different types of coherent states and different ways of approximating the full quantum dynamics. This certainly supports the robustness of the GFT cosmology approach overall.

Deparametrised approach
The idea of deparametrisation was introduced in GFT cosmology in [11] (see also [46] for related ideas) where an effective Friedmann equation similar to (28) was recovered in a simple toy model with a squeezing Hamiltonian, without requiring a mean-field approximation. In this model the Hamiltonian generates evolution with respect to a preferred matter clock. The general idea of deparametrisation, using some degrees of freedom as coordinates parametrising the others, was then used in full GFT in the Hamiltonian formalism reviewed in section II B 2 6 . Initially, coherent states were used to obtain approximate solutions for expectation values of operators of interest, but this limitation was overcome in [13] where the dynamics were solved for operators in the Heisenberg picture, without requiring a specific state. One can then choose specific states to compare expectation values with previous results; effective Friedmann equations could be found for many types of coherent states, e.g., based on the su(1, 1) algebra generated by the volume and Hamiltonian operators [48].
Since the Hamiltonian for the free GFT is quadratic, one can solve the Heisenberg equation (19) to find the (relational) time evolution of the number and volume operators (22). Then, again for a single Peter-Weyl mode, the expectation value of the volume operator latter can clearly be interpreted as the volume at χ = 0. On the other hand, a nonvanishing K(0) would imply an asymmetry with respect to χ = 0, i.e., different pre-and post-bounce scenarios. While (29) holds for any state, requiring a semiclassical description at least at late times again motivates the use of coherent states. Rather than specifying a state at all times as in (23), here one defines a coherent state at a given "initial" time χ 0 (usually χ 0 = 0), but does not assume that it will remain an eigenstate of the annihilation operator at later times. Again one convenient type of states is given by Fock coherent states defined bŷ where |σ is an eigenstate of the annihilation operator only at χ 0 . The expression for the expectation value (29) is of almost the same form as (27), but the analogue of the "GFT energy" is now fixed to m 2 j . From (30), we can see that |σ J | 2 ≡ N (0) represents the expectation value of the number operator at the initial value of our clock χ; the volume at later times is obtained from (29).
From (29), one can find an effective Friedmann equation similar to (28), One can now again define a critical energy density ρ c = m 2 j /(2V 2 j,ı ) and rewrite the last term inside the brackets in the suggestive form −ρ eff /ρ c . In general, the effective energy density ρ eff takes the form ρ eff = ρ χ (χ) + corrections, where the additional terms depend on the initial conditions, but where the matter energy density ρ χ (χ) is now defined by identifying the Hamiltonian H ≡ Ĥ J (cf. (20)) with the conjugate momentum of χ (as one would expect), i.e., ρ χ ( For our purposes, the main result is that (31) also reduces to the classical Friedmann equation in the late-time (or large-volume) limit. Moreover, it has corrections ∼ 1/V and ∼ 1/V 2 , very similar to the ones of (28) obtained by mean-field approximation in the algebraic approach.

Towards anisotropic cosmologies
We have seen that, using a restriction to a single mode, the GFT literature is rich of ways to derive effective cosmological dynamics. The results differ in numerical factors (e.g., the high energy corrections coming from the cross term 2V j, (29)), but represent the same general scenario. Whether we follow the algebraic approach (which is how the GFT cosmology programme began) or use a deparametrised point of view, we obtain comparable effective Friedmann equations, which in particular have the same large volume limit. They are also all similar to the effective dynamics of loop quantum cosmology, and describe an effective repulsive behaviour at high energies.
We now want to extend this discussion to anisotropic cosmologies. Classically, incorporating anisotropy means going from the FLRW Universe to the more general class of Bianchi models. We aim to define a GFT model that can generalise the results presented in this section at least to the simplest anisotropic cosmology (Bianchi I). To do this, we need to consider two new aspects: understanding how anisotropies modify the effective Friedmann equation, and understanding the dynamics of the anisotropies themselves. An anisotropic GFT model naively requires quanta of geometry which are non-equilateral tetrahedra. Moreover, in order to obtain a nontrivial time evolution for the anisotropies we will need to lift the single-mode restriction. Heuristically, this is because if the shape of non-equilateral tetrahedra is fixed there is no room for a dynamical notion of anisotropy. Multiple modes are required, so that the relative contributions of different shapes can change and a macroscopically "average" anisotropy can become dynamical.
Some preliminary work on anisotropic GFT cosmology was done in [16] where anisotropies were seen as perturbations. Furthermore, an anisotropic trirectangular tetrahedron was used as building block in [17] where a notion of "dynamical isotropisation" was found, but without details on the evolution of a (global) anisotropy parameter. This will in fact be the main novelty in our work: we define a measure for anisotropies given by a parameter (corresponding to a Misner-like variable β ± in classical cosmology) emerging from fundamental quantum gravity arguments.

III. RELATIONAL DEFINITION OF COSMOLOGICAL DYNAMICS
In this section we briefly recall how one can obtain a relational description of cosmological models in general relativity, starting from the Hamiltonian (Arnowitt-Deser-Misner, or ADM) form of the Einstein-Hilbert action In addition to the metric tensor of three-dimensional spatial slices q ab and its conjugate momentum π ab , there are four Lagrange multipliers: the lapse function N and the shift vector field N a . They multiply the Hamiltonian and (spatial-)diffeomorphism constraints, defined as where q = det(q), R is the Ricci scalar of q ab , D a is the spatial covariant derivative compatible with q ab and indices are raised and lowered with q ab . The constraints (33) need to vanish for physical solutions, as can be seen varying the action with respect to N and N a respectively.
We will be interested in gravity coupled to a (free, massless) scalar field χ with conjugate momentum p χ , which will serve as relational clock. Moreover we will only deal with homogeneous settings, therefore we do not have an energy gradient term for χ, and we can set N a = 0. The action then reads where now the total Hamiltonian constraint has a gravitational and a matter part (denoted C χ ), If all fields are assumed to be spatially homogeneous, the integral over space d 3 x just gives a constant, which we can set to unity. This does not play any physical role in homogeneous settings.

A. FLRW
We now specialise the above framework to a spatially flat, homogeneous and isotropic FLRW Universe. Given the metric the action is simply expressed in terms of the volume V (related to the familiar scale factor via V = a 3 ) and its momentum p V , as We want to describe dynamics in relational terms, using χ as time variable. Recall that the scalar field can act as clock because it is monotonic, since p χ is a constant of motion as {p χ , N C} = 0. The equation of motion for χ allows to choose the scalar field to be the time variable, fixing N : Now the equation of motion for the volume,V = {V, N C} = −12πGN V p V , can be reformulated in relation to χ. In particular, usingχ andV one can express p V as Substituting (39) in the vanishing of the Hamiltonian constraint (37), we find the relational Friedmann equation The large-volume limit of the isotropic GFT models introduced in the previous section is compared to (40). Note that (40)

B. Bianchi I
We now move our attention to a Bianchi I cosmology, with metric This generalises the flat FLRW metric to the case with separate scale factor a i (t) in each Cartesian direction i = 1, 2, 3. We introduce the Misner parametrisation (see, e.g., [18]) with a volume variable V = a 1 a 2 a 3 , The variables β ± represent anisotropy parameters and have their own momenta p ± . Using this parametrisation, equations (34) and (35) become where Notice the similarity between the anisotropy variables and the scalar field χ. Their contribution to the Hamiltonian is basically the same, expect for numerical factors; at least classically, a Bianchi I Universe is not different from an FLRW Universe with free massless scalar fields. Even though we still make use of a matter clock in this paper, this equivalence suggests that the anisotropy variables β ± could play the role of a clock in GFT cosmological (anisotropic) models. For now, as before, we use χ as a clock. Thus, from the equations of motioṅ we extract relational dynamics as follows. One can use the first two equations of motion to obtain p V as in (39), which can then be substituted into the vanishing of (44), to get This relational (generalised) Friedmann equation reduces to (40) in the isotropic case. Moreover, we now have anisotropy degrees of freedom, whose dynamics are obtained in a similar fashion. We use the last two equations of motion (45) to write the momenta Substituting these into the Hamiltonian constraint, we obtain (46) and (48) can also be combined into the relational equation This form has the advantage that it only relies on derivatives with respect to χ rather than canonical momenta. This equation will be used as a classical comparison for our anisotropic GFT cosmological model in the large volume limit.

C. Bianchi II
For the sake of completeness, we extend the above formalism to the Bianchi II case to see how curvature affects relational (classical) dynamics. We follow the same steps as before, but now we also have a curvature term appearing in (35) [18], and the Hamiltonian constraint (35) reads Without repeating all the steps seen in the Bianchi I case, we recall that we use the equations of motion to obtain the momenta p V and p ± . Plugging them into C = 0 one obtains and dβ ± dχ Finally, putting everything together one can write a Bianchi II generalisation of the (relational) Friedmann equation, which reduces to (49) when the last term (which comes from (3) R) is zero. Notice that in an expanding Universe this term becomes dominant at later times: the dynamics are initially close to Bianchi I and deviate only when the exponential expansion of the volume takes over. Because of this, and given that the Bianchi II scenario is the simplest Bianchi model involving spatial curvature, we will also compare our anisotropic GFT cosmology to these equations later.

IV. ANISOTROPIC GFT MODEL
In order to study anisotropic cosmologies in GFT, we need to define a notion of "anisotropy variables", analogous to the Misner variables β ± . We will study the (relational) dynamics of these variables and observe how they affect the effective Friedmann equation for the volume V (χ). We will compare these effective dynamical equations with those of Bianchi models, in particular with the simplest Bianchi I cosmology which would be the natural extension of the spatially flat FLRW Universe previously studied in GFT. As in this isotropic case, we will follow both the algebraic and deparametrised approaches. We will once again see that they give slightly different behaviours close to the bounce, but basically match otherwise.
We will study two different observables, representing the degrees of freedom of a classical Bianchi I cosmology: the volume V (χ) and "average anisotropy parameters" β ± (χ), defined as where N J (χ) is the expectation value for the particle number in the mode J. V (χ) is defined as in previous work, namely as the expectation value of (13) or (22). On the other hand, the expression for β ± (χ) is different from the usual structure of operators in GFT (cf. (12)). We think of anisotropies as determined by the shape of our geometric building blocks; these variables should be "intensive" and not simply grow with the number of quanta, therefore we divide by the total number N (χ) = J N J (χ) (see, e.g., [49] for a similar concern related to a possible scalar matter quantum operator). β ± (χ) is not an expectation value, and we do not propose any definition of operators β ± representing anisotropy; the overall 1/N (χ) could not arise from taking an expectation value. Instead, the definition (55) introduces a semiclassical notion of anisotropy only. At the end of section II, we already argued that this setup requires including multiple Peter-Weyl modes. Indeed, for only a single J in (55) β ± (χ) would clearly be constant, consistent with the idea that we would not incorporate any dynamical "shape" degrees of freedom. Moreover, we already know that the volume emerging from a single mode would give rise to an effective Friedmann equation of the form of (28) or (31), which corresponds to an isotropic Universe. These two statements agree with the classical observation that the Bianchi I relational dynamics (49) do not differ from FLRW (40) if the anisotropy parameters are constant, as they only appear through derivatives. Another way of seeing this is to observe that a Bianchi I metric for which the Misner variables are constant can be brought to FLRW form by rescaling coordinates. Hence, to make β ± (χ) dynamical one needs to allow for contributions coming from multiple shapes, i.e., modes with different values for β j,ı ± . In order to compute the sums (55), we need to specify the single-mode expectation values N J (χ) and the coefficients V j,ı and β j,ı ± . The volume eigenvalues V j,ı , as introduced in (13) and (22), are imported from LQG; a procedure to (numerically 7 ) compute them is well-known [50], and will be recalled in appendix A. The expectation values N J (χ) were discussed in section II, and we have explicit expressions for the two approaches given in (27) and (29), using that for a single mode ± , on the other hand, deserve some further elucidation.
A. Defining β j,ı ± : the trisohedral tetrahedron Unlike for quantities such as areas and volumes, there is no fundamental operator in LQG corresponding to Misner variables β ± . Our task is therefore to give a proposal for β j,ı ± as functions of such more fundamental geometrical quantities, defined at the level of each tetrahedron. These variables (areas and volume), in turn, are determined by the spins j and intertwiner ı. We think of the fundamental tetrahedra as embedded into a manifold with Bianchi I metric, and reconstruct parts of that metric from geometric quantities of a tetrahedron, as originally proposed for GFT in [9]. There is clearly some ambiguity in any such procedure, given that the required embedding is not part of the GFT formalism but additional input. Our proposal is a relatively direct extension of previous work in GFT cosmology: we follow the idea [10] that equal spins for a four-valent node (i.e., j = (j, j, j, j)) capture a discrete notion of isotropy; hence such a configuration should correspond to β ± = 0. Departures from this microscopic notion of isotropy are, in a sense, assumed to add up coherently to a macroscopic definition of anisotropy.
For any excitation associated to J = ( j, m, ı), the spins j determine the areas of faces of a tetrahedron (see figure 1) 8 . However, these areas are not sufficient to determine the shape of a tetrahedron; there is a two-sphere's worth of different tetrahedra with given face areas [50], so additional assumptions are needed to identify a given quantum state with a configuration of a 7 Given that these eigenvalues can only be computed numerically there is little or no hope to be able to solve the sums (55) in a closed form. We will have to truncate them and make some suitable choices for the modes. 8 This is again a statement imported from LQG, where the eigenvalues of the area operator AI associated to the I-th face of a quantum tetrahedron are given by l 2 0 jI (jI + 1), l0 being a fundamental length scale parameter. classical tetrahedron. The idea of [10] is that for equal spins we should choose a regular tetrahedron, whose edges are all of equal length 9 . To justify this assumption one fixes the intertwiner ı to maximise the volume eigenvalue, as this would represent a situation which goes as close as possible to the classical picture. More generally, we decide to focus on orthocentric [51] simplices (discussed in some detail in appendix A). Orthocentric tetrahedra maximise the volume for given face areas. The regular tetrahedron is a particular example of orthocentric tetrahedra. This motivates us to always choose the largest allowed volume eigenvalue for given face areas, and hence spins j; we will then interpret our states as orthocentric tetrahedra. For a few specific shapes, we explicitly show in appendix A that the largest volume eigenvalue for a given mode is indeed the closest one to the classical volume of an orthocentric tetrahedron for the same face areas. In this sense, a choice for ı is dictated by the comparison that we make between GFT models and classical geometry.
The sum over J then reduces to spins j and magnetic indices m only, since for each choice of spins the intertwiner ı is already fixed in all the sums from now on. Symbolically, we are left with since we always use the largest volume eigenvalues in such sums. Note that we still need to make sure that each combination of j I included in the sum allows for a nonvanishing intertwiner. A very minimal requirement for our definition of β j ± (now implicitly associated to the intertwiner fixed by the j) is that it should vanish for j = (j, j, j, j), but be nonzero if at least one of the four spins differs from the others. The simplest "non-isotropic" (but still orthocentric) building block one can think of is a tetrahedron we will refer to as trisohedral. As one can see in figure 2, this too is a quite particular shape, with three isosceles triangles (called "sides" with area A) and an equilateral one (called "base" of the tetrahedron, with area B).

FIG. 2:
The trisohedral tetrahedron has two types of areas, edges, or dihedral angles (between sides and between a side and the base), et cetera. It generalises the regular tetrahedron, still remaining rather specific. We propose to associate a notion of anisotropy to such a non-equilateral shape.
Following the analogy from before, we then make the assumption that such a tetrahedron is 9 Such a platonic solid is specified by a single number (e.g., edge length, height, distance between opposite edges) with fixed dihedral angles between faces, et cetera.
represented by modes of the form where the spin j a is associated with the area of the sides and j b with the base. Again, this assumption is partially justified by choosing an intertwiner that maximises the volume eigenvalue. It should be clear, however, that for given face areas this volume eigenvalue will not exactly match the classical volume of a trisohedral tetrahedron. We again refer to [50] for more discussion of this in the context of LQG, and give more details in appendix A. In order to specify the numbers β (ja,ja,ja,j b ) ± , we now think of the trisohedral tetrahedron of figure 2 as embedded in a locally rotationally symmetric (LRS) Bianchi I spatial slice. This is a Bianchi I geometry with only one preferred direction in which the expansion (or contraction) differs from the other two directions. The (spatial) line element follows from (41) and (42) Because of the symmetry of our building block our model only needs two different scale factors, so we have set β − = 0. Now consider a tetrahedron embedded in such a space, with one of its triangles lying on the x − y plane. The tetrahedron is chosen such that it would be regular with respect to the background "fiducial" coordinates x, y and z; but its physical geometry depends on the dynamical variables 10 . In particular, the areas of the equilateral base and the isosceles sides are We identify B and A with the area eigenvalues associated to the spins j b and j a (see figure 2). We then see that a configuration in which B and A are different would be interpreted as anisotropy (β + = 0) whereas B = A corresponds to β + = 0, as anticipated. Inverting either one of equations (59), one can express the anisotropy parameter as a function of face areas and volume, and define this to be the "anisotropy" associated to the tetrahedron. In the quantum theory, A, B and V are represented as possible eigenvalues determined by the spins. This then leads us to possible proposals for how to define β (ja,ja,ja,j b ) + in (55). In appendix A we compare different definitions for β (ja,ja,ja,j b ) + obtained from (59). They do not agree, because we fixed the volume of the tetrahedron to agree with the quantity V in the metric assuming that we have an orthocentric tetrahedron, but there is no quantum eigenvalue satisfying exactly the same relations between volume and areas. Only one definition is simple and satisfies our desired property β (j,j,j,j) + = 0; as one might have anticipated, this is a definition that does not use the volume eigenvalue at all. Indeed, from the ratio A/B in (59), one finds β + = −1/6 log (9A 2 − B 2 )/(8B 2 ) . Converting to LQG eigenvalues, we then define an effective local anisotropy associated to the mode (57) of the quantum tetrahedron by This gives zero if j a = j b and if we assume j ≥ 1 2 , (60) is always finite and well-defined, regardless of the volume eigenvalue for the given mode (this is not true for other definitions, see appendix 10 We fix the edge length l = √ 2 3 √ 3 such that the physical volume of the tetrahedron is equal to V in (58).
A) which may vanish for some spin configurations. For instance, V (1,1,1,3) = 0. Note that the anisotropy is conventionally negative when j a > j b . With this definition, we are finally in a position to evaluate (55). We now focus on the initial conditions needed for tackling the sums.

B. Initial conditions
Independently of the approach one follows, for a mode specified by (57) the kinetic term (25) is characterised by Since we would like to include a number of Peter-Weyl modes in (55) for which m 2 j > 0, the coupling M needs to be small relative to the "mass" m. For any given ratio M/m there will be maximum spins after which (61) becomes negative; for instance, choosing M/m = 0.1 means that 3j a (j a + 1) + j b (j b + 1) < 100 for m 2 j > 0. In the numerical analysis below, we will assume M/m = 0.1.
Given that we are not interested in the bounce or the connection between the contracting and expanding phases, we will be looking at χ-symmetric solutions, for which the minimum of the volume is at χ = 0. Thus, in the general expressions in section II, we set α + J = α − J ≡ α J in the algebraic approach (cf. (26)- (27)) and K(0) = 0 in the deparametrised formulation (cf. (29)). The expectation values for single-mode number operators then become N MF j, m (χ) = |σ j, m (χ)| 2 = 2|α j, m | 2 1 + cosh(2m j χ) = 4|α j, m | 2 cosh 2 m j χ for the algebraic method with a mean-field (MF) approximation, and for the deparametrised (DP) approach. In (62), α j, m is related to the number of quanta at χ = 0 as N MF j, m (0) = 4|α j, m | 2 . In the deparametrised approach, we then fix N DP j, m (0) = 4|α j, m | 2 so that the initial conditions are the same for both methods.
We then need to fix the range of the sums in (56). We are interested in modes of the form (57), specified by two spins j a and j b . In principle, the sums over j a and j b run from 1 2 to ∞; but in practice, they have to be truncated because the eigenvalues are only computed numerically. Furthermore, not all combinations are allowed by SU (2) recoupling theory, as not all of them allow for a nonvanishing intertwiner. Moreover, modes corresponding to zero volume are not really physically interesting and would not allow for a useful notion of anisotropy. In principle, all these considerations are to be taken into account. We will simplify these issues greatly by only keeping a few modes.
We also need to sum over magnetic indices m. None of the geometric observables we are considering depend on m, so the m index corresponds to a degeneracy factor for physically indistinguishable modes. Given this, we will assume that the initial condition parameters N j, m (0) do not depend on m, and so the coefficients α j, m ≡ α j in (62) and (63) are independent of m. The sums over m then just return additional multiplicative factors, We then still need to truncate these sums to a few chosen values for j a and j b . One motivation for keeping only a small number of modes (but more than one) is the increasing arbitrariness coming from the need to specify initial conditions for each additional mode. This question of dependence on initial conditions generally affects cosmological models derived from fundamental quantum gravity, which tend to depend on some specific choice of (class of) initial states [33,34].
There is then a choice of which modes to include, in particular whether j a and j b should be large or small. Here two main factors may come into play. At the kinematical level, the quantum properties of tetrahedra show semiclassical geometrical features for large spins [40,50,52]; large spins allow for more possible discrete (quantum) states, so that classical concepts such as orthocentric tetrahedra can be approximated closely. Indeed, the largest volume eigenvalues we choose are always smaller than the volume of the corresponding orthocentric tetrahedron, but the relative difference decreases as the spins grow. We discuss this feature in appendix A for two specific examples. It would suggest keeping modes associated to large values of (j a , j b ), since the geometric picture of orthocentric tetrahedra would be closer to the volume eigenvalues used.
On the other hand, we know that for the type of GFT dynamics considered in this paper, the lowest spins will dominate the sum eventually [25], and larger spins give less important contributions to the dynamics at late times. This is why in previous work the sum was often trivialised to only the smallest spins, or to a few values. This truncation is analogous to loop quantum cosmology where the spins are usually all fixed to be j = 1 2 (see [53] for an attempt at generalisation). Therefore, even though using large spins might make sense kinematically, when we include the dynamics we see that larger spins will quickly become insignificant for the evolution of the cosmological model. Finally, on general grounds one would expect that, for any fundamentally discrete (simplicial) approach to quantum gravity, a useful continuum limit is obtained for large simplicial complexes made of little simplices, rather than by magnified semiclassical building blocks. Thus, we will only consider a few relatively small spins.
We now need to specify initial conditions given by the initial number of quanta N ja,j b (0). These initial conditions are of course in general arbitrary, but for concreteness we assume they follow a Gaussian distribution where we always set N ja,j b = 0 for modes with vanishing volume. For continuous parameters, such distributions allow for analytical integrations, as we will show in a toy model later. Here the spins are discrete, so we are considering a "discrete Gaussian distribution", assuming that the initial configuration is dominated by values around (j a , j b ). By setting these away from j a = j b = 1 2 , we will see some initial contribution coming from larger spins, before the lowest ones take over dynamically. Moreover, given that we effectively set to zero all terms after an arbitrary spin, it is reasonable to choose initial quanta having occupancy numbers that gradually go to zero, as modes approach the last one. This motivates using a Gaussian also in the discrete case. All N ja,j b (0) are now determined by the peaks j a and j b , standard deviation σ and normalisation factor N .

C. Three modes with fixed base
We now focus on a simple model in which we assume the spin associated to the base of the tetrahedron to be fixed to its minimum value j b = 1 2 . We keep three modes associated to the lowest (allowed) three spin values for the sides of the tetrahedron, j a = { 1 2 , 3 2 , 5 2 }. This is a simple generalisation of the single-mode case, since we consider two modes which encode anisotropy (as defined in section IV A) and one associated with equilateral tetrahedra (see figure 3). Note that the first mode corresponds to the regular tetrahedron.
, and noticing that (2j b + 1) = 2 in this case, we truncate (55)  can be found in table I at the end of the paper. Our initial condition function (65) now only depends on j a . We fix the peak to be at j a = 3 2 ; σ is chosen to be around 1 so that the initial distribution is peaked at j a = 3 2 but also includes some quanta in the other modes. N is fixed by the requirement that the initial state should be reasonably semiclassical; given that we only focus on expectation values of operators, quantum fluctuations should not be too large (see [13] for a detailed analysis of relative fluctuations). As these fluctuations decrease with N 1 (this is a generally expected behaviour also found, e.g., in [54]), we demand that initially all modes have an expected particle number of at least 25. Figure  4 shows different σ choices according to these criteria.
In figure 5, we plot the effective Friedmann equation [V (χ)/V (χ)] 2 and the anisotropy evolution β + (χ) stemming from (66) and (67). We compare the effective Friedmann equation with the classical (relational) Bianchi dynamics reviewed in section III (namely, solutions of (49) and (54)), where we set β − = 0 because of the local rotational symmetry of our GFT model. While the Bianchi I model is a natural comparison, we also compare with Bianchi II as this is a less simple classical anisotropic cosmology, which deviates from Bianchi I only at late times (which is what the GFT description of the anisotropy does too, albeit in a different way). The other free parameters of the classical plots are fixed as follows: Newton's constant G is determined by demanding that at late times we follow the isotropic solution [V (χ)/V (χ)] 2 = 12πG; the slope dβ + /dχ is obtained from a linear fit of the initial part of the GFT anisotropy dynamics.
The value of p χ could be fixed from the fundamental GFT by using the quantity playing the role of a conjugate momentum of the scalar χ. In the algebraic approach, (28); in the deparametrised formulation the conjugate momentum of χ is given by the (relational) Hamiltonian p χ = H ≡ J Ĥ J as explained below (31). For our states with α + J = α − J , the quantities Q J are actually zero, but one could introduce an arbitrary phase into α + J or α − J to obtain a nonvanishing p χ . In the deparametrised approach, H = j m j N j (0) if we assume coherent states (30) and real α J . In either case this gives a relatively small p χ and in the plots Bianchi II would deviate from the linear Bianchi I behaviour very quickly, which is not what we see in the GFT model. In an attempt at obtaining a better fit, in figure 5 we fix p χ by assuming that the nonlinear behaviour of Bianchi II appears roughly when the function β + (χ) in GFT ceases to be linear.
Notice that the constant dashed red curve in the left panel represents Bianchi I but is indistinguishable from the asymptotic FLRW limit because the classical anisotropy backreaction (cf. (49)) would be very small: it scales as (dβ + /dχ) 2 ∼ 4 × 10 −5 if we take dβ + /dχ as the slope of the initially linear part of the GFT expression for β + (χ). Figure 5 shows the dynamical "isotropisation" already mentioned in the literature [17]. This late-time limit is inevitable given that, for a model involving multiple Peter-Weyl modes, the mode with the largest value that m j can take will end up dominating the dynamics [25]. In our case this is achieved by the smallest spins j a = j b = 1 2 ; this mode dominates at some point and will then dominate indefinitely (see figure 6). When this happens, the effective Friedmann equation reaches a constant plateau (described by the first term in (66)) while the anisotropy gradually converges to zero (the crossed-out term in (67) now dominates the average). The constant value taken by (V /V ) 2 in this asymptotic limit is then the one we would have obtained for a single-mode model corresponding to the FLRW universe. In the left-hand plot of figure 5 we label this constant value as "FLRW/Bianchi I" since, as we explained, the difference between FLRW and Bianchi I is too small to be noticeable in the plot.
The period before this asymptotic isotropisation (but after the bounce) can be compared with the classical Bianchi I model. In the left panel of figure 5 we see that the GFT curve does not show a constant [V (χ)/V (χ)] 2 (with value higher than the asymptotic FLRW value), as the classical dynamics (49) would suggest. Instead, the green line shows the transition between different modes, which will eventually stop when the last mode takes over and dominates. This behaviour can be made more or less evident by changing the standard deviation σ in (65). We will see an example with a manifest transition between modes later on.
In the right panel of figure 5 we see that β(χ), on the other hand, shows a much better agreement with classical (relational) cosmology. The anisotropy is approximately linear for a reasonably long time after the bounce, as would be expected from general relativity. Because of the tendency to isotropise, this agreement will obviously stop at some point as the anisotropy will approach zero; but before that point, β(χ) compares well with a classical Bianchi I model. It is worth reiterating here that in this construction we are ignoring GFT interactions (cf. section II C), which become important when χ becomes sufficiently large 11 . Hence these plots are not to be trusted for too late times, when the weak-coupling assumption breaks down. In particular, we might never reach isotropisation but simply remain in a phase where β(χ) is approximatively linear, before interactions take over and change the picture completely.
Regardless of when β + (χ) ceases to be linear, it is always monotonic; so in principle, as in 11 When the total particle number is large, one expects correlations between the GFT quanta to be non-negligible. classical cosmology, it could be used as a relational clock for GFT cosmology instead of the ad hoc introduced scalar field χ. Adding matter arbitrarily is often seen as an inescapable necessity plaguing any cosmological model coming from a background-independent quantum gravity theory, given the absence of (time) coordinates. In contrast, this model contains a gravitational degree of freedom which could be used as relational time, and perhaps help addressing issues such as the problem of time or clock dependence in quantum cosmology (for more on this see, e.g., [14]). So far, we have completely glossed over the question whether the time evolution of the particle number in each mode follows (62) in the algebraic (mean-field) approach, or (63) in the deparametrised approach. It turns out that, as one might have expected, the two approaches give identical results after a very short initial phase directly after the bounce, see figure 7. for both approaches and with the same parameters (fixed as in figure 5).
The discrepancy between the two approaches can be traced back to the different effective Friedmann equations for a single mode, as discussed in section II C (see (28) and (31)). After the 1/V (χ) 2 contribution dominates (giving rise to the bounce), and before the constant term of the Friedmann equation takes over, there is a difference in the 1/V (χ) term which in the mean-field approach depends on the arbitrary "GFT energy" E J (which is negative for (65) and our choice of parameters), whereas in the deparametrised approach it is fixed and positive. The plot for β + (χ) also shows a minor difference between the two approaches. Changing initial conditions changes the details of these differences, but the underlying features are always the same; and given that we are interested in comparing later times (larger volumes) with classical models, these differences do not play a role in our analysis. Already for these small values of χ we observe an almost-constant [V (χ)/V (χ)] 2 and quasilinear β + (χ).

D. Mean-field continuous toy model
As discussed, the general expressions (55) must be evaluated numerically, even in simple cases such as the model of section IV C. In order to obtain some analytical result, one can consider toy models which allow explicit solution due to some simplifying assumptions. The model presented in this section uses the solutions in the mean-field approximation of the algebraic approach (62), but drops the assumption of discrete spins j.
Recall that the classical relations (59) relate the two types of area (for the base B and the sides A of the trisohedral tetrahedron) with the volume and the anisotropy parameter, Hence, classically the pair (V , β + ) is in one-to-one correspondence with the face areas (A, B), and it would be easier to take these variables as the basic characterisation of a tetrahedron, expressing A and B as functions of them. We will do this here, and forget about the fact that in LQG (and GFT) the fundamental variables are discrete areas A, B.
If we focus on V and β + , we can make another simplifying assumption: we assume that the volume per tetrahedron is simply fixed to be V = V 0 , and that only β + varies continuously. We follow the picture coming from single-mode truncations of GFT that the evolution of the total volume comes only from adding or removing building blocks, rather than changing their "size"; but we still have a range of different "shapes". This avoids the previously inevitable mixing of effects coming from modes with both different volume eigenvalues and different values of β + .
We want to work with the same GFT kinetic term, i.e., the same expression (61) for the effective couplings for each mode. But given that we no longer have discrete spins, we need to rewrite this expression using V and β + . We will do this by using the LQG relation A 2 I = l 4 0 j I (j I + 1) for area eigenvalues (see footnote 8) and the classical relations (59). We can then write down a relation between the discrete and toy models for modes of the form (57), The symbol= means that we are equating quantum eigenvalues with classical quantities using the LQG interpretation given to the spins. This relation is then substituted into the time-dependent expression of the average particle number (62), with initial conditions fixed by a suitable choice of α which we again take to be a Gaussian, α β + = exp{−(β + − β + ) 2 /(2σ 2 )}. Given that we do not have discrete modes, we do not have to worry about a specific normalisation factor. This is now a continuous normal distribution, the analogue of (65). The "toy model counterparts" of (55) then read and Numerical evaluation of [V (χ)/V (χ)] 2 and β + (χ) shows qualitatively similar behaviour to the plots in figure 5 derived in the previous three-mode model; see figure 8 for an example. Again we see an asymptotic "isotropisation" leading to the effective FLRW limit at late times, and an approximately linear evolution for the parameter β + . We can now also strengthen these results by analytical approximations.
In order to find analytical solutions to the integrals (70) and (71) we assume small anisotropy, i.e., a Gaussian distribution with β + 1. This justifies approximating the cosh function under the integral up to second order in β + , where m := m 2 − 3M 2 (3V 0 /l 3 0 ) 4/3 . This approximation leads to integrals over β + that can be done immediately, and relatively simple expressions such as To simplify these further we then use the fact that M/m is small (see discussion below (61)) and expand β + (χ) and [V (χ)/V (χ)] 2 up to second order in M/m. This yields a simple analytical expression for the anisotropy, and the effective Friedmann equation where F(mχ) := sech 2 (mχ) (2mχ + sinh(2mχ)) tanh(mχ). We can see that the anisotropy (74) is essentially a linear function of χ very soon after the bounce, with a slope dependent on parameters of our Gaussian such as β + and σ. Moreover, noticing that F(mχ) |χ|→∞ → 2, we also obtain a constant late-time limit for the effective Friedmann equation, By comparing with exact numerical results we can see that these approximations cannot be trusted for too large χ, see figure 8. In the limit in which we "switch off" anisotropic contributions, β + → 0 and σ → 0, (74) vanishes and in (76) we have (V /V ) 2 → 4m 2 − 12M 2 (3V 0 /l 3 0 ) 4/3 which is nothing but the orange line in figure 8. In (69) this corresponds to 4j(j + 1) = 3(3V 0 /l 3 0 ) 4/3 when the spins are all equal.
As a final comment, recall that in the classical Bianchi I model the Friedmann equation gets a constant contribution (dβ + /dχ) 2 compared to the FLRW Universe. We can ask whether the terms in (76) that depend on β + and σ are related to the derivative of (74). But we see that As already noticed in the full GFT model, contrary to what happens in general relativity, in our toy model the presence of anisotropy decreases (V /V ) 2 . The two quantities we are comparing also are of different orders of magnitude, in particular different powers of the small ratio M/m. Our toy model could reproduce the main qualitative features seen in the full GFT analysis, in particular a nearly linear growth in the anisotropy for a range of χ and a negative contribution to the effective Friedmann equation.

E. Including more modes into GFT models
We now return to the setting of GFT, presenting results for models that go beyond the simple case described in section IV C. The easiest extension of what we showed before is to include more than three modes, but keep the assumption of a fixed base area (i.e., j b = 1 2 ). This means we add shapes to figure 3 for greater j a which are increasingly more stretched ("anisotropic"). We find that such an extension gives results that are not qualitatively different from the previous case, regardless of how many modes of this form we add. In figure 9 we show the case of five modes, letting j a range between 1 2 and 9 2 . shifts between modes before reaching the plateau for j a = j b = 1 2 , and β + (χ) can be approximated by a linear function for some time before going to zero. Initial conditions and parameter choices are shown in figure 10.
Since j b = 1 2 for all modes, we can see from (60) that the values of β + for the modes we consider all have the same sign, and are increasingly more negative as j a grows. This is why β + (χ) does not lose its monotonicity property in figure 9. From the geometrical point of view, the shapes we  (65)) for five and for eleven modes. As in figure 4, N is such that the number of quanta is 25 in the modes which are the farthest from the peak. Left: j b = 1 2 , j a ∈ { 1 2 , . . . , 9 2 }, j a = 5 2 and σ = 1. Right: (j a , j b ) ∈ { 1 2 , . . . , 5 2 }, j a = j b = 3 2 and σ = 1 2 . Combinations of spins with zero volume are not included.
are adding are progressively more stretched along one axis so that their local anisotropy never flips the direction in which it changes.
A more important generalisation of our model can be obtained by relaxing the assumption that j b is fixed to 1 2 . If j b can vary, the first minor novelty comes from the fact that we now have some combinations of spins which need to be removed; these are spin configurations for which volume eigenvalues are zero. For instance, when j b = 3j a the volume eigenvalue vanishes (see table I and appendix A) as one might expect from classical arguments (the tetrahedron of figure 2 flattens into a plane when B = 3A). A second key novelty comes from the fact that not all the β ja,j b + parameters have the same sign: in the dynamical evolution, we now no longer obtain a strictly monotonically increasing sequence of β ja,j b + values associated to the modes dominating at different times. Hence, we generically find a non-monotonic β + (χ). We show in figure 11 an example with eleven generic modes, defined by letting both spins vary in the range { 1 2 , . . . , 5 2 } and excluding the ones not allowed by SU (2) recoupling theory. See right panel of figure 10 for a depiction of the initial number of quanta in the allowed modes. One could change initial conditions such that β + (χ) is monotonic even with variable base spin j b , by choosing ad hoc modes such that the succession of dominant β ja,j b + values goes to zero monotonically. While any number of modes can be included without any computational obstacle, we do not report additional details on these many-mode scenarios because they are characterised by a larger arbitrariness encoded in further initial conditions, without introducing important novelties.

V. CONCLUSIONS
This paper is a first attempt at characterising homogeneous but anisotropic (Bianchi) cosmologies in group field theory (GFT), focusing on the simplest case given by a Bianchi I model with an additional local rotational symmetry. Previous work on GFT cosmology had included anisotropies perturbatively, showing that for a certain class of kinetic terms they undergo a "decay process" leading to isotropisation. However, what was largely missing was a precise definition that would allow to discern whether an effective cosmology in GFT is isotropic or not, and to quantify the amount of anisotropy. Our work fills this gap in the literature by proposing one particular measure of anisotropy and by studying its dynamics. This is done within the context of simplified models which share the main setup and basic premises (such as neglecting interactions and making assumptions on the form of kinetic term) of previous work in GFT cosmology.
Inspired by Misner's parametrisation of the Bianchi I model, which introduces anisotropy parameters β ± behaving as free massless scalar fields on a flat FLRW background, we defined an analogue of the β ± variables in GFT. Given that GFT inherits from loop quantum gravity notions of geometry which are discrete, the construction of analogue anisotropies requires some care. Our proposal is to define a quantity which to some extent follows the structure of expectation values of GFT operators, assigning a microscopic value to anisotropy at the level of each Peter-Weyl mode which is multiplied by the number operator for this mode, but then dividing by the total number of quanta in order to obtain an "intensive" quantity. The question is then what function of LQG eigenvalues for areas and volumes should be used to effectively describe the anisotropy of a quantum building block. This function is analogous to eigenvalues of LQG operators used, e.g., to define a total volume in the usual treatment.
Having discussed various options for such an effective notion of anisotropy, and defined one which satisfies the demand that a tetrahedron with equal face areas is considered isotropic, we then studied the dynamics of the expectation value of the total volume and of the newly introduced anisotropy observable (due to the additional rotational symmetry, there is only one anisotropy degree of freedom). We were mostly interested in a relative late-time regime far away from the GFT bounce, where one would hope to see reasonable semiclassical physics emerge. At very late times, previous results already suggested a process of isotropisation and reduction to effective FLRW cosmology, so the regime we are interested here is before this final stage.
For isotropic cosmology, restricting the analysis to a single Peter-Weyl mode is enough to obtain a simple cosmology that reduces to general relativity at large volume while resolving the singularity by a bounce. However, such a restriction is too drastic to allow for dynamically evolving anisotropies. We hence considered various simple models which include multiple modes.
One model includes three Peter-Weyl modes, one describing equilateral (isotropic) and the other two describing trisohedral (anisotropic) tetrahedra. Here we found partial agreement with general relativity: while the evolution of the volume does not show the anisotropy backreaction that one would expect classically, the dynamics of anisotropy match the linear evolution dictated by classical relativity quite well. The behaviour one would expect for the volume classically is a faster rate of expansion when anisotropies are present compared to the isotropic case; but this can never be matched in the type of GFT model we consider, with a certain kind of kinetic term, since the fastest rate of expansion is reached for equilateral (isotropic) tetrahedra. This is precisely also why asymptotically only equilateral tetrahedra dominate and we see isotropisation at very late times. Hence, something more drastic than simply changing initial conditions or, e.g., considering more modes would be needed to find agreement with classical relativity. For instance, one could use different types of kinetic term, or include the effects of GFT interactions. Interactions will always dominate at some point as the Universe grows to a certain volume, and so the asymptotic isotropisation seen in a noninteracting approximation may not be physically relevant. Despite all this, the anisotropy observable of our model shows exactly the classically expected behaviour for some period of time before isotropisation, which is promising.
We showed that one can include more modes into the model, and depending on the additional assumptions one makes this can lead to slightly different results. In particular, if one of the four spins (associated to the base of the tetrahedron) is fixed, including many modes gives qualitatively similar results; but if one relaxes that assumption, one can in general obtain a non-monotonic evolution for the anisotropy.
Our main results could only be obtained numerically, and the origin of the exactly linear evolution of the anisotropy was not transparent. To gain further insights we also studied an analytical toy model with further simplifications; assuming that anisotropy is a continuous parameter while the volume per tetrahedron is fixed, we could find expressions for the effective Friedmann equation and anisotropy dynamics which reproduce the main features of the discrete GFT models. Here we found that the toy model expression for β + shows a linear behaviour soon after the cosmological bounce, as required by general relativity and as seen in the more refined numerical analysis.
Let us mention possible developments stemming from this work. As already pointed out, one could ask whether the effective dynamics of anisotropies change if GFT interactions are taken into account. This would complicate the resulting calculations further, beyond the need for multiple Peter-Weyl modes (see, e.g., [13] for an already quite involved numerical study of a single mode with interactions). A different direction would be to analyse the detailed effects of anisotropy on the bounce phase. In this work we focused on a post-bounce regime because we aimed to compare with classical relativity, but in general one might ask whether the bounce itself could be spoiled (or in general modified) by anisotropies, as is often a main worry in bounce scenarios in which anisotropies dominate asymptotically on approach to the singularity. In our model, we have a massless scalar field as matter, which would classically prevent the domination of anisotropies. We saw that anisotropies are present at early times and disappear at late times, but singularity avoidance does not seem affected by the inclusion of anisotropies. The details of the bounce may still be altered by their presence. Another line of investigation, closely related to the previous one, would be to compare the exact details of the cosmological bounce with the similar singularity resolution of loop quantum cosmology. In the Bianchi I context, this would mean to investigate whether the anisotropic nature of the model has the same influence (if any) on the bounce in loop quantum cosmology [19] and in our work.
One restriction of this paper was that we were studying the GFT analogue of a Bianchi I Universe with an additional rotational symmetry, so that there is only one β variable rather than two. To lift this restriction, one could study more general types of tetrahedra rather than trisohedral ones, discuss different proposals for β ± variables in this more general context, and study their dynamics along the lines we have discussed. We would anticipate additional ambiguities in such a process on top of the ones mentioned in the appendix, and perhaps not many additional insights beyond the ones here, given that both β ± variables appear on the same footing in the Bianchi I model.
Finally, given the monotonic evolution of the new anisotropy observable β + in GFT, one may hope to use such a gravitational degree of freedom as relational clock. This would be similar to what happens at the classical level, where in the simplest Bianchi I Universe without matter the classical Friedmann equation can be written as so that β + can be a relational clock with no need for a separate matter field. One might hope to incorporate this idea into GFT, and describe relational evolution without coupling to the somewhat arbitrary massless scalar field. Such a relational formalism would appear to require using an expectation value as a clock parameter, perhaps along the lines suggested in [14]. The availability of multiple candidate clock degrees of freedom would also help identifying a GFT analogue of the symmetry of continuum general relativity under time reparametrisations. We leave these open questions for future work.
for the regular and the trisohedral tetrahedron of figure 2, respectively. In the quantum theory, we take the analogue of an orthocentric tetrahedron to be the SU (2) intertwiner with largest volume eigenvalue; we explain the concept of LQG volume eigenvalues in what follows. Following the terminology introduced in section II, we focus a single four-valent node, with spins labelling its links denoted by j = (j 1 , j 2 , j 3 , j 4 ) (see figure 1). To each representation j I , with I = 1, . . . , 4, we associate a vector space H j I that carries the action of the SU (2) generators J I . The Hilbert space of the quantum tetrahedron then reads H 4 = Inv [H j 1 ⊗ . . . ⊗ H j 4 ], and objects that live in this space are called intertwiners. A nonvanishing intertwiner can only exist if the j I sum to an integer. We introduce a basis labelled by k in the recoupling channel H j 1 ⊗ H j 2 , where the index k ranges between k min = max {|j 1 − j 2 |, |j 3 − j 4 |} and k max = min {j 1 + j 2 , j 3 + j 4 } in integer steps. The Hilbert space H 4 is d-dimensional with d = k max − k min + 1. States on this space can be understood as describing quantum tetrahedra, as firstly pointed out in [39,41].
We refer to the literature (see, e.g., [55]) for derivations of the geometrical eigenvalues in the context of LQG. Here we only recall that the area operator takes the form A = l 2 0 J · J, where l 0 is a fundamental quantum gravity length scale, usually taken to depend on the Barbero-Immirzi parameter of LQG. States in H 4 are eigenstates of the operator A I , which measures the area of the I-th face of the quantum tetrahedron, with eigenvalues l 2 0 j I (j I + 1). Finally, the volume operator introduced for a tetrahedron in LQG reads [56,57] We can focus on the radicand in (A4). Without showing all the details of the calculation (which can be found in [50,52]), one defines the operator Q := J 1 · ( J 2 × J 3 ) with matrix element Q k,k ≡ k| Q|k , where k and k label intertwiner states as explained above. One can then show that nonvanishing matrix elements are obtained only if k and k differ by 1. Thus one can denote a k := iQ k,k−1 , so that is a d × d matrix. With the above conventions, the matrix elements are found to be [58] a k = 1 4 (j 1 + j 2 + k + 1)(−j 1 + j 2 + k)(j 1 − j 2 + k)(j 1 + j 2 − k + 1) √ 2k + 1 × (j 3 + j 4 + k + 1)(−j 3 + j 4 + k)(j 3 − j 4 + k)(j 3 + j 4 − k + 1) √ 2k − 1 .

(A6)
To obtain a compact expression, (A6) can be cast in terms of Heron's formula (A2) as Hence, for given spins j computing the spectrum of the volume operator amounts to finding the d eigenvalues of (A5) (let us denote them q j,k ) and then computing, according to (A4), Note that if j = (j, j, j, j) the matrix elements (A7) simplify to a (j) while for j = (j a , j a , j a , j b ) (as in (57)) they reduce to These matrix elements can be used to find the maximal volume eigenvalues V max j and V max ja,j b for the quantum shapes corresponding to regular and trisohedral tetrahedra, respectively; see table I  for examples. For large spins, these volume eigenvalues show semiclassical behaviour (see, e.g., [40,50,52]), and so the "eigenvalue-counterparts" of equations (A3), (A11) become more accurate as the spins grow. Here V max means that we fix the intertwiner k so as to obtain the largest volume eigenvalue.
In figures 12 and 13 we show the possible volume eigenvalues V k j and V k ja,j b in units of l 3 0 (dots), compared with the maximal classically allowed volume (curve and surface) for the same face areas (given by l 2 0 j(j + 1)). We also show how the relative difference between this classical volume (for orthocentric tetrahedra) and the highest LQG volume eigenvalue decreases as the spins increase. Figure 12 focuses on the case of equal areas j a = j b , whereas figure 13 shows various other combinations of j a and j b . Notice that both the classical and quantum volume of a tetrahedron are no longer well-defined for j b > 3j a . In the quantum theory, this constraint comes about because the dimension d of the Hilbert space H 4 needs to be positive. The right panel of figure 13 shows, as an illustrative example, the relative distances between the maximum eigenvalues and the surface along the plane j a = 2j b , but the qualitative behaviour is similar for other ratios.
dots) and the classical volume of an equilateral tetrahedron as function of the area (curve). Right: relative difference between the classical volume and largest LQG eigenvalue. After two initial anomalies the mismatch decreases indefinitely: it becomes reasonably small (∼ 1%) when j ∼ 40, and goes down to ∼ 0.1% when j ∼ 350. ja,j b /l 3 0 (dots) and the classical volume of a trisohedral tetrahedron (transparent surface) as function of the two areas. We plot eigenvalues along specific planes for which j a = j b (blue), j a = 2j b (orange), j a = 5j b (green), 2j a = j b (red) and 3j a = j b (grey). The latter gives vanishing eigenvalues only and represents the limiting case as there are no tetrahedra (no nonzero intertwiners) for 3j a < j b . Right: relative discrepancy along the plane j a = 2j b .
We now finally turn to the question of how to define a notion of anisotropy at the level of quantum geometry, in terms of a function β ja,j b + . In particular, we show that though the definition (60) is not unique, it is the best candidate to play such a role.
Starting from (59), there are three ways to classically define the anisotropy of a tetrahedron in terms of two other geometric quantities: by inverting any one of (59) or combining them, one can obtain classical expressions such as β + (V, B), β + (V, A) or β + (A, B). Once we pick a favoured definition, the strategy is to replace classical expressions by LQG eigenvalues to obtain an effective (discrete) notion of anisotropy as a function of the spins associated to the faces of a tetrahedron.
For instance, inverting the first equation in (59), one can readily obtain where we simply replaced the classical area B with l 2 0 j b (j b + 1) and the classical volume with the maximal eigenvalue V max ja,j b . Notice that by definition β ja,j b + > 0 when j a < j b (naively, for volumes smaller than the equilateral ones). (A12) could be used as a definition, but we can identify a number of issues. First, it cannot be written more explicitly, since the volume eigenvalue V max ja,j b is only obtained numerically, as outlined in the previous section. Moreover, (A12) is not well-defined when the volume is exactly zero, and it is nonvanishing when j a = j b .
The same arguments would apply if we defined the anisotropy starting from the second relation in (59). In this case we would obtain an even more cumbersome expression β + (V, A), which could then be turned into a β ja,j b + as a function of j a and volume eigenvalues V k ja,j b . A third definition is obtained by taking the ratio of equations (59). The advantage of this choice is that one can get rid of the volume: We can now rearrange to obtain a simpler definition of the anisotropy parameter as a function of A and B, which straightforwardly turn into the spins j a and j b . Indeed, our candidate reads which is nothing but (60). This definition again gives a negative value for j a > j b , but it also gives β ja,j b + = 0 for "isotropic" configurations, as we would like to demand. Moreover, it is always finite and well-defined regardless of the volume eigenvalue associated to the same pair (j a , j b ).
Even though we consider the definition (A14) to be well-motivated for GFT cosmology, the ambiguity we have highlighted here introduces, at least in principle, an additional uncertainty into the cosmological interpretation of GFT models. To quantify this uncertainty we can compare the various definitions we have discussed. To do this we plot in figure 14 two-dimensional slices of the anisotropy dependence on the spins j a and j b , first for constant j b = 1 2 and along the plane j a = j b . The table below shows the maximum volume eigenvalues (A8), expressed in units of l 3 0 , and the corresponding (dimensionless) anisotropies (60) for given spins j = (j a , j a , j a , j b ). Even though we may use larger spins in some of our calculations, we only give values up to j a = j b = 7/2.