Analytic Analysis of Irregular Discrete Universes

In this work we investigate the dynamics of cosmological models with spherical topology containing up to 600 Schwarzschild black holes arranged in an irregular manner. We solve the field equations by tessellating the 3-sphere into eight identical cells, each having a single edge which is shared by all cells. The shared edge is enforced to be locally rotationally symmetric (LRS), thereby allowing for solving the dynamics to high accuracy along this edge. Each cell will then carry an identical (up to parity) configuration which can however have an arbitrarily random distribution. The dynamics of such models is compared to that of previous works on regularly distributed black holes as well as with the standard isotropic dust models of the FLRW type. The irregular models are shown to have richer dynamics than that of the regular models. The randomization of the distribution of the black holes is done both without bias and also with a certain clustering bias. The geometry of the initial configuration of our models is shown to be qualitatively different from the regular case in the way it approaches the isotropic model

1 Introduction matter distribution in the universe is discrete and therefore strongly locally inhomogeneous, it is natural to ask how this affects the overall expansion rate. The question is to what extent the dynamics of the real universe is described by the homogeneous fluid FLRW models. Due to the nonlinear nature of the Einstein equations this is a difficult problem often referred to as the averaging problem in general relativity [1,2]. It has been extensively discussed in the past but there is still no consensus about its solution (see e.g. [3]). A new avenue to address this issue was recently opened by considering a universe filled with Schwarzschild black holes as sources [4,5,8,9,12,21]. The idea is to first specify the sources as initial data on a time-symmetric spatial hypersurface with spherical topology. It is possible to do this using a finite number of Schwarzschild black holes as an exact solution of the Einstein constraints. This configuration of initial data can then in principle be evolved using the Einstein evolution equations. While the general evolution can only be done numerically, some subsets with sufficiently high symmetry can be evolved exactly or approximately by analytical methods. The line of investigation is still in its infancy with no clearcut answer to the averaging issue. However, we believe that this approach has a promising potential in this respect. It may also illuminate other related issues, one being the role of interaction energies in many body systems.
In previous work on such cosmologies the focus has been on models with a regular matter configuration of the 4-polytope type. This kind of configuration has been used to simplify the analysis. The regular 4-polytopes are 4-dimensional analogues of the (3-dimensional) Platonic solids. There are six different regular 4-polytopes with 5,8,16,24,120 or 600 cells. Previous work on such regular models placed identical sources at the cell centres. In this way the regularity of the model could be preserved. The edges of the cells in such a model have a high degree of symmetry in the sense that the gravitational field is locally rotationally symmetric (LRS) about the edge. The edge is then said to be an LRS curve [8]. This symmetry will allow an analytic analysis of the evolution up to some time away from the initial data hypersurface.
Our primary goal in this work is to relax the obviously unrealistic condition of regularity to see how this will affect the curvature structure and also the evolution. We wish to construct a model in which the sources can be distributed randomly while still preserving at least one LRS curve. It would not be possible to do this with a completely random distribution. However, the existence of LRS curves depends only on discrete symmetries. Each reflection symmetry has an associated symmetry surface (the mirror). If three or more such symmetry surfaces have a common intersection curve, then that curve is necessarily LRS [9]. Consider the regular 4-polytope with 16 cells as an example and fill one of its cells with a random distribution of sources. At an edge of that cell there are in total four cells that meet. One is the cell we have already filled and then there are three more cells meeting at that edge. We can now fill the other three cells by reflecting the given sources in the adjoining faces going around the edge. The newly filled cells can then be used as seeds to fill the rest of the polytope. Since there is an even number of cells at the edges, this procedure gives an unambiguous distribution of sources for the whole polytope. Each cell is then a mirror image of all its neighbours. However the distribution inside the cells can be arbitrarily irregular.
While the above construction works well it cannot be used for any other regular 4-polytope. The reason is that all the other polytopes have an odd number of cells meeting at the edges. In particular, it cannot be used for the 8-cell polytope which is the one that has been most extensively investigated in the past. For that reason we will use another model with eight identical cells which is also a tessellation of the 3-sphere but which is not of the regular polytope type. Instead, the cells can be characterized as solid lenses. All the cells meet at a great circle. If the number of cells is of the form n = 6 + 2k with k a non-negative integer, then the common great circle retains the LRS property when the above construction of the source configuration is used. The choice k = 1 gives us an 8-cell model which contains the desired LRS curve.

