{'fraction_non_alphanumeric': 0.027217206585236325, 'fraction_numerical': 0.007611966719773411, 'mean_word_length': 4.646239380309845, 'pattern_counts': {'":': 1, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Toward matter dynamics in spin foam quantum gravity

Any approach to pure quantum gravity must eventually face the question of coupling quantum matter to the theory. In the past, several ways of coupling matter to spin foam quantum gravity have been proposed, but the dynamics of the coupled matter-gravity system is challenging to explore. To take first steps towards uncovering the influence quantum matter has on spin foam models, we couple free, massive scalar lattice field theory to a restricted, semi-classical 4d spin foam model, called quantum cuboids. This model can be understood as a superposition of hypercuboidal (and thus irregular) lattices. Both theories are coupled by defining scalar lattice field theory on irregular lattices via discrete exterior calculus and then superimposing these theories by summing over spin foam configurations. We compute expectation values of geometric and matter observables using Markov Chain Monte Carlo techniques. From the observables, we identify a regime in parameter space, in which the spin foam possesses a finite total volume and looks on average like a regular lattice with an emergent lattice spacing dependent on the mass of the scalar field. We also measure the 2-point correlation function and correlation length of the scalar field in relation to the geodesic distance encoded in the spin foam. Our results are consistent with the correlation function of ordinary scalar lattice field theory defined on a fixed regular lattice with the emergent lattice spacing and the same mass. We conclude that in this regime of the model, the scalar field is not sensitive to the fluctuations of the spin foam and effectively behaves as if it is defined on a fixed regular lattice.


Introduction
The definition of a consistent framework of gravity and quantum matter is one of the strongest reasons to search for a theory of quantum gravity. Indeed, at the classical level the dynamics of gravity and matter are deeply intertwined and it is expected that the same holds in the quantum regime of both theories. Therefore, all approaches to quantum gravity must eventually face the question how to include quantum matter in their approach, in particular those that describe pure quantum gravity. The purpose of this article is to make progress in this direction in spin foam quantum gravity and present for the first time expectation values of geometric and matter observables in a restricted setting.
Spin foam quantum gravity [1,2] is a background independent approach to quantum gravity, frequently referred to as the covariant / path integral formulation of loop quantum gravity [3]. The starting point is Plebanski's formulation of general relativity [4], as a constrained topological quantum field theory, where we first discretize topological BF theory on a combinatorial 2-complex [5] and then impose the constraints. This theory decorates the discretization with group-theoretic, algebraic data encoding quantum geometric quantities, e.g. areas of triangles and (fuzzy) shapes of tetrahedra. In the path integral one then sums over these data to describe a sum over geometries. One of the most widely used and explored modern spin foam models for 4d gravity is the Engle-Pereira-Rovelli-Livine / Freidel-Krasnov (EPRL-FK) model [6,7,8], which has the particular advantage of connecting on its boundary to the kinematical loop quantum gravity Hilbert space.
Different scenarios on how to include matter have been explored in loop quantum gravity and spin foams, e.g. unification scenarios in which a large symmetry group encodes both gravity and matter [9] or scenarios in which matter degrees of freedom are used to deparameterize the system [10,11,12,13]. In many cases, also in this work, one considers adding matter "on top" of the pure quantum gravity theory, here the spin foam. The idea is to use the geometry encoded in a spin foam as the (fluctuating) space-time on which matter degrees of freedom and their interactions are defined. These works include gauge fields [14,15,16,17] and fermions [18,19]. Moreover, a spin foam model coupled to a massless scalar field exists [20], which was derived from loop quantum gravity coupled to a scalar field following similar works in loop quantum cosmology [21,22]. Defining such coupled systems is vital to investigate how quantum matter and gravity mutually affect each other, and has given rise to interesting questions in other approaches to quantum gravity. One example is the question whether there are (unexpected) constraints from quantum gravity on the possible matter content of the universe or whether quantum gravity effects significantly influence the matter theory, e.g. the mass of elementary particles. Both effects have e.g. been seen in asymptotic safety [23], see e.g. [24,25,26]. In discrete approaches, studying such systems is challenging but possible, e.g. recently in Causal Dynamical Triangulations [27] a scalar field was coupled to the gravitational theory and the authors observed a matter driven change of space-time topology [28,29]. In spin foams the dynamics of such coupled systems, in particular in the deep quantum geometric regime, is barely known, yet more humble goals might be within reach. Indeed, such matter-gravity systems should possess a regime in which one recovers an effective theory that can be understood as a quantum field theory defined on a fixed (emergent) background space-time. Such a result would be a key consistency check for the theory.
However, the dynamics of matter coupled to spin foams is challenging to explore for different reasons. One of them is the complexity of spin foam quantum gravity itself, e.g. the computation of its fundamental amplitudes. Fortunately, in recent years significant progress has been made to overcome this challenge through numerical means, e.g. by developing an algorithm to explicitly compute the vertex amplitude [30,31,32,33], the invention of effective spin foams, which utilize the semi-classical approximation of spin foam amplitudes [34,35,36,37], the usage of Lefshetz thimbles making Monte Carlo algorithms applicable to oscillatory spin foam path integrals [38] and the contributions from complex critical points in spin foams with multiple simplices [39]. Despite these advances, performing pure spin foam calculations beyond small triangulations is challenging, not to mention adding matter dynamics on top. Instead, we opt for a restricted, semi-classical spin foam model, called quantum cuboids defined in [40], which is accessible for numerical simulations and whose amplitudes are derived from the EPRL-FK model defined on a hypercubic 2-complex [41]. Essentially, one can understand this model as a superposition of flat, hypercuboidal (and thus irregular) lattices weighted by spin foam amplitudes. In the last few years, this model gave the first chance to study important questions in spin foam quantum gravity, e.g. renormalization [42,43,44] and the spectral dimension [45]: In [46,47] the renormalization group flow of cuboid spin foams was computed, which showed indications of a UV attractive fixed point and a second order phase transition. These results were later generalized to frusta [48,49], generalizations of cuboids to cosmological setting allowing for non-vanishing curvature, which confirmed these findings in a more general context [50]. Moreover, in [51] the spectral dimension, an effective dimension measure derived from a diffusion process, for quantum cuboids was computed, which can be lower than four due to the superposition and scaling of cuboid geometries of different size. Thus, while quantum cuboids are no realistic model of quantum space-time, they pertain features of spin foam models and provide an ideal test scenario to define and study suitable observables. In this work, we will e.g. define an observable, the 2-point correlation of the scalar field in a relational manner [52,53,54], where we measure both the scalar field correlations as well as their distance encoded in the spin foam state. As a matter system we consider a free, massive scalar field defined as a lattice field theory on the cuboid spin foam. To do so, we interpret a spin foam configuration as an irregular lattice on which we define scalar lattice field theory via discrete exterior calculus [55,56]. Then, by superimposing spin foam states we also superimpose lattice field theories defined on these configurations. At the level of the amplitudes, the spin foam is unaffected by the addition of the scalar field, while the latter is sensitive to the spin foam only via the geometry of the spin foam state it is defined on. We consider this a "minimal coupling", but it is certainly not free of ambiguities. Due to the non-oscillatory nature of quantum cuboid amplitudes, which even persists into the deep quantum regime [57], we can study the whole system with a Markov Chain Monte Carlo algorithm and compute expectation values of observables. Thus, this model gives us the first chance to study the dynamics of a matter theory coupled to a 4d spin foam model and learn valuable first lessons about such systems.
This article is organized as follows: we open in section 2 with a brief review of free scalar field theory defined on a regular lattice. In section 3 we generalize this field theory to irregular lattices using discrete exterior calculus. Section 4 briefly introduces spin foam models and the cuboid restriction. We define the coupled system, scalar lattice field theory on quantum cuboids, in section 5 and describe the Monte Carlo algorithm to compute its expectation values in section 5.3. The results are detailed in section 6 and we close with a summary and discussion in section 7.

Scalar field theory on a regular lattice
As an introduction to lattice field theory, we start by briefly discussing the free, massive scalar field defined on a regular, hypercubic lattice. We start from the continuum (Euclidean, Wick-rotated) action: and define the partition function as: To make this expression well-defined, we regularize it by introducing a discretisation, a regular lattice with lattice constant a. On the lattice, the continuum scalar field φ(x) is replaced by the field φ x located at the coordinates x of the vertices of the regular lattice. We discretize the derivatives by approximating them by difference quotients of fields sharing the same edge, but this definition is not unique and different choices differ by lattice artifacts of O(a). Instead, we consider the lattice Laplacian (by performing an integration by parts and dropping boundary terms due to periodic boundary conditions): where e µ denotes a unit vector in µ direction, and ±µ denotes that we are summing over positive and negative directions for each direction µ.
The integral over all points in the continuum is replaced by a sum over all vertices, d 4 x φ(x) → a 4 x φ x , where a 4 is the 4-volume of a lattice hypercube. Here x denotes the coordinates of the vertices of the lattice and φ x := φ(x). The functional integral in this discrete setting is defined as a product of regular integrals, D[φ] → x dφ x . Combining all these ingredients, we find the discrete lattice action.
where we have used the periodicity of the lattice to rewrite the kinetic part of the action. This action can be brought into a simpler form, by enumerating the vertices of the lattice and writing the field configuration as a single N -dimensional vector φ, where N is the number of lattice vertices. Then, we express the scalar field action as the contraction of a matrix K by two vectors of field configurations: Given this representation, computing observables such as the 2-point correlation function φ n φ m is straightforward and given by the inverse of K, K −1 nm . The correlations of scalar fields at different vertices drops off exponentially with the distance d(n, m) between them i.e., φ n φ m ∼ e − d(n,m) ξ , where ξ denotes the correlation length. For the massive scalar field, this correlation length is inversely proportional to the mass M of the scalar field, ξ = M −1 , i.e. the larger the mass the faster the correlations die off as we increase the distance. Conversely, we can estimate the mass of the scalar field from the exponential decay of its correlations.

Taking the continuum limit
In lattice field theory, the lattice plays the role of a regulator and particular care is necessary to make sure that the results are independent of this regularization. To this end, one investigates whether a continuum limit can be taken, i.e. the regulator can be removed. Naively, this limit corresponds to taking a → 0, but it can only be taken on a second order phase transition in the phase diagram of the theory. Only on such a transition scale invariance is guaranteed and the limit a → 0 is well defined. Note however, that such a phase transition is not possible for finite volume. The specifics of this procedure go beyond the scope of this article, and we refer to the literature of lattice field theory.
For the free, massive scalar field, the continuum limit can be readily taken. Note that the lattice spacing a also determines the dimension of analytic operations, i.e. dimensionful quantities can be rescaled by (powers of) a to turn both fields and coupling constants into dimensionless quantities. In 4d, we define dimensionless fields ϕ = φ a , and the scalar field action is written in terms of a new matrixK: where we have defined the lattice massM (a) = aM . Thus, the scalar field action contains no dependence on the lattice scale besides the lattice massM . Taking the continuum limit corresponds to taking the lattice massM → 0, while keeping M =M a constant. This leads indeed to a second order phase transition, since the correlation lengthM −1 , i.e. the physical correlation length expressed in lattice units, diverges. While this is a trivial example of taking the continuum limit, it demonstrates that free, massive scalar field theory has only one free parameter, the mass M of the scalar field. Hence, it should also be a suitable test case for coupling matter to a discrete approach to quantum gravity.

Scalar matter on an irregular discretization
In order to couple scalar lattice field theory to spin foams, we generalize its definition to irregular lattices, which includes both their geometry and combinatorics. Therefore, the construction we present is applicable e.g. to triangulations. To do so, we use discrete exterior calculus [55,56,58]: the idea is to discretize a continuum action expressed in terms of differential forms by smearing p-forms (their duals (d − p)-forms) over p-simplices (dual (d − p)-simplices) respectively. On a manifold M a real scalar field φ with mass M is a 0-form (φ ∈ Ω 0 (M)) with the action where dφ is the exterior derivative of φ and * denotes the Hodge dual.
Let us now consider a simplicial discretization K of M. In this discrete setting the role of a p-manifold is played by p-chains, i.e, formal linear combinations of p-simplices (σ p ) in K generating the vector space of p-chains c ∈ C p (K). We write an element of C p (K) as Similarly linear forms on K are represented by p-cochainsc ∈ C p (K) written as: This identification between the continuum linear forms and the discrete p-cochains assumes that the discretization process involves an integration of p-forms over p-simplices and in the discrete setting integration is replaced by the evaluation of a cochain on a chain ( c|c ). Thus, a p-form smeared on a single p-simplex may be written as Here ω|σ p is the integrated value of the p-form ω on σ p and contains information about the volume of σ p , whereas the coefficients ω σp are the discrete components of the p-form. Therefore the scalar field can be discretized by smearing on 0-simplices (σ 0 ) as: where the vertex volumes are assumed to be trivial in the last equality and the coefficients φ σ0 are the discrete components of the field. In order to write down the action for this discrete scalar field we need to define the concept of a discrete Hodge dual. Let us then consider the complex * K that is dual to K. The complex * K consists of (d − p) cells σ p dual to the primal cells σ p with orientation and cellular structure induced from the orientation and adjacency matrix of K. The duality between the primal and dual simplicial complexes induces a duality between the chains of the primal (c ∈ C p (K)) and dual complexes ( c ∈ C p ( K)) (as well as the co-chains of primal and dual complexes). That is, Along with the identification between p-forms and cochains on K, this allows the definition of the discrete Hodge dual of a p-form as At the level of coefficients the duality reads: For the smeared p-forms this implies: With these definitions, we can define a pairing of the p-chains on K and its dual on K as σ p | σ p = δ p,p and a completeness relation p∈K |σ p σ p | = id [58]. In order to define an action for the discrete field we need the definition of the discrete exterior derivative. The discrete exterior derivative is defined in order to satisfy a discrete version of the Stokes' theorem dω|c = ω|δc , such that the exterior derivative d is the natural dual of the boundary operator δ as in the continuum. The boundary operator is defined by its action on a p-simplex with vertices [v 1 ...v p ]: where the first equality is the alternating sum over all (p − 1)-simplices obtained by dropping a vertex in σ p (the parenthesis denotes the vertex that has been dropped) and the second equality is just a compact rewrite.
Thus we can write the discrete analog of the exterior derivative of a p-form ω as: We now have all the ingredients to write the discrete action for the scalar field. It takes the form: where the second equality arises after integration by parts and setting the boundary term to zero, e.g. when using periodic boundary conditions. Interested readers will find further details in [55,58]. Using equations (12), (16) and (19) repeatedly 1 the discrete action can be written as In 4d, the case we are interested in, V σ1 is the length of the edge σ 1 , whereas V σ1 is the volume associated to its dual 3-cell. V σ0 denotes the 4-volume of the 4-cell dual to a vertex. After a little rearranging we can write this is in the form S D = 1 Note that the different components of K ij associated to the Laplacian and the potential have the same scaling behaviour in terms of lengths as the scalar field theory action discretized on a regular lattice in section 2 2 . Indeed, (22) simplifies to (6) if defined on a regular lattice, where the dual lattice is simply given by the same (shifted) lattice. The definition of the scalar field action (22) can be used for general cellular complexes, but its definition is not unique, e.g. it is a priori not clear whether to smear the scalar field over vertices of the complex or over vertices of the dual complex. Additionally, there is freedom in constructing the dual lattice and thus defining the dual volumes.
In the next section we briefly introduce spin foam models and the restricted model we will couple with scalar lattice field theory.

Restricted spin foam models
Coupling a lattice field theory to spin foam models and examining this system numerically in full generality is a daunting task. Thus, in order to get a first insight into such spin foam-matter systems and to determine which methods are suitable to investigate such models, we will study a spin foam path integral restriced to so-called "quantum cuboids" [40]. In this section, we give a brief introduction to spin foam quantum gravity [2] and define this restricted model. Similar to lattice field theories, spin foam models are defined on a discretisation, more precisely a 2-complex, which is a collection of vertices v, edges e and faces f . Frequently this 2-complex is chosen dual to a triangulation, but it can be more general [41]; here we choose it to be dual to a (hyper)cubulation. The discrete geometry of the spin foam is then encoded in group theoretic data assigned to the 2-complex. To each face we assign an irreducible representation of the underlying symmetry group, while we assign an invariant tensor, called an intertwiner, to each edge.
Let us explain these concepts for the concrete case of SU(2) BF theory [5]: an edge e in the 2-complex is shared by several faces f ⊃ e. Each face f carries an irreducible representation j f , where we assign to the edge e a vector space V j f ⊗ V * j f for each f ⊃ e. Thus, each edge e is endowed with a tensor product of representation spaces from the faces by which it is shared, where we assign the various representation spaces and their duals to the source and target vertex of the edge 3 . However, not the entire vector space is permitted: to each edge we assign a projector onto the subspace invariant under the action of SU(2), called the Haar projector: P j1,...,jn m1,...,mn,n1,...,nn := where D j denote Wigner matrices in representation j and magnetic indices m, n with −j ≤ m, n ≤ j. These projectors can be expressed as a (finite) sum over orthonormal intertwiners ι d P j1,...,jn m1,...,mn,n1,...,nn = d |ι d ι d | .
As an example, if the edge is 4-valent as in a 4d triangulation, the intertwiner basis elements are labelled by an additional SU(2) spin arising from expanding the 4-valent intertwiner into two 3-valent ones using recoupling theory. Given this projector property and the fact that an edge connects two vertices, we associate one intertwiner to each vertex of the edge. Then, at the vertex the intertwiners get contracted according to the combinatorics of the vertex. This spin network evaluation defines the vertex amplitude. Additionally, spins and intertwiners permit a quantum geometric interpretation: we think of an edge to be dual to a polyhedron, with as many faces as the valency of the edge. The areas of the faces are given by the spins associated to them, while the intertwiner encodes the shape. For a tetrahedron, an intertwiner in the orthonormal basis is given by five parameters, four areas and one intertwiner label. The last label corresponds to the area of a parallelogram inside the tetrahedron according to the recoupling scheme. This does not uniquely determine the shape of the tetrahedron, which would require six parameters, e.g. six edge lengths. Thus, this is often referred to as a quantum tetrahedron [59], and similar arguments also apply to more general polyhedra.
Modern 4d spin foam models like the Euclidean EPRL-FK model [7,8] we are considering in this article use BF theory as their starting point. In order to break the topological nature of BF theory and arrive at a gravitational theory, simplicity constraints are imposed. To each polyhedron we assign a time-like normal 4 , such that -as seen from 4d Minkowski space -the polyhedron is entirely space-like. This normal singles out a subgroup of the full symmetry group (Spin(4) in Euclidean and SL(2, C) in the Lorentzian signature), which stabilizes the normal. This subgroup is isomorphic to SU(2), and we construct states of the 4d spin foam model by embedding SU(2) states into states of the 4d symmetry group. This leads to restrictions on the representation labels depending on the Barbero-Immirzi parameter γ compared to 4d BF theory. Furthermore, on the boundary we can connect to the kinematical Hilbert space of loop quantum gravity via projected spin network states [61,62].
For the Euclidean EPRL model for γ < 1, the symmetry group Spin(4) ∼ = SU(2) × SU(2) and its representations (j + , j − ) are labelled by two SU(2) representations j ± ∈ N 2 . They are related to the same SU(2) representation as follows: where j ∈ N 2 . This immediately requires γ ∈ Q for this map to be non-trivial, and is considered a pathology of the Riemannian model, which is absent in the Lorentzian theory.
To implement this construction for the intertwiners, we define a boost map Y γ consisting of two parts. For an SU(2) intertwiner ι, we isometrically embed each vector space V j into the unique component appearing in the Clebsch-Gordan decomposition of V j + ,j − V j + ⊗ V j − . We call this map β γ j as it explicitly depends on γ. Note that the resulting tensor is not necessarily a SU(2) × SU(2) intertwiner; we thus act with the Haar projector P. For an N -valent intertwiner Y γ reads: Eventually, the vertex amplitude A v is defined as the contraction of intertwiners associated with the vertex v, also called the vertex trace.
To complete the definition of the model, we must assign amplitudes to the edges e and faces f . The edge amplitude A e is given by the norm of the associated intertwiner ι e : To the faces we assign the dimension of the representations j + , j − up to an exponent α: We introduce this parameter α by hand to reflect an ambiguity in defining the face amplitude of the model. Common choices are either the SU(2) representation 2j f + 1 or the SU(2)×SU (2) representation The former choice has the advantage of being invariant under face and edge subdivisions [63]. Therefore, we understand the parameter α as a modification of the path integral measure. In the quantum cuboid model we introduce below, these face amplitudes can be rewritten as the volume of 4d hypercuboids to a power of a multiple of α. This is similar to path integral measure choices proposed in quantum Regge calculus [64]. The choice of α has a profound impact on the scaling behavior of the amplitudes and played a crucial role in the coarse graining [46,47] and the spectral dimension [51] of restricted 4D spin foams models. We will also see that it has a significant impact on the coupled matter spin foam system.
Finally, the spin foam partition function is defined as: where we sum over all spin foam states given by assignments of representations and intertwiners.

Cuboid spin foams
Spin foam amplitudes are well-defined, yet the complexity of the theory makes actual calculations demanding. This is due to the sheer number of configurations one has to consider even for small 2complexes as well as the challenge to compute its fundamental amplitude beyond asymptotic methods.
To at least explore a subset of the full spin foam path integral, restricted semi-classical models were introduced in [40,48,65], and simplify the theory in different ways: • Hypercubic combinatorics: We choose the 2-complex to be dual to a 4d cubulation. This has the advantage that the combinatorics is regular, e.g. each face is bounded by four edges and each edge is 6-valent etc. The initial motivation for this choice in [40] was to simplify renormalization / coarse graining [46,47,50], but here it will be helpful to relate our model to lattice field theory. Figure 1: The shape of a cuboid is encoded in the areas of its faces j i and the the outward pointing normal vectors.
• Intertwiner restriction: Spin foam models with hypercubic combinatorics are generically more complicated than their counterparts defined on (dual) triangulations. This is due to the fact that more data are required to specify the building blocks, e.g. a hypercubic vertex has eight edges and 24 faces compared to the five edges and ten faces of a vertex dual to a 4-simplex. This is counteracted by only allowing specific coherent intertwiners that are peaked on a cuboid [40] or a frustum shape [48]. These choices also restrict the representations, reducing the complexity of the model further.
• Semi-classical amplitude: The first numerical computation of the vertex amplitude for hypercubic combinatorics has only recently been achieved beyond smallest spins [57], however at great numerical costs. Thus it is currently out of reach to explore larger 2-complexes using the full amplitude; instead we will use its semi-classical approximation derived in [40,48].
• Continuous variables: The semi-classical approximation is typically derived in the large j-limit for coherent boundary data, where the representation labels of the vertex amplitude are uniformly scaled up. In this limit, the discreteness of representations is barely noticeable and we approximate them as continuous variables.
Most of these assumptions are drastic simplifications, restricting the full spin foam model to a subset of the gravitational path integral that is more readily accessible by numerical means. Note that this does not fix a background geometry a priori, as we sum over all configurations permitted by these restrictions. However, as we will see in this article, the system might be dominated by particular configuration (plus fluctuations around it), which one can interpret as the emergence of a dynamical background. Clearly, the restricted models do not capture all properties of spin foams, and these assumptions must be eventually removed in future research. The restriction of intertwiners is the final piece left to define the concrete restricted model. We choose the intertwiners to be coherent Livine-Speziale intertwiners [66] peaked on the shape of cuboids, dubbed cuboid intertwiners. Coherent SU(2) intertwiners (and also the boosted intertwiners) are parametrised by spins j i and normal vectors n i ∈ S 2 assigned to their faces, which agree with the areas and outward pointing normals of the classical polyhedron on which they are peaked. For cuboids, opposite faces carry the same area and opposite normals are anti-parallel, while each normal is perpendicular to the normals of its four adjacent faces, see also fig. 1. We define the intertwiner as follows: The states |j i , e i ∈ V j denote SU(2) Perelomov coherent states [67], i.e. maximum (or minimum) weight states on which we act with SU(2) representation matrices. The state |j, n is defined as g n |j, j =: g n |j, e z , where g n is the group element encoding the rotation e z → n. By e i we denote the unit vectors for Cartesian coordinates. The group integration, also called group averaging, ensures the invariance of the intertwiner under SU(2) transformations. Boosting these intertwiners for γ < 1 is straightforward: each coherent state |j, n gets mapped to a tensor product of coherent states |j + , n ⊗ |j − , n for the representations j ± . The SU(2) × SU(2) group averaging is then implemented as two group integrations, one for all "+" and one for all "−" states. Thus, the resulting intertwiner is a tensor product of two SU(2) intertwiners, one for "+" and one for "−" labels, peaked on the same coherent data { n}. See [68,69] for a more detailed derivation and explanation.

Semi-classical approximation of amplitudes
The vertex amplitude for cuboid spin foams is given by the evaluation of the spin network in fig. 2. For coherent boundary data, it reads: We denote the intertwiners in the amplitude by a, b, and have one group integration for each. Then, to each link connecting two intertwiners, we assign an inner product of coherent states labeled by their respective normal vectors: n ab denotes the normal vector to the face ab belonging to the polyhedron a and n ba vice versa. Moreover, we have used the property of coherent states (as maximum weight states) j, n 1 | · |j, n 2 = ( n 1 | · | n 2 ) 2j , where the latter expression is given in the fundamental spin 1 2 representation. For γ < 1 the amplitude factorizes over the "+" and "−" labels, such that we can write it as the product of two SU(2) BF theory vertex amplitudes for specific representations.
These group integrations are highly oscillatory integrals for large j ± ab . Thus, the vertex amplitude is approximated by performing a stationary phase approximation of the group integrals: to do so, the inner products of SU(2) coherent states are exponentiated to define an action: The vertex amplitude is dominated by the critical points of S ± if all spins j ± ab are large. These critical and stationary points are derived by varying the action with respect to the dynamical variables, here group elements, and enforcing that Re S ± = 0. Variation with respect to the group elements enforces closure of the coherent data b =a j ± ab n ab = 0 ∀a, while the reality conditions determine how the polyhedra are glued together. Again see [40] for more details.
With the critical and stationary points found, we approximate the vertex amplitude A ± v as follows: where we sum over the critical point solutions g a . The action S ± and the determinant of the Hessian matrix H are evaluated on these solutions. We get one factor of √ 2π for each integration; we integrate over seven SU(2) group elements (one integration can be absorbed due to gauge invariance), which are each three-dimensional. For cuboid boundary data, this formula simplifies drastically [40]: The action S ± has two critical points with a multiplicity of 2 7 due to a discrete symmetry. We also obtain one factor of (16π 2 ) −1 for each SU (2) integration. Note that the action evaluated on the critical points, which corresponds to the Regge action [70] associated with the 4d polytope, vanishes. This has two important consequences: firstly, for larger complexes all deficit angles identically vanish such that we only find flat geometries. Essentially we are gluing together flat geometries in a flat way, corresponding to a cubulation of flat space-time. Secondly, due to the vanishing of the Regge action the vertex amplitude is not oscillatory and determined by the scaling behavior encoded in the determinant of the Hessian. Thus, we can also pull the dependence on the Immirzi parameter γ in front of the amplitude. The non-oscillating behavior of the cuboid amplitude even persists in the quantum regime [57].
Due to the symmetry of the cuboid intertwiners, i.e. the fact that opposite sides of the 3d cuboids carry the same spin, the amplitude B v is a function of six SU(2) spins. This is peculiar since a 4d hypercuboid is determined by just four edge lengths; this mismatch comes about as follows. Each 3d cuboid is specified by three spins, which determine the six areas of its faces. As long as these areas are non-vanishing, we can uniquely translate three areas into three edge lengths for each cuboid. Depending on the choice of six spins, these edge lengths may however not agree, such that the faces shared between two cuboids have the same area but their shapes do not match. In larger 2-complexes, this can lead to torsion effects where a loop may not close due to the mismatch of lengths in neighbouring hypercuboids [40].
In the following we consider the subset of spin configurations for which the edge lengths, as seen by different cuboids, agree and shape matching is ensured. The non-shape matching is related to the so-called volume simplicity constraint, which is one of several simplicity constraints imposed at the discrete level. In a classical 4-simplex, expressed as a bivector geometry, the volume of this 4-simplex is spanned by any pair of bivectors associated to triangles that only share one vertex. Again, for a 4-simplex this constraint is automatically satisfied once all other simplicity constraints are implemented, which is the reason why it is not enforced in the original EPRL model. However, this no longer holds for more general 2-complexes [71,72,65], generically giving rise to non-matching shapes [73,74].
In [40] also the semi-classical expressions for the edge and face amplitudes were derived. The edge amplitude for the coherent cuboid intertwiners is a simple functions of its three spins: Similarly, the face amplitude reads: As a final point, due to the regular combinatorics of the hypercubic spin foams, we can combine vertex, edge and face amplitudes into a single amplitude associated to the vertex. Each edge is shared by two vertices, while each face contains four vertices; hence we define:

From spins to lengths
Let us briefly discuss the transition from spins to lengths. Firstly, we must restrict the spin configurations of the hypercuboid such that the areas of its rectangles arise from four edge lengths. To do so, we express two of the spins per hypercuboid as functions of the remaining four spins. However, since we are originally integrating over all spins in the path integral, we effectively gauge fix the redundant spins and account for this by including a Fadeev-Popov determinant. In the second step, we obtain a Jacobian for the change of variables from the four remaining spins to four edge lengths. Both steps are explained in detail in [47]. Eventually, due to the regular combinatorics, we combine face, edge and vertex amplitudes as well as Jacobian and Fadeev-Popov determinant into a common amplitude that we associate to the verticesÂ v , which is a function of four lengths of the hypercuboid and the parameter α 5 . This vertex amplitude is a homogeneous function of degree 24α − 14: The spin foam partition function is then given by: In this partition function we integrate over all possible lengths assignments to the 2-complex, which are weighted by the spin foam amplitudes. Note that the amplitudeÂ v is positive and non-oscillatory. Thus, after introducing upper and lower cut-offs for the lengths l i , we can use the partition function to define a probability distribution suitable for Markov Chain Monte Carlo techniques. This probability distribution then determines the likely length configurations of the spin foam.

Scalar matter coupled to spin foams
In the previous sections 3 and 4 we have defined the two models individually that we couple in this work, i.e. a free, massive scalar field defined on a general, irregular lattice and cuboid spin foams, which describe superpositions of discrete flat space-times. Our idea for coupling both theories is as follows: the discrete scalar field action is a function of various (dual) volumes of the underlying irregular lattice. Given a specific spin foam state, here a configuration of lengths {l i }, we define the volumes in the scalar field action as functions of these lengths, such that the spin foam describes the space-time the scalar field is defined on 6 . A similar idea is also used in [14] to couple Yang-Mills theory to the Barrett Crane spin foam model [75]. Then, we integrate / sum over all field configurations and lengths configurations to define the partition function.
In this setup, we solely modify the matter part of the system without altering the gravitational side in analogy to general relativity. There, the Einstein Hilbert action remains unchanged, while the matter action has an explicit dependence on the metric in contrast to its definition in flat space-time. Thus we will refer to our ansatz as "minimal coupling": the spin foam amplitudes remain unchanged, while we include an explicit dependence on the geometry, here given as a function of the edge lengths, in the matter action. Explicit modifications of the gravitational part are conceivable, but they are not subject of this work. However, the coupling of matter and spin foams is not uniquely defined. One important choice concerns "where" the matter degrees of freedom are placed: the actual discretisation (here the cubulation) or the dual 2-complex. We can either decide to place the scalar field on the vertices of the cubulatation or on the vertices of the 2-complex, dual to each 4d hypercuboid. In the example we are considering here, where the combinatorics of the dual 2-complex is the same as the lattice and we are imposing periodic boundary conditions, we do not expect this choice to have a major impact on the dynamics. Nevertheless, these choices lead to different theories because in general the matter theories are not dual to themselves. Moreover, if the space-time possesses a boundary, this choice affects the boundary Hilbert spaces we associate to them.
In this work, we choose to place the scalar field on the vertices σ 0 of the cubulation, and infer the volumes from the lengths describing the spin foam. Volumes of the lattice itself, e.g. the length V σ1 , are computed directly from the lengths. Analogously, we directly compute the dual volumes from the dual lengths; the lattice dual to a hypercuboidal lattice is again a hypercuboidal lattice. Since the spin foam model we are considering is restricted to hypercuboidal geometries, we simply place the vertices of the dual lattice in the center of the hypercuboids. Thus, the dual lattice is again a hypercubulation with the dual edges connecting the centers of neighbouring hypercuboids, from which we straightforwardly compute the dual volumes. All of these definitions are simple and straightforward due to the rigid right angles of the geometries given by cuboids. Note that in general the construction of a dual lattice is not unique, which thus influences the definition of the matter action, see e.g. [76] for further details. Finally, we define the partition function of the coupled system as: Here we denote the edges of the cubulation by i and the vertices of the cubulation by a, b. σ 4 are the hypercuboids and σ 0 are the vertices of the cubulation. M denotes the mass of scalar field.

Qualitative behavior from regular lattices
It is difficult to extract qualitative behaviors from the composite system by analytical methods alone. To gain some initial insight into the qualitative features of the model, we set all lengths of the spin foam equal (yet dynamical), such that we consider a superposition of regular lattices of different sizes. As discussed above, the spin foam vertex amplitude has a scaling behavior of l 24α−14 . Additionally, for fixed l we perform the Gaussian integrations of the scalar field, which leads to a factor in the integrand of the following form (in four dimensions and N lattice sites): a i denote numerical factors which are not relevant for our argument here. Thus, the probability distribution for the length l is given by: Thus, we observe two counteracting behaviors of the scalar field and the spin foam (for sufficiently large α): the scalar field leads to a polynomial suppression of large lengths, which gets stronger as we increase the mass M . Conversely, the spin foam amplitude eventually favors larger lengths as we increase α. Thus we can readily identify two qualitative regimes of the model: • Finite M and small α: α is too small to counteract the suppression of large lengths by the scalar field. Thus, the amplitude diverges as we approach l → 0 and the smallest lengths dominate. To explore this regime in the actual simulations we introduce a lower cut-off on the lengths. Note that such small lengths are outside the validity of the semi-classical approximation of spin foam models.
• Finite M and large α: α is so large that it dominates over the suppression of large lengths by the scalar field. Then, the amplitude diverges for l → ∞ and largest lengths dominate. Thus we introduce also an upper cut-off in the simulations.
Note that for fixed lattice size the values of α when either the smallest or the largest lengths dominate are independent of M . Moreover, we expect these two α-dependent regimes to also exist in full spin foam models, albeit for different parameter values. Spin foam vertex amplitudes are (in the semi-classical limit) polynomially suppressed under uniform scaling of their representations [69,68], while the scaling behavior of the scalar field in terms of lengths is encoded in its action. Thus, while other effects will appear under lifting the cuboid restrictions, e.g. oscillatory spin foam amplitudes, the respective scaling behaviors persist, which depends on the combinatorics of the 2-complex. The form of the denominator additionally suggests that an intermediate region between small and large α exists, where the length l is finite. This is due to the various summands in the denominator: qualitatively, for small lengths the term with smallest exponent dominates while clearly the opposite is true for large lengths. If we choose a value for α such that the scaling behavior of the spin foam amplitude is between these two extremes, the probability distribution should have a peak at finite lengths, whose position and width is a function of α and the mass M . This is indeed the case as we show in fig. 3.
From the simple form of the probability distribution we can straightforwardly derive the interval in α for which we expect finite lengths. For small length l the probability distribution scales like l N (24α−15)−1 , whereas for large length it scales as l N (24α−16)−1 . Hence, for N 1 the interval of interest is 5 8 = 0.625 ≤ α ≤ 2 3 = 0.66. Of course, all of these considerations are strictly only true if all lengths are the same, but they provide us valuable guidance and insights into how the dynamics change once all lengths in the spin foam are dynamical, i.e. we consider a superposition of irregular lattices.

Observables of the coupled system
The expectation values of observables are defined analogously to other path integrals or statistical systems. Given an observable O, which depends on the configuration variables {l i } and/or {φ a }, its expectation values are straightforwardly defined as: Here we have suppressed the dependence of the spin foam amplitudes and the matter action on the configuration variables to keep the expression readable. In this work we are considering observables of either the spin foam model or the scalar field. For the spin foam, we examine simple geometric observables like the total 4-volume V or individual lengths l i and their variances respectively, while for the scalar field we are interested in the correlations φ a φ b of the scalar field located at vertices a and b. These observables are interesting for the following reasons. The geometrical observables, like the total volume of space-time, are relevant as they are a direct indicator for the scaling behavior of the coupled model. Cuboid spin foams by themselves usually have diverging edge lengths / 4-volume for values of α above ∼ 0.58; then the amplitudeÂ v is a homogeneous function of positive degree and largest lengths dominate. However, we expect that the scalar field instead favours smaller edge lengths: the kinetic part of the action scales with ∼ l 2 , while the potential part with ∼ M 2 l 4 . Hence the action in the fields is smaller if the lengths are smaller, which gets more pronounced as the mass is increased. This interplay between the spin foam and matter part is intriguing.
On the other hand, a free, massive scalar field is solely determined by its 2-point correlation function; any correlations of higher degree can be understood as arising from 2-point correlations. Moreover, these correlations drop off exponentially with the distance between vertices and depend on the mass of the scalar field. Thus, it is the ideal observable to study here, yet, due to our setup, the interpretation is more subtle. Consider the correlator φ a φ b for fields located at vertices a and b in lattice field theory on a regular lattice. Here, the lattice inherits the notion of distances as a discretisation of the flat background space-time, such that the vertices a and b are space-time events. We readily know the geodesic distance between these vertices and use this information to compute the correlation length from the correlations between all the vertices of the lattice. However, in a spin foam and more generally in a background independent setup, the vertices a and b have no inherent meaning; they are merely unphysical coordinates/labels. Their distance is not fixed by the discretisation, in fact we are integrating over superpositions of distances between these two vertices. Nevertheless, in addition to the correlations of the two fields located at a and b, we also need the distance the vertices a and b are apart to obtain a physically meaningful relational observable [52,53,54].
In our Monte Carlo simulations, we use the following approach: for each sample of lengths and fields, we compute the correlations of fields for all pairs of vertices and the distances between these vertices for a specific sample. Then we use the collection of data of correlations and distances to determine the correlation length as a non-local, coarse grained observable to learn how the matter theory effectively behaves on the spin foam background. Note that in this approach we cannot control the distance at which we probe the correlations; we can only access the correlations for probable spin foam and scalar field configurations. A suggestion to measure the correlations at arbitrary geodesic distance is described in the introduction of [27] 7 .

Setup for numerical simulations
The partition function for the cuboid spin foam model with minimally coupled scalar degrees of freedom on the vertices of the cubulation is given by (41). In order to study this system we compute observables of the coupled system using Markov Chain Monte Carlo integration. The general idea of Monte Carlo integration is to approximate an expectation by the sample mean of the function of simulated random variables [77]. We are interested in computing the expectation values of observables of lengths and fields as defined in the previous section. We approximate this expectation value with the sample mean of O(l i , φ a ), which is a function of the random variables (lengths and fields) drawn from the distribution We use the Metropolis-Hastings algorithm to sample the lengths and fields. This algorithm is based on a theorem that states that the Markov chain with transition probabilities arising from the Metropolis-Hastings algorithm converges to the target distribution if the chain is ergodic and satisfies detailed balance [78]. In our setting each configuration of lengths and fields is an element of the Markov chain. As this convergence can take time we need to discard some initial samples (called the burn in) in order to "forget" the starting configuration. Although the marginal distribution of the Markov chain converges to the target, it does not mean that the chain converges to a chain of identically independently distributed draws from the target distribution. This is because of the auto-correlation introduced by the proposal distribution leading to a bias in the observed sample means. Thus it is important to sample the chain at intervals longer than the auto-correlation length. This auto-correlation also introduces a trade-off between auto-correlation and acceptance rates. If the acceptance rate is too high the auto-correlation length increases as each consecutive sample represents only a small change in the target distribution. On the other hand a high rejection rate results in slower convergence to the target distribution. For this coupled system it was particularly difficult to balance the acceptance-rejection ratio. Changing both the lengths and fields together leads to high rejection rates. This forced us to settle for smaller changes in the configuration to maintain a reasonable acceptance rate. We chose to either change all the lengths associated to a hypercuboid or the field at one vertex. New lengths were proposed by rescaling where the rescaling factor r i is different for each length of a hypercuboid. Choosing this rescaling factor between 0.5 and 2. gives healthy acceptance rates (40 − 60%) and also allows us to explore several orders of magnitude of edge lengths. Additionally we also restricted the configuration space by imposing an upper (10 4 ) and a lower cutoff (5 × 10 −5 ) on the lengths in order to avoid numerical overflow as the spinfoam amplitude diverges as the lengths go to 0 or infinity 8 . Detailed experimentation showed this did not have any qualitative impact on the results presented in this article. New field values were proposed using where δ was chosen to maintain an acceptance rate between 40−60%. The consequence of these incremental changes is high auto-correlation resulting in a longer burn-in period and the need for more intermediate Monte Carlo steps between samples. In order to determine if the chain had thermalized we plotted the amplitude and the volume for each Monte Carlo step in fig. 4 for different distinctive regions of the model. During the burn-in period we 7 The idea is to define the correlator as a function of geodesic distance. Essentially one integrates the scalar field correlation function over all pairs of space-time points over all metrics and enforces a fixed geodesic distance by a delta distribution. 8 The physical meaning of these cut-offs is different. For the lower cut-off one might choose the minimal non-vanishing area eigenvalue of a face, however also vanishing areas are permitted in the state sum. Typically, the upper cut-off needs to be removed as it breaks gauge invariance, unless one is considering a model incorporating a cosmological constant, see e.g. [79,80,81,82,83]  expect a trend (increasing or decreasing) in the amplitude plots but after thermalization we expect the amplitude to be uniformly distributed over some range.
After thermalization, expectation values of observables can be computed by drawing samples from the Markov Chain at regular intervals (determined by the auto-correlation length). In fig. 5 we plot the auto-correlation over the samples for different parameters. There are four main sources of error in this estimate: • Estimating the state of thermalization of the chain incorrectly.
• Highly correlated samples • The Markov chain takes too long to explore some parts of the distribution and thus these states are not well represented in the sample (the "missing mass" problem).
• Not sampling enough The first two errors are easily controlled for. The first by identifying suitable observables to test for thermalization. We use the full amplitude and the total volume of the spin foam for this. The second error can be mitigated by sampling the Markov Chain at intervals greater than at least twice the autocorrelation length for the observable of interest, where the autocorrelation length is computed after thermalization. This is known as sub-sampling. The "missing mass" problem is harder to identify and control for. A dirty fix is to run multiple simulations each time starting with a different configuration. We tested our algorithm by running multiple simulations for randomly chosen initial configurations for lattices with 81 vertices. But the time cost was prohibitive for larger lattices.
The last problem stems from the fact that the error in the estimate of the expectation value of the observable is dependent on the number of samples. Since sub-sampling minimizes the auto-correlation between consecutive samples drawn from the Markov Chain, we can treat these samples as independent. Then the true expectation value of the observable of interest in terms of the sample mean of a sample of size N can be written as: where O n is the observable value calculated using the n th sample. The second term is the error term in our estimate of O , where σ is the standard deviation of O. Thus the error term can be reduced simply by increasing the number of samples. We must thus balance the time cost of generating the samples with the error in estimating the expectation value by the sample mean. The time complexity of the algorithm scales linearly with the number of Monte Carlo steps and as a polynomial of the number of vertices. Each time we sample the lengths we need to recalculate the amplitude of the coupled system. This is the most expensive operation as calculating the spin foam amplitude is O(N ) for a cubulation of N vertices and the calculation of the scalar field action involves a matrix multiplication which in the worst case results in O(N 3 ) and even in the best case scales much faster than N 2 . We optimize the algorithm to some extent by directly computing the change in the scalar When the observable is only a function of the lengths it is possible to sum over the scalar fields analytically and then apply the Metropolis algorithm to evaluate the sum over the lengths. Sampling over just the lengths is easier as the burn-in period is shorter and the Markov chain can be sampled more frequently since the auto-correlation is lower. Details are provided in the next section.
For this system, the observables we are interested in are detailed in Section 5.2. The geometric observables can be calculated both by sampling the lengths and fields as well as by analytically summing over the fields and then sampling the lengths. We compute the 4-volume of the spinfoam to investigate the interplay between the scalar field amplitude which favours small lengths and the cuboid spinfoam amplitudes, which favours large lengths as α is increased. We find that there exists a range of α where the opposing effects of the scalar field amplitude and the spinfoam amplitude balance resulting in finite values of the 4-volume. These results are detailed in Section 6.1.
The most important observable of the fields is the 2-point correlation function. We compute the correlation function φ a φ b by averaging the product of the field values over all pairs of vertices a and b. We find that when plotted as a function of the geodesic distance between vertices a and b this correlation function falls off exponentially with distance, much like the 2-point correlation function on a fixed lattice, for the range of α where the 4-volume remains finite.

An alternative algorithm
In addition to the algorithm described above, in which both lengths and fields are treated by Monte Carlo techniques, an alternative algorithm is possible. Independent of the shape of the underlying spin foam, the action of the scalar field is always of the form φ a K ab φ b . Thus, if the length are fixed, the integrals over the fields are Gaussian and can be straightforwardly performed. Hence, we can rewrite the partition of the coupled system purely as an integral over the lengths: In this representation, we would use Z −1 1 ) as the probability distribution to sample the lengths. The advantage is apparent: this distribution includes the full dynamics of the scalar field and thus its impact on the gravitational theory. This is particularly helpful for parameter ranges for which the scalar field is difficult to sample, namely small and large α. Moreover, sampling just over the lengths is less sensitive to how new length configurations are proposed, auto-correlations decay more rapidly and the burn-in period is shorter as well. Therefore, this algorithm is well suited to study geometric observables, such as the expectation values of lengths or the total volume. However, computing observables of the field is more involved. Of course, the 2-point partition function for fixed length is simply given by the inverse matrix of K, such that we obtain for its expectation value: This expression is clearly different than for geometric observables, since the factor √ det K −1 is missing. To sample the lengths for this observable, we could use Z −1 v∈ σ4Â v ({l i } i∈σ4 ). However, this is actually not a probability distribution, since the Z is the partition function of both matter and spin foams. Moreover, we would actually only use the spin foam amplitude to sample length without considering the effect of the scalar field, which leads to dramatically different results.
Alternatively, we can attempt to use the previous probability distribution by expanding the integrand with √ det K √ det K = 1; the denominator is part of the probability distribution, while we absorb the numerator into the observable. Unfortunately, this idea is not numerically feasible. Generically det K is a large quantity, even for relatively small lattice sizes. Hence, it completely dominates the expectation value and convergence of the results is slow. Therefore, we refrain from using this representation for Monte Carlo simulations.

Results
The simulations were performed at the Ara-Cluster at FSU Jena. For each set of parameters (α, M and N ), we randomly generated an initial configuration of fields and lengths followed up by a burn-in of 10 8 Monte Carlo steps. This was more than sufficient for the system to thermalize. Then, depending on the lattice size we sampled the Markov chain at varying intervals to ensure that the samples were not correlated. We chose the following interval lengths by explicitly measuring the autocorrelation lengths, i.e. the number of steps until two consecutive samples are uncorrelated and can be considered independent: • N = 81: 5 · 10 4 Monte Carlo steps, • N = 256: 10 5 Monte Carlo steps, • N = 625: 2 · 10 5 Monte Carlo steps.
Due to the increase in auto-correlation length with lattice size and the increase in the time per iteration of the algorithm, we reduced the number of samples taken for lattice size N = 625 to 2000 while 5000 samples were generated for lattice sizes N = 81 and N = 256.
The algorithm and the data generated for and presented in this article are freely available. The algorithm can be found on Github 9 , whereas the data is available on Zenodo 10 .

Geometric observables
In section 5.1 we obtained some qualitative insights for the scenario when all lengths are constrained to be equal. To check whether these qualitative insights also apply to the full model, expectation values of coarse geometric observables provide a good point of comparison. A suitable observable is the total four dimensional volume of the universe and the square root of its variance (normalized with respect to the volume expectation value). In the figs. 7 and 8 we plot both expectation values as a function of the parameter α for lattices of size N = 81, N = 256 and N = 625 for different values of the scalar field mass M .
We observe the same qualitative behavior in all cases, independent of the mass and lattice size. For α < 0.63 the total volume is small, close to the volume given by the minimal value of the lengths prescribed by the lower cut-off. On the other hand, for α > 0.67 the volume is large and close to the volume for largest allowed edge lengths (the upper cut-off). These two regions are clearly cut-off dependent, such that our simulations do not reflect the actual dynamics in these regions. But they suggest that for small α the lattice lengths are as small as possible, while they are as large as possible for large α consistent with our previous observations. The interesting regime lies between these two extremes: the boundary of the intermediate region is marked by steep increases of the average total volume, from the lower cut-off to a finite volume around α ≈ 0.63 and from a finite volume to the upper cut-off around α ≈ 0.67. The finite average volume is a function of α and the scalar field mass M . Increasing α gives a larger volume, while a larger mass drastically impacts the average volume. Crucially, this region is independent of the choice of upper and lower cut-offs on the lengths and thus accurately captures the dynamics of the system. Due to the high symmetry of the cuboid model, we thus expect it to look on average like a regular lattice with an emergent lattice scale. To confirm whether this average is a good representative of the geometry, further observables must be examined. Note that these results are in good agreement with our previous considerations of the equilateral system in section 5.1. Moreover, the cut-off independent existence of a regime of finite volume is of note, as this does not exist for cuboid spin foams without matter [40].  To learn more about the typical geometry of the spin foam, if such a notion makes sense, we also examine the (square root of the) variance of the volume presents an interesting feature: for smallest and largest volumes as well as the finite plateau region, the variance of the total volume is small. This has different meaning in the different regions: when the total volume is determined by upper or lower cut-off, the variance is small due to the cut-off. E.g. at the lower cut-off, the probability distribution would favor even smaller lengths, however such Monte Carlo steps are forbidden by construction. On the other hand, larger lengths are highly improbable. In contrast, the interpretation in the finite volume region is different: there all length variations are a priori permitted, yet large deviations away from the average configuration are not probable. Thus, we can conclude that the spin foam is peaked upon an emergent regular lattice with small deviations around it.
However, at the boundaries of the plateau, where the system transitions from a finite to a "cut-off" volume, the variance shows pronounced peaks. This is because the system transitions from being peaked at a finite volume to being peaked on the lower / upper cut-off. Therefore, the volume is permitted to vary in some cases over several orders of magnitude. The height and position of the peaks might be cut-off dependent, as we limit the minimal / maximal volume the system is allowed to take. This regime is interesting, as these peaks might indicate a phase transition of the system with the possibility of taking a continuum limit. More investigations are necessary though, e.g. considering larger systems.

Comparing different lattice sizes
In addition to the previous plots where we studied the total volume and its variance for different masses and the same lattice size, we briefly compare these quantities for the same masses and different lattice sizes. In fig. 9 we compare lattice size N = 81 (3 independent edge lengths) and N = 256 (4 independent edge lengths) for a mass M = 5. and lattice N = 81 and N = 625 (5 independent edge lengths) for a mass of M = 0.5.
Qualitatively, there is little difference between different lattice sizes. In the plateau region, the larger lattice size typically has a larger total volume, which indicates that the typical lattice lengths are similar for the same parameters. Moreover, the beginning and end of the plateau region is slightly shifted towards smaller α for larger lattices. This matches our expectations from our regular lattice considerations in section 5.1. We can confirm these observations also by the position of the peaks in the variance plots.

Distribution of edge lengths
From the finite expectation value of the total volume in the plateau region it immediately follows that the spin foam edge lengths are finite as well. Moreover, the low variance (compared to the volume expectation value) implies that the deviations are rather small. These results are compatible with our analytic considerations in section 5.1, where we fixed all edge lengths to be the same. Due to the high degree of symmetry of the system a finite total volume suggests that effectively the coupled system behaves like a regular lattice, whose lattice length is emergent and determined by the parameters of the model (α and M ), plus small deviations. However, the total volume is a coarse observable and might average out other features. To complement our understanding of the system, we consider the distribution of lengths.
To do so, we plot histograms of the length distribution for different values of α (in or close to the plateau region), M and the lattice size in figs. 10, 11 and 12. In these histogram we show how often all lengths, i.e. in all space-time directions, were sampled in an interval. From a qualitative perspective, the plots are similar for any lattice size, mass M and α: we see a (fairly) smooth distribution with a clearly defined peak at a finite length; some plots look more jagged, which is due to the finite sample size (in particular for lattice size 5). The peaks are asymmetric, with a sharper decline to shorter lengths and a long tail towards larger lengths in good agreement with our analytical examinations in section 5.1. In short, this distribution shows that in the plateau region all lengths are sharply peaked around the same value, which confirms our expectation that the system is peaked on an emergent regular lattice. Only small deviations are permitted, such that we are confident to characterize the effective spin foam geometry as a regular lattice plus perturbations. In the following, we will discuss how this effective geometry changes as we change the parameters of the model. While all plots look fairly similar, there are strong quantitative differences when changing the parameters. Most striking are changes in the mass M : increasing it by an order of magnitude moves the peak of the length distribution one order of magnitude lower for all lattice sizes and values of α (in the plateau region), and suggests a linear relation between the mass and the position of the peak. However, while the absolute variance of lengths decreases, relative to the peak of the distribution the variance appears to be the same. In this sense, a larger mass "shrinks" the distribution inversely proportional to the mass increase.
In contrast, increasing the value of α, while keeping the mass and lattice size fixed, moves the peak of the length distribution to larger lengths. Moreover, while we always see a sharp peak in the plateau region, the tail towards large lengths becomes longer as we increase α suggesting a larger variance. This   is also confirmed by the variance of the total volume. However, compared to our analytical study of superimposed regular lattices, the tail appears to be slightly shorter. However this impression could be misleading, since these lengths are less probable and might require larger sample sizes to properly reflect the size of these tails.
Finally, we consider the changes of the lengths distribution for different lattice sizes. Firstly, if we keep the mass M and α fixed, the peak of the distribution appears to be more or less constant, that is the typical edge length is mostly determined by M and α and not by the lattice size. Clearly, it is necessary to go to larger lattice sizes to confirm that. Additionally, it seems that the distribution becomes more sharply peaked for larger lattice sizes, in particular the tail is shorter. Again, though more samples might be necessary to properly represent the distribution at such low probability, there is another reason why larger lattice sizes can be more restrictive.
Due to the high degree of symmetry of quantum cuboids, the entire spin foam is characterized by 4N edge lengths, where N is the lattice size. If we vary a single edge length in the spin foam, we affect N 3 hypercuboids and many scalar field interactions via the change of edge lengths directly and also the changes in (dual) volumes. Hence, for larger lattice sizes edge length changes affect more interactions and degrees of freedom compared to smaller lattices, which might lead to more disfavored proposals. This could be an explanation why for larger lattice sizes the tails in the length distribution might be shorter. Obviously, this is a peculiarity of the cuboid spin foam model and does not hold for general triangulations. It will be interesting to see whether a similar effect also occurs in such cases.
In summary, the sharply peaked length distributions suggest that in the plateau region, the spin foam on average looks like a regular lattice with average lattice spacing determined mainly by the scalar field mass M and the parameter α. Deviations around this regular shape do occur. In the next section we investigate the 2-point correlation function of the scalar field on this spin foam, and whether it is sensitive to the fluctuations around the regular lattice.

Matter Observables
The relevant matter observables of the coupled system are the 2-point correlation functions φ a φ b of the scalar field, where a and b denote vertices in the cubulation. Technically, this is not a pure matter observable: in addition to the correlation of the fields, we also measure their distance given by the specific spin foam configuration. This is crucial for the physical interpretation of the correlator, since the labelling of the vertices by itself has no physical significance and is merely a choice of coordinates. Of course, the combinatorics of the lattice plays a role, e.g. whether two vertices are neighbours and thus the associated fields interact directly. However, we do not track this information, as the distance appears to be the determining factor for the correlators.
In the following, we discuss the 2-point correlation functions in the "plateau" region for different lattice sizes. Investigating the correlators outside this region has little merit, for different reasons. Primarily, the results depend on the chosen upper and lower cut-offs and therefore cannot be trusted. Even if the cut-offs could be removed, sampling fields for infinitely small or large lattices is not possible or reasonable. Not to mention that the semi-classical approximation breaks down for spin foam configurations with tiny lengths.
In fig. 13, 14 and 15 we plot the samples of the 2-point correlators φ a φ b for all pairs of vertices a, b over the distance that vertices a and b are apart (in the respective sample), while also dropping the labels. That way we consider scalar field correlations as a function of distance. We define the length between two vertices as the length of the shortest path between them, i.e. we sum up the edge lengths along this path. Note that there are several paths connecting two vertices due to periodic boundary conditions. Thus we consider field correlations as a function of distance between the fields, similar to lattice field theory, yet these distances are dynamical and determined by the gravitational field encoded in the spin foam. Hence, we are measuring a relational observable [52,53,54]. Following our discussion in the previous section, let us first discuss the distances of vertices in the correlation plots. For lattice sizes N = 81 and N = 625, see figs. 13 and 15 respectively, we see that the correlations are are measured at almost equidistant distances with barely any spread, which is difficult to resolve at the plotted distances. This observation is in agreement with the distribution of lengths, and shows that effectively the spin foam behaves like a regular lattice. Comparatively, for lattice size N = 256 in fig. 14 the correlations appear in "bunches", but in a distinct pattern: closest and furthest neighbours appear only in a single line, the intermediate one in roughly three lines. Since this consistently occurs in all simulations for lattices with 4 independent edge lengths but not for lattices with 3 and 5 independent edge lengths, and we find qualitatively similar length distributions for all three lattice sizes, this is most likely a peculiarity of lattices with an even number of independent edge lengths. Due to the cuboid symmetry and periodic boundary conditions, there exist 16 different paths (of different lengths) between two vertices, essentially the two directions per dimension. On odd lattices, there always exists a combinatorially shortest path, which immediately translates to (almost) regular lattices. On even lattices, however, certain pairs of vertices can have multiple paths of the same combinatorial distance between them. In these cases, the fluctuations around regular lattices become crucial and thus might allow for more deviations of the distance. Hence, the pairs of vertices which have an intermediate distance between them show the largest spread in their distance, while furthest away pairs are sharply peaked. Unfortunately, we cannot confirm this for a larger even lattice, e.g. lattices with 6 independent edge lengths per dimension, as the numerical costs grow drastically. These results furthermore support our hypothesis that on average, the spin foam appears as a regular lattice with fluctuations around this regular shape.
For the correlations themselves, we make two main observations. Firstly, the correlations fall off quickly and are typically below 1/4 already for the smallest non-vanishing distance. This is mainly due to magnitude of this distance, which is dynamical. As discussed above, while smaller edge lengths are permitted by the spin foam model, they are highly improbable and thus it is extremely difficult to measure the correlations at small distances. On the other hand, because these configurations are so unlikely, we can assume that the correlations of the scalar field for smaller non-vanishing distances are highly suppressed, but more research is necessary to explore this question further. The second observation is that there is significant spread in the correlations also at large lengths. At such distances, the system actually is difficult to sample as only extremely small values of the scalar field are viable. However, due to these large deviations on the correlations and finite sample size, the correlations might on average appear to be larger than expected for such distances in ordinary lattice field theory. This might affect the estimates of the mass from the exponential decline of the correlators.
Qualitatively, the results look strikingly similar for different masses M , but careful attention is   necessary: For all sets of parameters, we qualitatively observe that correlations exponentially decay as the distance between vertices grows. This observation is true for any lattice size, mass of the scalar field M and value of α (in the plateau region). However, α and M play starkly different roles. α is a parameter that only appears in the spin foam amplitude and it is directly related to the lengths: larger α implies a larger probability of larger lengths, meaning that the distance between vertices a and b increases, which in turn leads to lower correlations. Conversely, a larger mass implies a shorter correlation length, however at first sight the correlations do appear to decay similarly to the cases with lower mass. This is due to the change in length; a larger mass favours smaller lengths, such that the vertices are on average closer together. Thus, the measured correlation length is smaller, due to larger M , yet the fields are similarly correlated since the vertices are effectively closer. Hence, the scalar field and in particular the value of its mass has a strong impact on the geometry of space-time. On the other hand, the 2-point correlation function on average behaves like the 2-point correlation function of a scalar field theory of the same mass on a regular lattice, whose lattice spacing agrees with the average lattice spacing of the spin foam. This also holds for any lattice size, where essentially a larger lattice size adds pairs of vertices to the system that are further apart. To quantify our impression that the 2-point correlator of the scalar field is barely affected by the dynamical spin foam, we compute the correlation length, i.e. the inverse of the "effective" mass, from the exponential decay of the correlator. In scalar lattice field theory, the 2-point correlation function exponentially decays with the distance between the vertices, where the decay rate is called the correlation length. For a free, massive scalar field, this correlation length is equal to the inverse mass of the scalar field; a larger mass implies a shorter correlation length. In our setup, we can thus disentangle changes of the correlations and lengths by fitting the exponential decay of the correlation function and infer an effective mass of the scalar field coupled to the spin foam. Note that this is also a coarse, non-local observable derived from the entire system.
In figs. 16 and 17 we plot the measured effective mass for a fixed mass M of the scalar field as a function of the parameter α and different lattice sizes. First of all, we see that in the plateau region the effective mass is at a similar order of magnitude as the field mass M , yet it changes with α. For α around 0.64, we see a good agreement between the fitted and actual mass. However, the effective mass decreases as we increase α, to around half the value of the field mass close to the end of the plateau. Outside the plateau region, the mass estimates from the fits cannot be trusted, because either the lattice is close to the lower or upper cut-off respectively. However, the decline of the mass in the plateau region towards larger α remains to be understood.
The decline of the effective / fitted mass is universal for different lattice sizes and masses. Since this decrease is tied to α, the cause might be related to changes in the spin foam, e.g. the length distribution. For larger α, we observe a broader spread in the length distribution with a longer tail towards larger lengths. However, it is not clear how this would translate into correlations that die off later and thus result into a smaller effective mass. Moreover, for larger lattice sizes this broadening of the length distribution is less pronounced, yet the decline in mass is comparable. Instead a different explanation appears more likely.
As discussed above, we observe a significant spread of correlations also for large distances, even though the correlations are exponentially suppressed. This is due to the fact that scalar lattice field theory is difficult to sample in those regions. Due to this systematic issue of Monte Carlo simulations, the correlations might appear on average larger than they should be, from which we obtain a larger correlation length / smaller effective mass from the fit to the measured 2-point correlation function. To our best knowledge, this is the most plausible explanation for the observed effect, which makes it unlikely that this difference between mass parameter M and effective mass is a physical effect. Certainly, this question for the cuboids should be investigated in more detail. Furthermore, it will be interesting to see whether this result still holds for more general spin foams, which allow for curved geometries.
To summarize, the observables of the scalar field coupled to cuboid spin foams are remarkably similar to studying a scalar field on a regular lattice, giving further evidence to our hypothesis that our model can be effectively regarded as a scalar field on a regular lattice with dynamical lattice spacing. Changing the mass of the scalar field mostly affects the spin foam, while the correlation function remains similar to one on a regular lattice. Estimating the scalar field mass from the correlation function does reveal a decline with α, yet it seems unlikely that this is a physical effect and rather rooted in the difficulty of sampling the scalar field for large lattice lengths.

Summary and discussion
In this article we have defined the coupling of free, massive scalar field theory to a 4d (semi-classical) spin foam model and for the first time computed expectation values of this coupled matter gravity system in the restricted case of cuboid spin foams [40]. The idea of the coupling mechanism is to consider a spin foam configuration as an irregular lattice on which we define lattice field theory via discrete exterior calculus [55,56,58]. Then, we sum over spin foam configurations and thus superimpose lattice field theories weighted by spin foam amplitudes, akin to the ideas developed in [14]. Our ansatz is not restricted to cuboids and can be straightforwardly adapted to other semi-classical approaches like frusta spin foams [48] or effective spin foams [35] 11 . The particular advantage of cuboid spin foams is that its amplitudes are non-oscillatory, even in the quantum regime [57], such that we can use Markov Chain Monte Carlo techniques to compute observables of the matter spin foam system.
In this paper, we have studied two types of observables, geometric ones, like the total volume and edge lengths distribution, and matter ones, like the 2-point correlation function. From these observables, we identify two counteracting mechanisms of the coupled matter gravity system: the scalar field leads to a polynomial suppression of large lengths, derived from integrating out the scalar fields, which is more pronounced for larger mass. On the other hand, the parameter α in the spin foam face amplitude favours large lengths as it is increased. From the geometric observables, we then identify three regimes: for small α, the lengths of the system are as small as possible, while for large α the lengths are as large as possible. These two regimes are connected by an extended region in parameter space at which the total volume and the edge lengths are finite, and the system on average / effectively looks like a regular lattice. The average lattice size is a function of the scalar field mass M and α and appears to be independent of the lattice size. Note that this lattice size is dynamical and not related to the fundamental discreteness of area variables in spin foams; due to the semi-classical regime we have assumed the lengths variables to be continuous. Additionally, we studied the 2-point correlation function of the scalar field, which effectively behaves as if the scalar field is defined on a regular lattice. While our estimates show some deviations in the effective mass / correlation lengths, we attribute those to the difficulty of measuring correlations for large lattice distances. Interestingly, changing the mass of the scalar field mostly appears to affect the spin foam, whose average edge lengths change inversely proportional to the change of mass.
Our results have several interesting implications for spin foam quantum gravity. Primarily, the effective emergence of scalar lattice field theory defined on a dynamical regular lattice is first evidence that spin foam quantum gravity coupled to matter might possess a regime that effectively describes quantum field theory defined on a fixed background space-time. Additionally, the effective lattice spacing is also dynamical and given as a function of the parameters of the theory, in particular the scalar field mass. This could hint towards a new mechanism as to how new intermediate scales can emerge in spin foams / loop quantum gravity coupled to matter. It furthermore raises the questions as to how other types of matter will influence the spin foam, and whether there are limits to how many matter fields we can couple and still recover a quantum field theory from our coupled system. This is indeed a question raised in many approaches to quantum gravity, e.g. asymptotic safety [26], in particular whether the standard model of particle physics is compatible and whether there is a mechanism that limits the amount of matter fields one can couple to quantum gravity.
Additionally, in this work we gain first insights into how to extract physical observables from matterspin foam systems. In absence of a background space-time, no overall scale like a lattice spacing is a priori available to define a correlation lengths. One can resort to using combinatorial information or simply the labels of vertices, but this is a coordinate distance instead of a physical one. Instead we measure both the distance (here path distance) between two vertices as well as the correlation of the scalar field located at these vertices and drop their labels. That way we store physical, diffeomorphism invariant information, namely correlations in relation to the geodesic distance between two points, in the spirit of relational observables [52,53,54]. Note however, that in our Monte Carlo approach, we cannot probe the correlations at arbitrary distances, only at probable distances. In the review on causal dynamical triangulation [27] a scalar field 2-point correlation function is discussed that defines the correlator at arbitrary geodesic distances.
At this stage it is not clear how to define a continuum / refinement limit [42,44] of the coupled spin foam matter system. From the observables we studied in this work, none appears to be an order parameter indicating a second order phase transition. While the variance of the total volume shows a divergent behavior at the edges of the plateau region, these divergences do not grow for larger lattice sizes. Moreover, it would be necessary to evaluate in more detail whether the variance is dependent on the chosen cut-off. Beyond the conclusions we can draw from our current work, renormalization / coarse graining of the matter gravity system, similar to [84], is indispensable to address this question and to evaluate what influence various ambiguities might have on the dynamics, e.g. the definition of the scalar field theory. In our concrete case, we could attempt to define a renormalization group flow where the system defined on a finer lattice effectively looks like the refinement of a coarser system, while keeping the 2-point correlation function of the scalar field fixed. Essentially this amounts to matching the average total volume by tuning the parameter α and keeping the mass M fixed. In the plateau region, where the 2-point correlation function and correlation length are well understood, this strategy of identifying observables will likely only work for a few coarse graining steps. On the other hand, on the edges of the plateau region, where the lengths start to fluctuate, it might be possible to relate the average geometric observables across different lattice sizes by tuning α. Further work is necessary to explore this idea, in particular how the scalar field correlator behaves for such parameters on larger lattices.
While the results of this article a encouraging and intriguing, it is clear that many assumptions are made that must be lifted in order to check whether they also hold in physically more relevant situations. On the side of the spin foam this requires us to go beyond quantum cuboids and eventually beyond semi-classical models. A suitable first step could be to use effective spin foams [35] defined on a triangulation, which is numerically more accessible than the full quantum theory and allows for vastly more configurations with non-vanishing curvature and non-metricity. However, due to the oscillatory nature of the amplitude, Monte Carlo methods will no longer be applicable and other numerical methods must be explored. A first attempt might be to consider a scalar field coupled to effective spin foams expanded around a flat triangulation of hypercubes similar to the work [85]. Extending the matter coupling to the quantum regime of spin foams is another interesting challenge that must be addressed, as well as coupling other types of matter and studying their dynamics, e.g. gauge field [14] and fermions [18,19].