The Initial Metric of the DI Model
Following [4] (see also [5]) we consider an initial instantaneously static spatial hypersurface with metric whereĥ is the metric of a 3-sphere and ψ is a conformal factor to be specified below. In hyperspherical coordinates the metric has the form with coordinates in the ranges 0 < χ < π, 0 < θ < π and 0 ≤ φ < 2π. For some purposes it is useful to embed the 3-sphere in a 4-dimensional Euclidean space with Cartesian coordinates (w, x, y, z). The relations to the hyperspherical coordinates are w = cos χ x = sin χ sin θ cos φ y = sin χ sin θ sin φ z = sin χ cos θ .
The point where χ = 0 is referred to as the north pole of the hypersphere. We will proceed as in [4] and solve the Gauss-Codazzi equations: where K ij is the extrinsic curvature of the 3-space, K is its trace and R is the Ricci scalar of the (physical) 3-space. The vertical line in equation (3) denotes the (spatial) covariant derivative of the metric h. Since the hypersurface is instantaneously static, the extrinsic curvature vanishes, K ij = 0, and the Gauss-Codazzi equations are then reduced to the very simple form For the metric (1), this condition reduces to the Helmholtz equation on the 3-spherê whereR is the 3-sphere Ricci scalar and∇ 2 is the 3-sphere Laplacian. Since the Helmholtz equation is linear in ψ we can form a superposition of single black hole solutions. A conformal factor corresponding to N Schwarzschild masses on the background 3-sphere can be defined as where n and n k are unit vectors in the embedding space. The location of each source is provided by the vector n k while n = (cos χ, sin χ sin θ cos φ, sin χ sin θ sin φ, sin χ cos θ) .
The full hypersurface metric is then obtained by inserting the conformal factor (6) in (1).

Locally Rotationally Symmetric Curves
Locally Rotationally Symmetric (LRS) spacetimes have been defined and investigated rather extensively by several authors following [10]. A spacetime is LRS if the tangent space of every point has a direction about which there is no preferred perpendicular direction defined by the geometry. Here we consider spacetimes which are not LRS everywhere but contain curves of LRS points. Our DI models are constructed with symmetry surfaces whose intersections are LRS curves. Such curves are useful since they allow a significant simplification of the dynamics. Specifically, the dynamics for points on the LRS curves reduces to a system of ODEs driven by a non-autonomous term which can be evaluted (in principle) by a recursion procedure [9]. In practice, it has only been possible to evaluate the recursion a few steps before it becomes unwieldy. Nevertheless, even the undriven ODE system gives results that match current numerical work to an accuracy of about one percent [12]. Suppose that there are n surfaces which intersect along a curve and that the geometry admits reflection symmetry about each of the surfaces. Now arrange an orthonormal frame (e 1 , e 2 , e 3 ) with e 1 parallel to the curve. There will then be a discrete set of rotations about the curve such that the frame vectors are transformed asẽ 1 = e 1 e 2 = e 2 cos φ q − e 3 sin φ q (8) e 3 = e 3 cos φ q + e 2 sin φ q while leaving invariant all tensors picked out by the geometry The rotation angles are given by φ q = 2πq/n where n ≥ 2 is the number of intersecting surfaces and q = 1, . . . , n − 1. For n ≥ 3 it can be shown that curves associated with such a discrete symmetry group have the LRS property [8] while for n = 2 this does not follow (see discussion in [5] and [9]). Vectors and tensors satisfying (9) have the following properties when evaluated along LRS curves: 1. For a vector T α evaluated in the orthonormal frame we have T 2 = T 3 = 0 while T 1 remains arbitrary. Therefore vectors respecting the discrete symmetry group are always parallel to the LRS curve.

2.
A rank-2 tensor with orthonormal components T αβ must be diagonal along an LRS curve with components satisfying T 22 = T 33 .

A Piecewise Statistically Uniform Distribution of Sources in a 3-Sphere Universe
We wish to construct a universe in which the sources are distributed as randomly as possible but still contains at least one LRS curve. We can then use the existence of the LRS curve to get estimates of the dynamics of the model. To satisfy the LRS requirement we need a model which has a symmetry subgroup containing at least three distinct reflection isometries. The symmetry surfaces in such a model leads to a lattice of cells which are identical up to reflections. Every cell face is then part of a symmetry surface for a reflection. For simplicity, consider a lattice with identical cells. One could then try to define a piecewise random source configuration as follows. First distribute N sources randomly in one of the cells. Then for each face of that cell reflect the given source distribution to a neighbouring cell using the reflection symmetry corresponding to that face. Further reflections can then be used to distribute the given source configuration to all the cells in the model. For such a construction to work, the relevant symmetry subgroup of the model must be preserved when all the cells have been filled. Otherwise, the LRS property will not in general be preserved. It is straightforward to see that a necessary condition for this is that the number of cells meeting at an edge must be even. The number in question is the third Schläfli number. Of the regular models in [4], only the 16-cell satisfies the requirement with four cells meeting at the edges. It can be shown explicitly that the above construction does work for the 16-cell. However, we prefer to use a different model that contains only eight cells and for which there is a simple algorithm that can be used repeatedly to populate the cells.

Tessellating the 3-Sphere
Similarly to [4], we consider a 3-sphere representing a momentarily static hypersurface in the universe. It is tessellated into a set of empty and identical cells. However, unlike in [4], the tessellation in our case will not correspond to any of the regular polytopes but will instead be taken to resemble hyperspherical versions of "orange pieces". Suppose we have an orange and intend to share it by cutting it into a number of wedge pieces. Typically one might make the cuts where the black curves are in figure 1 which would result in several manageable pieces. If one ensures that each piece is of the same size, the number of pieces will be determined by the "equatorial width" ∆φ.
In terms of a general 2-sphere instead of an orange, the surface of each orange piece corresponds to a cell, while the knife cuts are cell boundaries. The cell boundaries are located at fixed φ-coordinates and the angular distance to the nearest cell boundary is given by ∆φ, which must be equal for all cells.
In addition, all cell boundaries intersect at two points -the north and south pole.
Keeping this image in mind, we now make our move to the 3-sphere. In this transition surfaces become volumes, curves become surfaces and points become curves. Thus our "orange" becomes a 3-sphere, with 3-dimensional cells, each bounded by two 2-dimensional spherical surface segments. These cell boundaries intersect along a curve C LRS (instead of at merely two points) which is a great circle passing through the 3-sphere's north and south poles located at χ = 0 and χ = π respectively. In Cartesian coordinates (w, x, y, z) the curve C LRS is parametrized by The cell boundaries ∂V q are a series of spherical surface segments which enclose the 3-dimensional cells V q . Using hyperspherical coordinates (χ, θ, φ), the cell boundaries ∂V q are given by The integer p corresponds to the number of intersecting surfaces (cell boundaries). The form of each cell can be viewed as similar to a biconvex spherical lens (cf. [20]).

Filling the cells
After tessellating the 3-sphere as described above, the resulting cells are to be filled with a number of Schwarzschild black holes in a manner which ensures that the curve C LRS actually remains LRS. We do this by first designating one of the cells as a starting cell (this choice is completely arbitrary since all cells are identical). Then fill the starting cell with one or more points (corresponding to the location of sources, specifically black holes) in any desirable manner. We choose to distribute the points randomly in the starting cell by using an algorithm invented by Marsaglia [6]. The algorithm is designed to generate random uniform distributions on the 3-sphere.
1. First select four points V 1 , V 2 , V 3 and V 4 randomly from the interval [−1, 1] and ensure that 2. Then the following point is picked randomly from a uniform distribution on the unit 3-sphere: 3. If it lies outside the starting cell -discard it. Otherwise -keep the point. Initially we only want points inside our starting cell.
4. Repeat steps 1-3 until a desirable number of points lie inside the starting cell. After the starting cell has been filled with a random configuration of points, the remaining cells are filled with a point configuration such that every cell becomes the mirror image of its neighbour. This is most easily done in the stereographic projection of the 3-sphere, where the cell boundaries are flat surfaces instead of spherical -see figure 2. Since the cell boundaries are symmetry surfaces, the configuration of the starting cell can easily be mirrored into the neighbouring cells. Once every cell is a mirror image of its neighbour, an inverse of the stereographic projection can be made to return to the original 3-sphere. The cells will still remain mirror images of each other since the projection is conformal. Furthermore since two reflections correspond to a rotation there will also be a discrete rotational symmetry around the common cell boundary C LRS [7].

The Specifications of our Model
We have chosen to investigate the dynamics of a DI orange piece model consisting of 8 cells corresponding to p = 4. This is a model with few cells which is also comparatively simple to handle technically when defining the source distributions. In addition, having 8 cells allow us to compare with regular models in [8] in which the total number of masses investigated were 8, 16, 24, 120 or 600 (corresponding to 1, 2, 3, 15 or 75 masses per cell respectively).

Finding the Proper Mass of Schwarzschild Black Holes in the DI Model
In paper [4], two types of mass were discussed, viz. the effective mass and the proper mass. The effective mass m eff corresponds to the mass parameters m k in (6). It is the mass the source would have if it were alone in an asymptotically flat universe. However when multiple masses are present, m eff will also contain interaction energies (which are negative). In this paper, the effective mass is set to a constant. The reasoning behind this is simplicity -it is easier to calculate proper mass from a fixed  Figure 3: Mass ratios m p /m eff plotted against the number of masses n obtained from more than 200 randomly generated configurations. The red curve shows the median mass ratios and the error bars indicate the quintiles. For 8 sources, there is the possibility to have one source at the center of each cell. This is the closest our DI model can approximate the regular 8-source model in [4] and is indicated by the blue dot. The green curve depicts the regular ratios found in the aforementioned paper. effective mass rather than the other way around. The actual value of m eff will then simply correspond to an overall scaling of all the masses. We therefore present the data as the ratio m p /m eff .
The proper mass is defined as in [4] and is analogous to the concept of bare mass in [14]. Physically it can interpreted as the "local" mass as measured in a vicinity of the black hole. The difference between the two mass concepts can be interpreted as due to the interaction energies between all the sources in the configuration. The proper mass of a Schwarzschild black hole at χ = 0 (north pole) is calculated by expanding the conformal factor (6) where c 1 and c 2 are constants depending on the mass configuration. The proper mass is then given by [19] The proper mass of sources at other positions can be calculated by first rotating the coordinate system (or the mass configuration) to place the north pole at a given source.
Referring to the expansion (13) it is clear that the proper mass depends on the configuration of sources -different configurations will yield different proper mass values. Since any source configuration is randomly generated, a wide range of proper masses can be expected. We proceed by randomly generating many (over 200) configurations and calculating the proper mass of the sources.
However, not all sources in a given configuration need to be accounted for. Any source in one cell and its counterpart in every other cell are equally "massive", because each cell are identical up to a rotation or reflection. It is therefore not necessary to study all proper masses in a given configuration, only those inside a single cell.
We have chosen to characterize the distributions by their median, and upper and lower quintiles (which separates the highest, respectively the lowest, 20% of the data from the remaining 80%). These are then plotted as a function of the number of sources n, as done in figure 3. We can also compare with the mass ratios obtained for the regular models in [4]. Clearly the mass ratios for the DI and the regular model converge for large n [4,8] and become more massive. It therefore appears that the overall interaction energies of the DI model are approaching those of the regular configuration in the limit of large n. 3 Initial Behaviour on the LRS Curve

Initial Length of the LRS Curve
The initial length of the LRS curve is simple to calculate, since the curve is a great circle which can be parametrized by χ in hyperspherical coordinates. Using the metric given by the expressions (1) and (6), the LRS curve length is then given by the integral A factor of two arises since both halves of the LRS curve are identical (due to reflection symmetry) and the above integral only takes half of the curve into account. This can be compared to the spherical matter-dominated FLRW model (with no dark energy present) at the moment of maximum expansion: In this geometry, the length of a great circle with the same parametrizations as mentioned above is where M is the total (proper) mass in the universe. Figure 4 plots the ratio between the LRS curve length in the DI model and the FLRW model as a function of the number of sources. The length in the FLRW model is computed according to equation (17) where M is the total proper mass in the DI model (obtained using the method in section 2). The large dots indicate the median, while the accompanying error bars indicate quintiles. For histograms detailing the length distributions in detail, see [7]. By including data for the regular models from [4], we see that both model types approach the FLRW limit as more sources are added. Our results are so far consistent with [23], where a thorough analysis is made of the continuum limit. However, while the regular models approach the FLRW limit from above, the irregular models approach it from below.
There are thus two different ways to approach the limiting homogenous case presented by the FLRW model.
A possible cause for these two types of behaviours can be found by examining figure 8 in [24]. In that work, a cosmological model with clusters of black holes and a tessellation identical to that in [4] is studied. There the authors also observe an approach to the FLRW limit from above. This observation suggests that the choice of tessellation of the 3-sphere will affect how the FLRW limit is approached.

Initial Curvature Along the LRS Curve
The DI models, being vacuum universes only support Weyl curvature. This is in contrast to FLRW models which have only Ricci curvature as a consequence of their conformally flat nature. The Weyl curvature tensor C µνρσ can be invariantly decomposed in two parts, a gravito-electric part E µν and a gravito-magnetic part H µν . They are defined by the relations where u µ is a 4-velocity field representing an observer family and * C µρνσ u ρ u σ is the dual Weyl tensor. The latter is given by * where η µνρσ is the Levi-Civita tensor. The tensors E µν and H µν are symmetric and traceless by their definition. Their definition also implies the relations E µν u ν = 0 and H µν u ν = 0. This shows that they are spatial tensors lying in the rest spaces of the u µ observers. At the initial momentarily static hypersurface, the gravito-magnetic field is zero while the gravito-electric field equals the traceless part of the intrinsic Ricci 3-curvature of the initial hypersurface, E µν = (3) R µν (see e.g. [8]). We follow the convention [11] and define the following five variables representing the orthonormal components of any symmetric tracefree spatial rank-2 tensor (using the gravito-electric tensor E αβ as an example): However, taking e 1 to be directed along the LRS curve, tracefree symmetric rank-2 tensors are diagonal with E − = 0. Therefore the quantity E + is the only nonzero curvature component along the LRS curve. We will investigate the behaviour of (E + ) 0 along the LRS curve on the initial static hypersurface with the zero subscript indicating evaluation at t = 0.
The LRS curve, as we have defined it, is located at θ = 0, π. However, as is evident from Eq. 1, the determinant of the metric is zero where θ = 0, π implying that the entire curve is located at a coordinate singularity. Therefore, to calculate the curvature, we first rotate the configuration such that the LRS curve instead lies at χ = θ = π/2 and is parametrized by φ ∈ [0, 2π]. Therefore in the ensuing graphs, the position on the LRS curve is given by the hyperspherical coordinate φ. Figures 5 -6 show some examples of how the quantity (E + ) 0 behaves along the LRS curve. It can be seen that in some cases the sign of (E + ) 0 changes. To compare with the regular models of Ref. [8] we can consider the sign of (E + ) 0 along the edges of the cells in that paper, noting that it is negative. Further remarks about the sign distribution can be found at the end of this section. We will also return to discuss a physical interpretation of the sign of (E + ) 0 when treating time evolution in Sect. 4 Now let us take a look at each case individually. From figure 5 we see a unique behaviour exhibited when only 8 sources are present, namely the curvature is constant! This is always the case for any configuration consisting of only 8 sources: the curvature is constant with (E + ) 0 < 0. The reason behind the constant curvature is most likely due to the topology of the 3-sphere.
For 16 and 24 masses, the behaviour of (E + ) 0 becomes more mixed -maxima and minima with occasional changes in the sign of (E + ) 0 are present. Furthermore the behaviour varies with different configurations. However, something which is consistent among different configurations is the number of local maxima and minima. For 16 sources the number of maxima found on the LRS curve is two, which is the number of masses per cell. Similarly for 24 sources the number of maxima is threealso the number of masses per cell. This is evident in figure 5. Due to the relatively small number of sources, the effect of each black hole is distinctly visible.  This average is not unique for a given number of sources n, but depends on the configuration. The plot above presents the result after generating many random configurations, which yields a distribution. As usual in this paper, dots and error bars indicate median and quintiles of this distribution, respectively.
For 120 and 600 masses, the (E + ) 0 behaviour is even more dramatic and irregular. Figure 6 displays positive peaks and negative troughs and Riemann flat points at (E + ) 0 = 0. The latter points are characterized by a vanishing Riemann tensor. The spike-like features are most likely due to a nearby black hole horizon, cf. figure 13 in [8].
Finally we investigate the sign dominance of (E + ) 0 on the LRS curve. The sign of (E + ) 0 at a given point affects the time evolution of any point along the curve. Since there are different signs along the same curve, it follows that different regions might evolve differently [8]. A very rough estimate of the sign dominance is provided by the sign of the mean (E + ) 0 . Figure 7 presents the median and the location of quintiles obtained from (E + ) 0 , as a function of the number of masses. Negative values are clearly the most prevalent, approaching (E + ) 0 =0 when there are many sources in the model.

Equations Governing Dynamics
With the initial conditions studied, we now proceed to evolve them. We restrict attention to the LRS curve for which the evolution can be studied by analytical methods. The evolution is governed by the system of ordinary differential equations [8,9] where C + (t) = − 3 2 n a n b curlH ab .
Here n a is a unit vector along the LRS curve and and D a is the spatial covariant derivative. The quantities a and a ⊥ are the scale factors parallel and perpendicular to the LRS curve. Although the curlH term depends on spatial derivatives it can be computed recursively term by term as a Taylor series in t (for given initial time-symmetric data) using the full Einstein vacuum equations (see [9] for details). Therefore, the system (20)- (22) can be considered as an autonomous ODE system driven by the explicitly time-dependent curlH term. The initial values of the scale factors on the time-symmetric hypersurface (at t = 0) can be taken to be (a ) 0 = (a ⊥ ) 0 = 1 and the time-symmetry impliesȧ =ȧ ⊥ = 0. The initial values of the gravitoelectric field is given for each spatial position by (E + ) 0 = − 3 2 (3) R 11 0 . Considering first the vertex points in the regular models (see [8]) we note that they are points of local spherical symmetry. Then any tensor object picked out by the geometry can only have a scalar part at the vertex. This implies in particular that E ab = H ab = 0. Therefore, the vertex points can be characterized as being Riemann flat. It then follows from Eqs. (20)-(21) that the scale factors are constant for all time, a (t) = a ⊥ (t) = 1.
Next we note that curlH ab is an odd function of t and that the Taylor expansion of (23) takes the form where the first nonzero term is of order t 3 [9]. In [8], Clifton et al. analyzed the autonomous part of the system (20)-(22) (i.e. with C + (t) assumed to vanish identically) for all the regular lattice configurations. By performing a numerical integration of the full Einstein system, Korzyński et al. [12] showed that the error incurred by neglecting the curlH term in (22) is about 1% at t ∼ 3m p when evolving a and a ⊥ at an edge midpoint in an 8-cell lattice. An analytic expression for the first term in (25) was given in [12] as a function of the conformal factor with respect to flat space. Covariant expressions for the first two terms were given in [9], the first one being where (E 2 ) ab = E ac E c b . To analyse our model, we will exploit the analytical solution of the undriven part of the system (20)- (22) given in [8]. In Sect. 6, we will estimate the error incurred by neglection of the driving term. From the physical point of view, this approximation amounts to neglecting the gravitomagnetic effects, i.e. the part of the gravitational interaction which emanates from the relative motion of the masses.

Evolving E + on the LRS Curve
The equations given in section 4.1 allow us to study the curvature, measured by E + , as a function of time. A first example is given in figure 8, where we evolve a model with 8 sources. In that case, since the curvature is constant along the LRS curve (see figure 5), it is sufficient to evolve a single point on the curve. The value of E + increases monotonically until a coordinate singularity appears along the entire LRS curve simultaneously, beyond which our current formalism provides no information.
We have evolved the models with 16 and 24 sources for which initial data was plotted in figure  5. The result is displayed in figure 9. The time and the curvature E + in those figures are scaled by appropriate powers of the median proper mass m p . All points on the LRS curve will eventually encounter a singularity unless initially Riemann flat (i.e. (E + ) 0 = 0). The sign of the curvature does not change along the part of the evolution that we can follow.
Referring to the evolution of models with 8, 16 and 24 sources (see figure 9) it appears that coordinate singularities appears before curvature singularities. This turns out to be the most common behaviour for the models that we have investigated. See figure 10 in which the first appearance of either a coordinate or curvature singularity is illustrated.

Dynamic Contractions and Expansions of the LRS Curve
We now proceed to study the length of the LRS curve (i.e. the LRS circumference) as a function of time. It is calculated by which is essentially equation (15) multiplied by a scale factor. Note that a (t) has an implicit position dependence as described in section 4.1. In particular, depending on the sign of (E + ) 0 , the scale factors will a have qualitatively different time evolution. The LRS circumference is plotted as a function of time in figure 11 for 500 randomly generated models containing 8 sources. The lengths are scaled by the circumference (L F LRW ) 0 of a great circle in a spherical FLRW model with the same total proper mass M p . They are evolved until a singularity appears. The similarity in the form of the curves reflects the scaling of time caused by different initial    values of E + in (28). The density distribution of the curves can be attributed to both the dependence on the position of the (single) source in the cells and to the nonlinear (square root) dependence on (E + ) 0 in (28).
A richer variety of dynamics can be seen for the remaining cases shown in figures 12. Looking at them, we can identify two general qualitative categories. The first category consists of curves which decrease monotonically with time, such as the FLRW case and models with regular configurations. A majority of the DI models exhibits this behaviour. The second category consists of curves which, at some point in time start to increase. This represents a turnaround in which the space in a neighbourhood of the LRS curve undergoes a change from contraction to expansion. This behaviour drastically deviates from the FLRW model and the regular cases. We cannot say if this is representative of the entire universe or if it is a semi-local anomaly.

Hubble and Deceleration Parameters
There are two additional parameters we can measure which enables further comparisons with FLRW cosmology. These are the Hubble and deceleration parameters which we define as where L is the length of the LRS curve. These two parameters are plotted for all the models we consider in figures 13-14. The black dashed curve depicts as usual the spherical FLRW case. In the 8 source model ( figure 13), all curves show the same qualitative behaviour due to the position independent curvature. Even though this case is very different from the FLRW model in terms of size (see figure  4), it is qualitatively very similar.
In the remaining models ( figure 14) the dynamics are more complex. Turning our attention to the deceleration parameter, we can broadly divide the behaviour in two categories. First we have all the curves for which the deceleration parameter is always positive. An example is provided in figure 15a. The curve starts at positive infinity before approaching zero. The FLRW curve and the regular models belong to this category, as well as a majority of the DI models.
The second category is characterized by the presence (at some point in time) of negative deceleration parameters. These curves are usually positive initially. An example is shown in figure 15b. In rare   (31)). These are the most common types. We encounter occasionally deceleration parameters which are always negative, as depicted in c).
cases can we find models with always negative deceleration parameters, as in figure 15c. Regardless, the appearance of negative deceleration parameters is interesting, since in FLRW cosmology, negative values correspond to the presence of some exotic matter such as dark energy [13]. We show however that acceleration can occur across large regions of space in models with only ordinary mass (specifically Schwarzschild black holes) for some matter configurations. Another interesting observation can be discerned from figures 13-14. From the Hubble parameters we see that the fraction of orange curves which terminates below the FLRW curve (black dashed curve) decreases as the source count increases. This provides a crude indication that the mean behaviour of the Hubble parameter becomes less Friedmann-like at late times as we increase the source count (if the mean was Friedmann-like, we would expect an equal number of curves to terminate above and below the FLRW case).
A valid question is whether or not our results would change if we would have included the curlH term. However, we know that the error starts from zero (being of order ∼ t 3 initially), and as demonstrated in [12] and to be illustrated in section 6, the gravitomagnetic effect is very small.

Fitting the Friedmann Equation
The Friedmann equation for a late time flat cosmology consisting of dark energy and ordinary matter can be expressed as [13] ȧ a where Ω M 0 and Ω Λ 0 are the density parameters corresponding to mass and dark energy respectively. Furthermore the relationship between the density parameters are By fitting expression 32 to the blue curve in figure 16, we obtain an indication of what kind of energy content is necessary for an approximate DI-like behaviour. We perform this fit on a single randomly generated DI model with 600 masses exhibiting generic dynamics in the LRS length evolution. The quantities with a 0-subscript in equation (32) indicate "present day" values, obtained through observation/measurement carried out at the "present time" t 0 . We define t 0 to be the time when an observer residing in the DI universe carries out measurements or observations, indicated by the vertical dashed line in figure 16. The observer's time is considered to progress from right to left in the figure. Consequentially, the observer's future is the initially time-symmetric hypersurface at t = 0 and the observer's past is located at large t. The notions of future and past are, with respect to the time axis in the figure, flipped in order to represent an expanding universe. Only parts of the curve where t > t 0 (regions to the right of the dashed line) will be fitted against using equation (32), since the observer possesses only information about the past.
By allowing the present time t 0 to progressively increase, i.e. the observer will be moved further away (temporally) from the initial time-symmetric hypersurface (located at t = 0), there will be less and less portions of the simulated data available for our fit. This will alter the density parameters obtained as a function of t 0 . During the fitting, Ω M is allowed to vary freely but is restricted to non-negative real numbers. Figure 17 depicts the density parameters Ω M 0 and Ω Λ 0 as a function of t 0 . For small t 0 , i.e. close to the initial hypersurface, we have Ω M 0 1 which therefore results in negative Ω Λ 0 . The inset is an enlargement on a region where both parameters are positive simultaneously. Within this region we can have that Ω M 0 = 0.3 (indicated by dashed red line), which corresponds well with current observations (see for instance [16]).

Calculating the curlH-term
In this section we will estimate the curlH term and see how it evolves with time in order to determine the magnitude of the resulting error. The curlH-term in equation (22) can be evaluated recursively in a truncated Taylor series over time. The first and second order terms are zero, permitting us to stick with the third order term only, given by [9] C (3) where the tensor E 2 ab is defined as Figure 18 shows an example of how the magnitude of the ODE expression (see equation (22)) compares to the curlH-term (34) as a function of time for 24 masses on a random point on the LRS curve. A general feature depicted in figure 18 is the initially zero curlH term, which makes the differential equation (22) valid for small t . We can make some statistics of when the ODE and curlH-term are equal in magnitude. The result is shown in figure 19, where we see the median time expressed in terms of the median proper mass m p . We managed to calculate expression (34) for configurations consisting of 8, 16 and 24 masses only. For 120 and 600 masses, the evaluating curl curl E 2 ab t=0 becomes a heavy computational endeavor and is therefore omitted. However a clear pattern can be distinguished: as we increase the number of black holes, the time it takes for the ODE and curlH term to become equal increases. If we compare figure 19 with figure 10, we notice that for 600 sources it appears as if coordinate singularities appear before the curlH term becomes dominant.  7 Evolution with Clusters of Black Holes

Generating Clusters
In our universe, galaxies are not uniformly distributed but form instead wide networks of filaments, clusters and voids. To simulate this behaviour to some extent, we have distributed black holes on the 3-sphere in such a way to encourage the formation of clusters. The effect of discrete mass objects with some form of structure or clustering have already been investigated in for instance [22,23,24]. Our approach however is stochastic and will use an algorithm inspired by the two-point function used in galaxy surveys to quantify clustering [15]: 1. Focusing on one cell, distribute the desired number of black holes in a random uniform way as described in subsection 1.3.
2. Select one of the points randomly, P. The probability of selecting any point is chosen to be equal.
3. Measure the distance between P and every other point inside the same cell. The probability of finding a mass at point P is given by the expression where n is the number of masses, d i is the distance on the 3-sphere between P and another point P i , while r 0 and γ are parameters. This gives the probability of finding a mass at point P due to other masses located at points P i . This algorithm is repeated an arbitrary number of times. In our case, we will only generate clusters in configurations containing 600 sources and repeat the algorithm 6000 times. The parameters we used are γ = 3.5 and r 0 = 0.05. We will henceforth refer to the DI model without any artificial clustering as uniform.

Dynamics for Configurations with Clusters
Since clustering is only meaningful for many sources, we have exclusively generated clusters in configurations with 600 black holes. The LRS circumference dynamics are plotted in figure 20 in purple for 50 random configurations. A general tendency one can notice is that the purple curves are located further away from the homogenous limit given by the FLRW model, than the orange curves. The median initial length for the uniform case is 0.96(L F LRW ) 0 with a 60% confidence interval ranging from 0.92(L F LRW ) 0 to 0.99(L F LRW ) 0 . This is closer to 1 than the clustered case, with a median of 0.82(L F LRW ) 0 and a confidence interval ranging from 0.76(L F LRW ) 0 to 0.87(L F LRW ) 0 . In addition, the clustered models have in general a longer "lifespan". The median time for the uniform model is 0.19M p , but 74% of the clustered models encounter a singularity later than this. A simple explanation is that clustered groups of black holes are treated as an effective discrete mass source. This causes a "grainier" mass distribution and moves the model away from the homogenous limit.  [4]. Figure 21 plots the Hubble and deceleration parameters. Turning our attention to the Hubble parameters, 60% of the clustered configurations in purple lie above the FLRW curve at t = 0.2M p , compared to 98% of the uniform configurations in orange. Comparing with figure 14, decreasing the number of sources is accompanied by a decrease in the percentage of curves lying above the FLRW value at t = 0.2M p , from the aforementioned 60% down to 12% for 8 sources. Thus such behaviour is more closely associated with configurations containing fewer sources than 600. This further strengthens our interpretation that a group of clustered black holes are considered as a single effective discrete mass.

Concluding Remarks
In earlier work on cosmological models with discrete sources, the focus has been on source distributions adapted to regular tessellations of spatial hypersurfaces. In such models, the sources reflect the regular nature of the models by being positioned at lattice centers. The LRS curves present in those models allow for a simplification of their dynamics due to the LRS symmetry. In this paper we have presented the DI cosmological model consisting of discrete sources in the form of Schwarzschild black holes having an irregular distribution. To maintain the possibility to do simplified dynamics, we have used a special tessellation of S 3 in which there is a single LRS curve. This work therefore represents a generalization of previously discussed discrete models by allowing sources to be distributed irregularly while retaining the possibility of analytic analysis of the dynamical evolution.
Our simulations suggest that initial conditions, such as the relation between total mass and size of the universe, approach the FLRW limit as the number of masses increases (see figure 4). However, the approach is slower than that of the regular models. Also, the approach in our models is from below rather than from above in the regular case. This is possibly related to the different types of tessellations used. When it comes to the dynamics, the DI models, as expected, exhibit a richer variety of behaviour than that seen in models with a regular configuration of mass sources. In fact, we find that the system is quite sensitive to the initial mass configuration. Other interesting behaviour is the appearance of negative deceleration parameters.
The DI model has proven to be a simple and flexible model which permits the investigation of a large number of discrete sources distributed in various configurations. By using DI models with a larger number of sources one could determine more precisely the relationship between initial mass configuration and the resulting dynamics. It may also be possible to improve the DI model by using  Figure 21: The above and bottom plots depict respectively the evolution of the Hubble parameter and deceleration parameter for four different models with 600 sources. The purple curves correspond to the DI model with clustered configurations while the orange curves correspond to the uniform DI model. The black dashed curve and the green curve correspond respectively to the matter dominated spherical FLRW model and the regular model with 600 sources found in [4]. geodesic slicing in order to avoid coordinate singularities. One question which is raised by the DI model is whether the expanding acceleration sometimes observed is shared by other regions beside the LRS curve. Occurrence of large scale acceleration could be due to structure formation, often referred to as cosmic backreaction, as discussed by several authors (see for example [17]). Our results for the DI models with built-in clustering show that their dynamics generally deviates more from that of the FLRW model compared to models without clustering even though the particular results for acceleration are less clear. Illuminating these issues could improve our understanding of how mass inhomogeneities affect the large scale geometry and dynamics of the universe.
One obvious drawback of our DI model is that it has a certain amount of intrinsic anisotropy. This is due to the presence of a preferred great circle as part of the construction leading to a preferred global direction. This drawback can actully be overcome by noting that besides the DI model used here it is also possible to construct DI models based on a regular lattice by the same method. To preserve the LRS property in that case when performing reflections with irregular configurations, it is necessary that the number of cells meeting at an edge is an even number. The only regular model which fulfills this requirement is the 16-cell. However, all the regular models could actually be used for this purpose if one relaxes the regularity slightly. This can be done by noting that each 3-cell in a regular tessellation is intersected by a number of symmetry planes. These planes divide the 3-cell in subcells called chambers [18]. Any two chambers in a given regular tessellation are either identical or mirror images of each other. Since adjacent chambers are related by a reflection, one can pick one chamber as the seed cell and then distribute its configuration by reflection to all the other chambers. By such a construction based on the chambers instead of the regular cells themselves, the total number of sources corresponding to a given number of sources in a seed cell can be increased significantly. For example, a regular tetrahedron has 24 chambers and so a 600-cell has 24 × 600 = 14400 chambers.