The magnetic part of the Weyl tensor, and the expansion of discrete universes

We examine the effect that the magnetic part of the Weyl tensor has on the large-scale expansion of space. This is done within the context of a class of cosmological models that contain regularly arranged discrete masses, rather than a continuous perfect fluid. The natural set of geodesic curves that one should use to consider the cosmological expansion of these models requires the existence of a non-zero magnetic part of the Weyl tensor. We include this object in the evolution equations of these models by performing a Taylor series expansion about a hypersurface where it initially vanishes. At the same cosmological time, measured as a fraction of the age of the universe, we find that the influence of the magnetic part of the Weyl tensor increases as the number of masses in the universe is increased. We also find that the influence of the magnetic part of the Weyl tensor increases with time, relative to the leading-order electric part, so that its contribution to the scale of the universe can reach values of ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}1%, before the Taylor series approximation starts to break down.


Introduction
equations that govern Newtonian cosmologies [1], but also has strong support from rigourous constructions that are based on perturbative expansions of Einstein's field equations [2,3]. Nevertheless, there are good reasons to be interested in relativistic effects in cosmology. These include the apparent existence of Dark Energy, as well as the dawn of the age of precision cosmology. Put bluntly, relativistic effects need to be understood in order to have faith in the cosmological models that we use to interpret observational data.
Geometrically, the electric part of the Weyl tensor is sufficient to determine the Newtonian part of the free gravitational field. The magnetic part of the Weyl tensor, on the other hand, describes other aspects of the relativistic gravitational field [4] (both of these objects are defined in Sect. 3, below). Well known relativistic effects such as frame-dragging and gravitational radiation require a non-vanishing magnetic Weyl tensor in order to exist, but the effects of the magnetic part of the Weyl tensor on the large-scale expansion of the Universe are still largely unknown. This is partly due to an absence of realistic cosmological models that have non-zero magnetic Weyl curvature, and that could therefore be considered to be non-silent [5]. This situation is, however, changing.
It has recently been demonstrated that cosmological models that contain regularlyarranged discrete masses can generate a non-zero magnetic Weyl tensor, even if no such curvature existed in their initial data [6]. The authors of this study considered the effect that this tensor has on the large-scale expansion of a universe that contains eight black holes, using both the leading-order term in a Taylor series expansion, and by numerically integrating the evolution equations of the initially silent space. We extend their study by calculating both the leading-order and next-to-leading-order parts of the relevant series expansion, and by using these results to determine the effect of the magnetic part of the Weyl tensor on the large-scale expansion of all cosmological models that contain regularly arranged discrete masses in a closed cosmology. The phenomenology of such models is interesting, as they can be considered a first approximation to the type of universe within which we actually live. More mathematically, they provide a tractable way to formulate n-body cosmology as a relativistic initial value problem. They therefore provide an ideal arena within which to study relativistic gravitational effects in cosmology.
Of course, the magnetic part of the Weyl tensor is a frame-dependent object, and its magnitude will change when different sets of observers are considered. However, in cosmology it is natural to consider sets of time-like curves that are both geodesic, and (in some sense) comoving with the objects that exist within the space-time. It is with respect to observers following such a set of curves that one can talk about "cosmological expansion", and it is with respect to the same set of curves that we will talk about "the magnetic part of the Weyl tensor". The influence of the magnetic part of the Weyl tensor on the cosmological expansion is of special interest because virtually all cosmological solutions of Einstein's equations are derived under assumptions that force it to be zero (for cosmologically interesting congruences of curves). Such situations are, however, highly unrealistic, as observers that follow the geodesic set of time-like curves that describe the world-lines of real galaxies will certainly experience the consequences of this part of the curvature of space-time, which is in general non-zero. We must therefore precisely quantify the effects of the magnetic part of the Weyl tensor on the large-scale expansion of space, if we are to fully understand the recessional velocity of galaxies in the real Universe (and all of its associated consequences).
In Sect. 2 we describe in detail the cosmological models that we will be using in this study. They consist of regularly-arranged discrete masses in a closed space, and are formulated as a relativistic initial value problem. In Sect. 3 we then introduce the compact formalism that will be used in the rest of the paper. This starts by using a 1 + 3 decomposition of the geometric and kinematic quantities involved, and finishes with a 1 + 1 + 2 decomposition that can be efficiently used to exploit the symmetries of the problem. The explicit form of the symmetry restrictions are then determined in Sect. 4, where it is shown that the influence of the magnetic part of the Weyl tensor on the relevant evolution equations does not necessarily vanish. Section 5 then contains some lengthy calculations to determine the coefficients of a Taylor series expansion that can be used to include the effect of the magnetic Weyl tensor on the large-scale expansion of space. In Sect. 6 we present the numerical results for each of the lattice universes that we will be considering, before concluding in Sect. 7.
Throughout the paper we use the first half of the Latin alphabet to denote spacetime coordinate indices, and the second half to denote spatial coordinate indices. Greek letters are reserved to denote spatial frame indices, where they are required.

Lattice models of the universe
The space-times we wish to consider are those that have come to be known as "lattice models", in some parts of the cosmology literature. These models have a periodic structure, and are constructed from a number of regularly arranged cells, each of which are identical to one another. Such models have been studied in a number of different contexts over the past few years, as they offer a simple enough setting to provide concrete answers to questions of direct physical interest.
The particular type of lattices we wish to consider are those in which the topology of the cosmological region is a hypersphere. Furthermore, we want to consider situations where the only mass present is in the form of black holes, without angular momentum or electric charge. We position one single black hole at the centre of every lattice cell, so that the cosmological region as a whole contains a regular array of these objects. The geometry of the space-time can then be treated as a vacuum solution of Einstein's equations, with its scale and expansion being prescribed accordingly.
In the rest of this section we will describe the specifics of the six different lattice configurations that we will be studying, as well as the solutions to the constraint equations, and the existence of curves within the space-time that display special rotational symmetries. This recaps results from previous work [7][8][9].

The six closed lattices
Before solving for the geometry of space-time, we need to understand how to divide a hyperspherical universe into n identical cells. If we choose each cell to be one of the five possible Platonic solids (i.e. one of the five possible regular convex polyhedra), then there are six different ways that we can arrange these cells to form a closed lattice. Each of these configurations corresponds directly to one of the six regular, convex polychora that exist in four dimensions. They are [10]: This lattice consists of five tetrahedral cells. Three cells meet along every  cell edge, and there are a total of five vertices (the points where cell corners meet). This is the smallest lattice that can be constructed from regular convex polyhedra. -8-cell. This structure is composed of eight cubes, arranged so that three cells meet along every cell edge, giving a total of sixteen vertices. It corresponds to the space-time studied numerically in [6,11]. -16-cell. This lattice is made from sixteen tetrahedral cells, with four cells meeting at every cell edge. It is dual to the 8-cell (meaning that interchanging vertices and cell centres gives an 8-cell lattice). -24-cell. The 24-cell is constructed from twenty four octahedral cells. Three of these cells meet along every cell edge, and there are a total of twenty four vertices. The 24-cell is self-dual, as it has the same number of cells and vertices. -120-cell. This structure contains one hundred and twenty dodecahedral cells, with three cells meeting at every edge. There are six hundred vertices in this structure, meaning that it is dual to the 600-cell, described below. -600-cell. This is the largest of all lattices that fits our specifications. It is built from six hundred tetrahedral cells, with five of these meeting along every edge.
In the rest of this section we will outline how Einstein's constraint equations can be solved for these lattices, as well as how curves that exhibit local rotational symmetry can be identified.

Initial data
If we treat general relativity as an initial value problem, then we can split Einstein's equations into constraint and evolution equations. The former of these can then be written as a Hamiltonian constraint and a momentum contraint, such that, in vacuum, where (3) R is the Ricci scalar of the geometry on the initial hypersurface, Σ, and K i j is its extrinsic curvature in the embedding space-time. The D i derivative is used here to denote a covariant derivative within Σ, and K is the trace of K i j . Equations (1) and (2) are, in general, very difficult to solve. However, we can simplify the problem dramatically if we choose our initial hypersurface carefully. For example, if we choose it to be symmetric under time reversal then we automatically have K i j = 0. The Hamiltonian constraint then reduces to while the momentum constraint is identically satisfied. This equation is much easier to solve, and is the basis of the study of geometrostatics (the gravitational analogue of the field of electrostatics [12]). We will assume initial data of exactly this type in what follows, and find solutions to the constraint equation (3) that describe the maximum of expansion of our lattices. These results were presented for the first time in [7]. A metric ansatz for the initial geometry can be chosen as follows: where dΩ 2 = dθ 2 + sin 2 θ dφ 2 is the line-element of a 2-sphere, and where ψ = ψ(χ, θ, φ) acts as the square root of the scale factor. The time-symmetric Hamiltonian constraint then reduces to∇ wherem is a constant. This function has a pole at χ = 0, and looks very much like the gravitational field one might expect for a point-like particle located at that position. Of course, there is nothing special about the point χ = 0, and we can generate any number of similar functions by rotating the coordinate system so that {χ, θ, φ} → {χ i , θ i , φ i }, and then replacing χ by χ i in the denominator. This gives the solution for N such particles as where χ i should be understood as the value of χ in a coordinate system that has been rotated so that the ith particle appears at χ i = 0. This is just what we need to construct our lattice cosmologies, as it is a solution that corresponds to N arbitrarily located black holes, in a universe that is at its maximum of expansion.
In the sections that follow, we will use the solution given by Eqs. (4) and (6) as initial data for our cosmology. This initial data is exact, but has some subtleties to its interpetation. The parameterm i that appears in Eq. (6) looks like it corresponds to the mass of a particle positioned at χ i = 0, but it is not the proper mass of the ith black hole. This must be found by looking at the geometry of the space in the limit χ i → 0, and comparing this to a time-symmetric slice through the Schwarzschild solution. One then finds that the proper mass of each of the black holes is a function of the location andm i of each of the other black holes. The parameterm i can then be interpreted as an "effective mass", which includes information about the gravitational field of the gravitational potential energy, as well as that from the proper mass of each of the individual black holes [7,13]. As we want all of the black holes in our lattice models to have the same mass, we choose m i = m for all values of i.

Curves with a rotationally symmetric tangent space
There exists a number of curves, within our lattices, that display local rotational symmetry (LRS) [8,9]. Perhaps the simplest of these curves to visualize is the example of a cell edge. In each of the lattice structures listed above there were between three and five identical cells meeting along every cell edge. This means that the geometry of space-time in the vicinity of these curves must be invariant under a rotation that interchanges the cells. If we pick a set of orthonormal frame vectors for our 3-space {e 1 , e 2 , e 3 }, and align e 1 with this special curve, then the rotational symmetry implies [8] where V α and T αβ are the spatial frame components of any vector V, and any symmetric tensor T. Furthermore, the commutators of the basis vectors are restricted by [8] γ 1 21 = γ 1 31 = 0, γ 2 12 = γ 3 13 and γ 3 12 = γ 2 13 , where [e α , e β ] = γ αβ e . The results in Eqs. (7) and (8) are only valid in the tangent spaces of the points that make up these LRS curves, but they are significant restrictions. As well as cell edges, there are other curves in our lattices that also exhibit rotational symmetry [9]. They are given by the lines that connect the horizon of the black hole at the cell centre with the cell corner, and the line that connects this horizon with the centre of a cell face. Such curves are illustrated in Fig. 1, for the particular case of a cubic cell. They will be the curves that we solve for in Sect. 6. The curves that connect the centre of the cell with the cell corner, and the centre of a cell face, are formally infinitely long in the initial data, which is why we have chosen to cut them off at the black hole horizon. They then have finite proper length. Also displayed is a curve that connects the centre of a cell edge with a cell corner (red). These are the three curves whose evolution will be solved for in Sect. 6 (colour figure online) The location of the horizon will be determined by assuming the horizon to be initially non-expanding. In this case, the location of the horizon along LRS curves is found by looking for solutions to the following Eq. [8]: where e 1 is again aligned with the direction of the LRS curves, and E 11 is the frame component of the electric part of the Weyl tensor. The location of the horizon at subsequent times is found by propagating null geodesics along the LRS curves, from the initial location (as can be seen to be required by the null Raychaudhuri equation [8]).

Covariant formalism and field equations
In order to exploit the symmetries of a space-time it is often convenient to decompose the geometric objects that exist within it. This can be achieved by picking out vector fields that are either aligned with, or orthogonal to, certain invariantly defined subspaces, and then by defining projection operators associated with these vectors. Within the field of cosmology, the most widely applied decomposition of this type is the 1+3decomposition pioneered by Ehlers [14], and developed by Ellis and others (see e.g. [15,16]). If there also exists a preferred space-like direction within the space, as in LRS space-times [17], then the 1 + 1 + 2 decomposition of Clarkson becomes a useful tool [18]. Going further, and decomposing using a full set of four mutually orthogonal unit vectors, one is then led to the orthonormal frame approach [19]. In this section we will outline the 1 + 3 and 1 + 1 + 2 approaches, as relevant for our study.

The 1 + 3 decomposition
If there exists a time-like unit vector field, u, then we can define tensors that project parallel and othogonal to this vector, respectively, as where δ a b indicates the Kronecker delta. We can also project all derivatives along and orthogonal to u a , respectively, by defininġ where S a...b c...d is any tensor. It can immediately be seen thatḣ ab = 2u (a u b) and D a h bc = 0, which can be considered the evolution and constraint equations for h ab .
The covariant derivative of u a can then be decomposed into its irreducible components, such that whereu a = u b ∇ b u a is the acceleration vector, Θ := u a ;a is the expansion scalar, σ ab :=u (a u b) + ∇ (a u b) − 1 3 Θh ab is the symmetric trace-free shear tensor, and ω ab := u [aub] + ∇ [a u b] is the skew-symmetric vorticity tensor. These last three quantities describe the rate at which the integral curves of u a expand, squash and twist around each other. It is straightforward to verify that U ab σ bc = U ab ω bc = 0 = U abu b , meaning that all of these quantities are already projected orthogonal to u a .
The fundamental object of interest for describing the free gravitational field, and which is required to complete a description of the vacuum scenario, is the Weyl tensor, C abcd . The existence of the time-like vector field u can be used to split the Weyl tensor into "electric" and "magnetic" parts, defined respectively by E ab := C acbd u c u d and H ab : where abc = η abcd u d is the spatial skew symmetric permutation tensor (defined such that η 0123 = − √ |g|). The electric part of the Weyl tensor contains all information about the tidal forces due to gravity. The magnetic part, on the other hand, contains all other information about the Weyl curvature. We also have˙ abc = η abcdu d and D a abc = 0, which can be thought of as the evolution and constraint equations for abc .

The vacuum field equations
In order to write the field equations in the most concise form possible, it is convenient to introduce the definitions of the covariant divergence and curl, which when operating on rank-2 tensors are given by It is also convenient to introduce the projected, symmetric, trace-free part of a rank-2 tensor by using angular brackets, such that Note that this means σ ab = σ ab , E ab = E ab , and H ab = H ab , as these tensors are already projected, symmetric and tracefree.
With the definitions of all of these quantities and operators in hand, it is now possible to write the vacuum evolution equations for a geodesic (u a = 0), irrotational (ω ab = 0) flow as [20] with the constraint equations taking the form The bracket operators in these equations are defined as [S, T ] a := abc S b d T cd , and S, T ab := S c a T b c . This extremely compact way of writing the field equations will greatly aid the manipulations we perform in Sect. 5. It also has the distinct advantage that all quantities involved have a direct physical interpretation, allowing for a better intuitive understanding of the physical set-up.

The 1 + 1 + 2 decomposition
Just as in the previous sections, and following Clarkson [18], if there exists a space-like unit vector field, n, then we can decompose all geometric objects into parts that are parallel and orthogonal to this field. The goal here is to first perform the decomposition with respect to u a , and then to take the tensors projected into the space orthogonal to u a and decompose them further with respect to n a . To this end, we introduce the following projection tensors: where h ab is the projection tensor from above. The tensor N a b can operate on any geometric object, and projects it along n a . The tensor f a b , on the other hand, only makes sense as a projection tensor when it operates on objects that have already been projected using h a b . If n a is surface forming, then f a b projects such objects onto the 2-dimensional "sheets" orthogonal to n a .
Under these definitions, we can irreducibly decompose any spatially projected vector, V a = h a b V b , into a scalar part V := V a n a , and a transverse vector part The projected vector can then be written as shows that V a has no parts in the sheet, and ⊥ V a has no part parallel to n a . The other geometric objects we wish to decompose with respect to n a are the symmetric trace-free tensors defined above. Any such tensor can be irreducibly decomposed into a scalar • S := N ab S ab , a transverse vector † S a := f a b n c S bc , and a transverse and tracefree tensor ‡ S ab : The full symmetric and trace-free tensor can then be written in terms of these fundamental objects as These three objects are each symmetric, and trace-free with respect to h ab . They are the only three such objects that can be constructed from a single power of • S, † S a and ‡ S ab , and as such they specify a unique decomposition of S ab .
In the next section we will consider the restrictions that reflection symmetry imposes upon these objects, and why this means that H ab can have an effect on the large-scale expansion of space.

Reflection symmetric surfaces
The presence of reflection symmetric surfaces can greatly simplify the formulation of problems in general relativity. Indeed, it is the presence of a reflection symmetry in the time-like direction that allowed the constraint equations to be solved for a lattice universe [7]. Reflection symmetry in the space-like directions normal to certain 2 +1-dimensional surfaces within the space-time also lead to the notion of "piecewise silence", as gravitational waves are forbidden from travelling between lattice cells, or the various chambers within cells [9]. In this section we will investigate the properties of reflection symmetric surfaces further, by formulating their consequences for geometric objects in the 1 + 1 + 2 decomposition, and by considering the situation where multiple reflection symmetric surfaces intersect (as often occurs in our lattice universe construction).

Polar tensors and axial tensors
When determining the consequences of reflection symmetry for tensorial quantities it is important to distinguish between polar and axial tensors. If a tensorial object can be defined without reference to a specific frame, then that object is said to be a "polar tensor". By contrast, if such an object is defined with respect to some fixed direction or combination of directions, then that object is said to be an "axial tensor". The 3-dimensional Levi-Civita tensor, abc , is the archetypical axial tensor. If {e 1 , e 2 , e 3 } constitute a set of orthonormal spatial frame vectors, then the frame components of this tensor are given by 123 = 1, and all permutations. In other words, the Levi-Civita tensor defines a preferred orientation of the spatial frame. A reflection, or any other improper orthonormal transformation, changes the orientation of the spatial frame. Reflections therefore change the sign of abc , meaning that it should be properly identified as a pseudotensor, rather than as a tensor. This distinction is extremely important when determining the consequences of reflection symmetry on geometric quantities.
To illustrate this in a formal way, we note that an inversion can be represented by the orthonormal transformation matrix I = −1, where 1 denotes the 3-dimensional unit matrix. Applying the standard tensor transformation formula for an orthogonal transformation matrix O a b would then give If we take S abc... = abc and O = I then this law would give us O( abc ) = − abc . The sign change indicates that the transformed object, O( abc ), has opposite parity to the untransformed object, abc . This is an unwanted feature if we require the parity of the Levi-Civita tensor (right or left-handedness) to be invariant under any orthogonal transformation. The standard way around this is to modify the transformation rule for pseudotensors, so that it becomes [21] O(S abc . . We can then transform tensors according to Eq. (28), whilst transforming pseudotensors according to the expression in Eq. (29). The extra factor of det(O) in this latter equation means that any quantity that contains an odd number of Levi-Civita tensors picks up an extra minus sign under any improper orthogonal transformation, such as a reflection.

Implementing reflection symmetries
Consider a reflection symmetry in a spatial hypersurface. Such a symmetry has fixed points that form an invariant symmetry surface, S . Let {n, k, l} be an orthonormal frame adapted to the symmetry in such a way that n is orthogonal to S , and consequently that k and l are parallel to S . The reflection symmetry in the tangent space of a point in S can then be represented by the orthonormal frame transformation matrix R a b = diag(−1, 1, 1) for a polar tensor, and byR = −R for an axial tensor. The restrictions on a tensor Q, imposed by invariance under a reflection symmetry, is then determined by the equations depending on whether Q is a polar or an axial tensor, respectively. We note that, for the specific fields that appear in Eqs. (17)-(24), the set {Θ, σ ab , E ab } are polar, while H ab is axial (as it involves the Levi-Civita pseudotensor). The only operation in the field equations that changes the parity of an object is the curl, which means that curlH ab is polar, while curlσ ab and curlE ab are axial tensors. Decomposed tensor parts, such as † E ab , always have the same parity as the full tensor from which they are derived.
To calculate the action of the reflection R, let us first consider the case of scalars. We note that the only change possible for a scalar is the change of sign that affects pseudoscalars under an improper transformation, according to the pseudoscalar version of Eq. (29). This implies that all pseudoscalars must vanish on S , leading to the restrictions • (curlσ ) = • (curlE) = • H = 0. To see the effect on vectorial and tensorial objects we can express the transformation matrix in the form The action on the vector part of a tensor, S ab , is then given by 1 Therefore, in order to fulfil the symmetry condition given in Eq. (30), the vectorial part, † S a , of a polar tensor, S a , must vanish on S . It then follows that † σ a = † E a = † (curlH ) a = 0 on S . Applying the reflection transformation in an analogous way to the tensor part of S ab gives R( ‡ S ab ) = ‡ S ab . Hence, the polar tensor part is not affected by the reflection, and is therefore unrestricted. The axial tensor part, on the other hand, is required to vanish, due to the determinant appearing in Eq. (29). This implies that ‡ (curlE) ab = ‡ H ab = 0. In summary, the tensor parts remaining on a reflection-symmetric surface, after the symmetry restrictions imposed by Eq.

Intersection of two symmetry surfaces
We now consider the symmetry restrictions imposed at points where two reflection invariant surfaces meet. Let us call these two surfaces S k and S l , where the subscripts k and l indicate their normals. The intersection of these two surfaces forms a curve, C , and reflection symmetry about each surface ensures that they meet at right angles to each other. In this case we can choose to adapt the frame in such a way that n is parallel to C . The quantities • (curlσ ), • (curlE) and • H now have a different interpretation, but they are again pseudoscalars. We therefore have that To analyze the restrictions on vectorial and tensorial objects, we note that the reflection transformations associated with k and l can be expressed in the forms Applying R k to the vectorial part of the 1 + 1 + 2 decomposition gives Applying the symmetry condition from Eq. (30), for a polar tensor, then leads to k a † S a = 0 on S k . A similar argument for the reflection symmetry R l gives l a † S a = 0 on S l . It follows that † S a = 0 on C , as both of its components must vanish. This implies the restrictions † σ a = † E a = † (curlH ) a = 0.
For the vector part of an axial vector we get In this case, the symmetry condition from Eq. (30) takes the traceless symmetric form Only the vectorial component of C ab is non-zero, which means the symmetry condition in Eq. (35) then reduces to † C a = 1 2 This tells us that the component of † S a in the l-direction vanishes on S k . Similarly, the symmetry R l implies that the component of † S a in the k-direction vanishes on S l , and consequently that the vectorial parts of all axial tensors vanish on S . This gives Let us now consider the consequences of the symmetry conditions, from Eq. (30), on the tensorial part of a polar tensor. We find that the reflection operation gives where we have used the relation ‡ S ab The symmetry condition in Eq. (30) then implies ‡ S ab k a l b = 0. The tensorial parts ‡ σ ab , ‡ E ab and ‡ (curlH ) ab are therefore diagonal, with ‡ S ab k a k b as their only independent component. As the right-hand side of Eq. (37) is symmetric in k and l, there are no further restrictions implied by the symmetry action of R l . Finally, we consider the action of the reflection symmetries on axial tensors. They giveR so that the symmetry conditions in Eq. (30) then become ‡ S ab = 2( ‡ S cd k c l d )k (a l b) .
This is the condition that ‡ S ab k a l b is the only nonvanishing independent component of ‡ S ab . This restriction applies to ‡ (curlσ ) ab , ‡ H ab and ‡ (curlE) ab .
In summary, at the intersection of two reflection symmetric surfaces we have that the only non-vanishing tensor parts are given by Θ, • σ , • E and • (curlH ), the single independent component of the diagonals of ‡ σ ab , ‡ E ab and ‡ (curlH ) ab , and the single independent off-diagonal components of ‡ (curlσ ) ab , ‡ (curlE) ab and ‡ H ab 2 .

Intersection of three or more symmetry surfaces
When three or more reflection symmetric surfaces meet the restrictions on geometric objects are much stricter. Such situations arise along LRS curves when the symmetry along the intersecting surfaces is threefold or higher, as was studied in some detail in [8]. Here we adapt our frame to the symmetry curve in the same manner as in the twofold symmetry case, considered above. As we know that all the vectorial and tensorial parts of every geometric object must vanish in this case [8], we only need to consider the scalars. As before, the axial scalars must once again vanish, implying that • (curlσ ) = • (curlE) = • H = 0 . It then follows that curlσ ab = curlE ab = H ab = 0. The only remaining parts are therefore given by • σ , • E and • (curlH ).
If we now re-cast Eqs. (17)- (19), in terms of the quantities that remain after symmetry restrictions, then we finḋ These equations can be further simplified by defining the expansion scalars parallel and perpendicular to n a : H := 1 3 Θ + • σ and H ⊥ : • σ , respectively. The evolution equations for these two quantities are then given bẏ while the evolution equation for the source term, • E, is given by The scalar • (curlH ) therefore acts as a source for • E, which is itself the source of the expansion along S . This completes our study of the restrictions imposed at reflectionsymmetric surfaces. In the following sections we will study the evolution of the nonvanishing parts of the covariant scalars, and in particular the effect that • (curlH ) has on the expansion of LRS curves with threefold, or higher, reflective symmetries.

Taylor series expansion of curl(H)
We want to solve the evolution equations (43)-(45), with the initial conditions specified in Sect. 2. In terms of the variables that result from the 1 + 3 decomposition, these initial conditions can be specified as together with the following conditions on E ab : The first set of initial conditions in Eq. (46) all follow immediately from the requirement of time-reversal symmetry in the initial data [8]. The second set, involving time derivatives, follows from Eq. (17), and from the fact that the Cotton-York tensor on the initial hypersurface is given by * C ab = −Ḣ ab [19]. The conditions in Eq. (47) then follow from the evolution equation (20), and from the constraint equation (23). If one were to take • (curlH ) = 0 along LRS curves, as was done in [8], then Eqs. (43)-(45) form a closed system, and the initial data specified in Sect. 2 is sufficient to determine H , H ⊥ and • E at all future times (up to singularities). If • (curlH ) = 0, however, then this is no longer true. The value of • (curlH ) must be determined at future times, in order to evolve the system of evolution equations (43)-(45). This can be done by either numerically solving the full set of Einstein equations, or by performing a Taylor series expansion of • (curlH ) at the initial hypersurface [6]. In this paper we opt for the latter of these two methods, and in the remainder of this section we determine the coefficients of such an expansion using the compact notation from Sect. 3.
In order to perform such an expansion, we first note that the time derivative of any spatial unit vector can be written aṡ where Ω α := 1 2 αβγ e β aė γ a is the angular velocity of the spatial unit vectors, with respect to a set that are Fermi propagated along u a . If we choose our flow so that it is geodesic (u a = 0), and take our unit vector n a to belong to a set of basis vectors that are Fermi-propagated (Ω α = 0), then we haveṅ a = 0 at all points along the evolution of all of our LRS curves.
Letting (curlH ab ) (n) stand for the nth covariant derivative of curlH ab along u a , this means, in particular, that for any number of time derivatives. We can therefore Taylor expand curlH ab , and simply contract the result with n twice, in order to find the series expansion of • (curlH ).

First order
The Taylor series expansion we wish to investigate is given by where (curlH ab ) (n) should be understood to be evaluated on the initial hypersurface. The time parameter, t, is the proper time that would be measured by an observer following the integral curves of u a . The n = 0 term in Eq. (50) is zero, as H ab must vanish at all points in the initial hypersurface. In fact, we can immediately see that every term that corresponds to an even value of n must vanish from this equation. This follows directly from the requirement of time-reversal symmetry. It now remains to evaluate the odd values of n. At first order we get where we used the commutation relation between the curl operator and time derivatives (Eq. (A18) from [20]): where S ab is any projected symmetric tracefree tensor, S ab = S ab . The initial conditions in Eqs. (46) and (47) can then be seen to imply at all points on the initial hypersurface. This means that the leading-order term in Eq. (50) cannot occur at any lower value than n = 3.

Third order
In order to find the third-order term in the Taylor series expansion, from Eq. (50), we must differentiate each term in Eq. (51) twice. For the first term we get where we have again used the commutation rule from Eq. (52), as well as the following identity (Eq. (A16) from [20]): where f is any scalar function. Applying the Leibniz rule, and plugging in the initial conditions, we then get −(curl(Θ H ab )) | t=0 = 0. For the second term in Eq. (51) we get where we twice applied the identity from Eq. (52), and where we have also used vertical bars around indices to exclude them from symmetrization operations. The initial conditions stated in Eq. (46) then imply that (curl(σ c a H b c )) | t=0 = 0, as all terms contain quantities that vanish at t = 0. The second time derivative of the fourth, fifth and sixth terms in Eq. (51) can also be seen to vanish, using similar logic.
It is therefore only the third term on the right-hand side of Eq. (51) that can be nonzero, which means (curl(H ab )) | t=0 = −(curl curl(E ab )) . Commuting the curl operator with the time derivatives, using Eq. (52) and the initial conditions from Eqs. (46) and (47), then gives This can be simplified further by using the following identity (Eq. (A4) from [20]): where (S 2 ) ab := S c a S bc . We then have, after using Eqs. (57)-(58) and the evolution equations (18)- (19), that −(curl curl(E ab )) | t=0 = 4curl curl(E 2 ) ab . This is the only term that is non-zero, after twice differentiating the right-hand side of Eq. (51).
The final result, for the third derivative of curlH ab on the initial hypersurface, can therefore be written as This is the first non-zero coefficient in the Taylor series from Eq. (50). Explicitly evaluating (curlH ab ) , in terms of functions that appear in the metric, is likely to result in a very long expression, for all but the most trivial geometries. It is therefore remarkable that the same quantity can be written down in such a simple way in Eq. (59). This is a result of the compact notation that we have used, which hides an enourmous amount of underlying complexity.

Fifth order
The third-order coefficient, found in the previous section, is expected to be the leadingorder term in the Taylor series expansion from Eq. (50), as it is the first one that does not vanish. In order to determine the regime in which the Taylor series provides a good approximation to the full geometry, however, it is necessary to find the next-toleading order term (when the leading-order and next-to-leading-order terms become the same size, then we expect the series expansion to break down). We therefore need to calculate the fifth-order term in the series expansion from Eq. (50). The fifth-order term will be given by performing four differentiations on each term in Eq. (51). We will do this by following the same procedure used in the previous section, where simplifications were obtained using the identities in Eqs. (52), (55) and (58), and the initial conditions in Eqs. (46) and (47). Each of the six terms in Eq. (51) will be considered separately, and then summed to give out final result.

First term
To find the fourth derivative of the first term in Eq. (51), it is useful to first determine the second derivative without applying the initial conditions. This is given by where we first expanded curl(Θ H ab ) using Eq. (55), and then applied the commutation relation from Eq. (52) twice. It can be seen from this equation that, if we differentiate twice more, then each term will contain at least one factor of a quantity that vanishes on the initial hypersurface. We therefore have so that the first term of Eq. (51) does not contribute to the fifth derivative of curlH ab .

Second term
Applying Eq. (52), it can be seen that the second derivative of the second term in Eq. (51) is given by This allows us to see, by inspection, that it is only the first term in this equation that will survive two further differentiations. After commuting these derivatives with the curl operator, and disregarding terms that vanish on the initial hypersurface, we then arrive at 3[curl(σ c a H b c )] | t=0 = 3curl ((σ c a H b  c ) ). The only part of (σ c a H b c ) that does not vanish in the initial data is given byσ c a (H b c ) . Using the evolution equations for σ ab and H ab , and the result of the previous section, in Eq. (59), then gives This is the simplest expression we can find for this term.

Third term
This is the most lengthy term to evaluate. We find, by repeated application of the commutation relation (52), that The complete evaluation of the third term requires two further derivatives, and we can see that the only terms that survive this operation are given by where we have made use of the initial conditions. In order to evaluate this term we note that, after some manipulation, we can write where we have again made repeated use of Eq. (52), as well as the following identity (Eq. (A12) from [20]): We have also made the definition (S 3 ) ab := S ac S cd S db . Substituing Eq. (66) into Eq. (65), and performing some further manipulations, then gives This is the lengthiest of the six terms being evaluated.

Fourth and fifth terms
The fourth and fifth terms, from Eq. (51), can both be seen to give This is because H ab and Θ, together with both their first and second derivatives, must vanish on the initial hypersurface. Two further time derivatives would therefore be required before either of these terms become non-zero at t = 0.

Sixth term
The final term evaluates to where we have used the commutation rules from Eqs. (52) and (55), as well as the evolution equations (18), (19) and (20). Putting all of these results together leads us to our final expression for the fifth order term in our Taylor series expansion of curlH ab : This equation is considerably less compact than the result for the third-order derivative, presented in Eq. (59), but is still remarkably short, given the complexity it conceals.
In what follows, we will use the expressions calculated in Eqs. (59) and (72) to evaluate the first two non-vanishing terms in the Taylor series from Eq. (50). This will allow us to not only determine the leading-order effect that H ab has on the expansion of LRS curves, but also to estimate when the series expansion approach breaks down. We will perform this analysis for all three classes of LRS curves, and in all six of the lattices that were described in Sect. 2.

The effect of curl(H) on cosmology
The effect that H ab has on the expansion of space can be estimated by considering its consequences on the scale factors, a and a ⊥ , defined implicitly by The function a is then the scale factor along our LRS curves, and a ⊥ is the scale factor in all perpendicular directions. Both are a function of proper time, along the integral curves of u, as well as position in space.
Following the method used in [6], and extending it to higher order, we can use Eqs. (43)-(45) to write the correction to the scale factors, from including the curlH ab term in Eq. (45), as and where all quantities on the right-hand side of these equations should be taken to be evaluated at t = 0. These expressions result from simultaneously Taylor expanding a , a ⊥ and • (curlH ) about the initial hypersurface, and substituting the results into Eqs. (43)-(45). The value of • E is then given by noting that on a time-symmetric hypersurface, and in vacuum, the electric part of the Weyl tensor is equal to the Ricci tensor of the 3-space, E ab = (3) R ab . The quantities involving derivatives of curlH ab are given, in terms of E ab and its derivatives, by Eqs. (59) and (72). We therefore have all the information required to evaluate all of the terms present in Eqs. (74) and (75). The LRS curves that we will consider in this section are those displayed in Fig.1. The red curve in this figure extends halfway along a cell edge, from the centre of an edge to a vertex. The green curve connects the horizon of the mass at the centre with the centre of a cell face, and the blue curve connects the horizon and a cell vertex. These curves are depicted in Fig. 1 for the special case of a cubic cell, which is the basic constituent of the 8-cell. However, similar curves also exist for the case of tetrahedral, octahedral and dodecahedral cells. We will not present specific diagrams of these cells, but instead rely on the reader to visualize the corresponding curves in each of these cases.
The proper length of each of our LRS curves, in a hypersurface of constant proper time, is then given by the integral where χ is the coordinate from Eq. (4), and where we have rotated the configuration so that the LRS curve under consideration is at constant θ and φ. In what follows, we will present the intial values of the functions • E, as well as the derivatives • (curlH ) and • (curlH ) , at every point along the three LRS curves depicted in Fig. 1. We will then calculate the proper length of each of these curves, as a function of proper time along the integral curves of u, until the Taylor expansion given in Eq. (50) breaks down. In calculating these results, we have substituted the expansion from Eq. (50) directly into the evolution equations (43)-(45). The indicator that we use to determine when the Taylor series approximation breaks down is when the first two non-vanishing terms on the right-hand side of Eq. (50) become equal. This happens at different times along each of the curves, depending on the spatial position that one considers. We choose the two endpoints of each of the curves, as depicted in Fig. 1, in order to follow the magnitude of the ratio of these two terms. This is a convenient choice as the endpoints of the curve are uniquely picked out by the geometry of the problem. We also expect these points to give a fair reflection of the validity of the Taylor expansion at intermediate points, as they are in the two most extreme environments available along each curve (in the case of the green and blue curves, one end touches a black hole while the other is equidistant from black holes). If the value of • (curlH ) or • (curlH ) vanishes at any one of these points, then we simply take the ratio of the relevant terms a very small distance away (1% of the distance along the curve, where all quantities are non-zero).

Curves along cell edges
The first set of curves we wish to consider are those that lie along the edges of our cells, as depicted by the red line in Fig. 1. In this study we will take the curve in question to begin at the cell vertex, and extend halfway along the length of the edge. We will then calculate the length of this curve, using the method described above, until the Taylor series approximation has broken down at both of its ends (i.e. at both a point near the cell vertex, and at the centre of a cell edge).
First let us present the relevant information about the electric and magnetic parts of the Weyl tensor along these curves, and on the initial hypersurface. The form of the electric part of the Weyl tensor, contracted twice with the space-like unit vector tangent to the edge, is displayed in Figs. 2 and 3. In these plots we have chosen the cell vertex to be located at χ = χ 1 , and the centre of the cell edge to be located at χ = χ 2 . This information is presented in two different plots, as there are two different functional forms for • E; those that have a non-zero derivative at χ = χ 1 , and those that have a vanishing derivative at χ = χ 1 . These two cases correspond to lattices in which the cell edges are non-contiguous and contiguous, respectively. The former of these two cases contains the 5-cell, the 8-cell and the 120-cell lattices, while the latter contains the 16-cell, the 24-cell and the 600-cell lattices. As well as the electric part of the Weyl tensor, we also need to know about the magnetic part of the Weyl tensor. This tensor vanishes on the initial hypersurface, but its derivatives do not. In Figs. 4, 5, 6, 7 we present information about the value of the  using Eq. (72). We again display these results for the lattices with contiguous and noncontiguous edges separately, as they have different functional forms. Interestingly, although the curves have similar shapes to those of • E, they do not all have the same sign (the curve for the 16-cell can be seen to be negative in Fig. 5, and positive in Fig. 7). One may also note that the magnitude of • (curlH ) and • (curlH ) are both much smaller than • E, in the cases where the lattices have contiguous edges (and in the chosen units).
Next, we use the information presented in Figs. 2, 3, 4, 5, 6 and 7 to evolve Eqs. (43)-(45). The lengths of the curves being considered in this section can then be calculated by finding the scale factors defined in Eq. (73), and integrating them as prescribed in Eq. (76). If the scale factor becomes zero at any point along the curve, we simply remove this point from the integral, as described in [8]. The result of all of this is shown in Fig. 8. The proper length of the curve is evolved in time, in this plot, and displayed using units of the total proper mass in the cosmology (i.e. the proper mass of each cell, multipled by the number of cells in the lattice). This choice of units allows us to display all six lattices together with just a single FLRW solution that has the same total mass as each of the lattices. The particular FLRW solution we choose to compare with our models is a dust-filled cosmology with positive spatial curvature.
Also indicated in Fig. 8 is the region where the Taylor series expansion of • (curlH ) breaks down. We indicate this on the plot by the solid line turning into a dashed line. We start the dashing when one end of the curve fails our convergence criterion (i.e. when the second non-zero term in Eq. (50) becomes as large as the first). We end the dashing when the other end of the curve also fails this condition. The dashed region therefore corresponds to a domain where some of the points on the curve obey the convergence criterion, but not all. In this region the Taylor series approximation is breaking down. It can be seen that the length of the solid and dashed region, for each line on the plot, is different from the others. Roughly speaking, the Taylor series expansion seems to break down more quickly as the number of cells in the lattice increases.
As well as considering the evolutions of the overall length of the cell edge, which are quite similar to those found in [8], it is also interesting to consider the magnitude of the effect of including • (curlH ) in Eq. (45). This effect is displayed graphically in Fig. 10, where we plot the difference between the curves shown in Fig. 8 and the curves that would have existed if we had neglected • (curlH ) (as was done in [8]). Once again, we used solid and dashed lines to indicate where the Taylor series expansion starts to break down. The effect of including • (curlH ) is most dramatic in the smallest lattices. In the 5-cell and 8-cell the difference can reach values of |δL| 0.01. In the larger lattices this difference is much less, although the Taylor series expansion breaks down much sooner. One may also note, however, that at the same cosmological time (measured in units of n × m) the value of δL is larger in the bigger lattices.
As the curves in Fig. 9 are difficult to distinguish, we have read off particular values for each of the lattices, and displayed these separately in a log-log-plot in Fig. 9. The three curves in this plot show the value of δL at t = m (black points), at t = 2m (red points), and at t = 3m (blue points), where m is the proper mass in any given cell. We have excluded points that correspond to locations on the lines in Figs. 8 and 9 that become dashed, so as to exclude regions where the series expansion of • (curlH ) breaks down. Although the black curve dips in the middle, the general trend seems to show that the effect of • (curlH ) becomes smaller as the number of cells in the lattice is increased (at least, when comparing the lengths of the curves at these values of t). The effect of • (curlH ) does, however, increase with time, reaching levels of ∼1%, in the smaller lattices at t = 3m.

Curves through face centres
Let us now set aside the curves at the edges of cells, and instead consider those that extend from the horizon of a black hole to the centre of a cell face. These are the Fig. 10 The difference between including • (curlH ) in Eq. (45) and neglecting it, as inferred from the length of a cell edge, and as a function of the number of cells in the lattice, n. This is shown at t = m (black points), t = 2m (red points), and t = 3m (blue points), where m is the proper mass contained within any one cell in the lattice. Points are excluded when the curves shown in Figs. 8 and 9 become dashed, rather than solid (colour figure online) Fig. 11 The value of • E at t = 0 along a curve that goes from the cell centre to the centre of a cell face, for the 5-cell (red), the 8-cell (orange), the 16-cell (yellow), the 24-cell (green), the 120-cell (blue), and the 600-cell (purple) (colour figure online) type of curves depicted by the green line in Fig. 1, for the particular case of a cubic cell. Such curves exhibit local rotational symmetry, in exactly the same way as a cell edge. Once again, we will evaluate • E and the derivatives of • (curlH ) along these curves. The former of these quantities is shown in Fig. 11, while the latter are shown in Figs. 12, 13, 14, 15. In these plots we have taken χ 1 to correspond to the centre of the cell (rather than the location of the horizon). The curves then extend through the horizon, which is initially located at the minimum of • E [8], and out to the edge of the cell, where they meet the centre of a cell face at χ = χ 2 . If the proper lengths of these curves were to be calculated from χ 1 to χ 2 , in the initial hypersurface, then we would find the result to be divergent. We therefore restrict our calculations of the proper length to the region of space exterior to the horizon (a finite quantity). To find the location of the horizon Fig. 16 The length of the curve that goes from the cell centre to the centre of a cell face, normalized by the maximum value of a corresponding curve in an FLRW universe that is filled with dust, and has positive spatial curvature. The six different lattice models displayed as before, black dotted line is the corresponding FLRW solution (colour figure online) at all times after t = 0 we simply propagate a null geodesic outwards from the initial location, along the LRS curve [8].
The form of • (curlH ) along each of these curves is shown in Fig. 12. The larger lattices are difficult to make out on this scale, so we also show them separately in Fig. 13. The form of these curves in more complicated than was the case for the cell edges. In the two smallest lattices it is apparent that there is a single maximum in the function, located at the centre of the cell face. For the larger lattices, the maximum is located somewhere between the horizon and the cell face centre, with the cell face centre itself becoming a local minimum. The behaviour of • (curlH ) is even more complicated, and is shown in Fig. 14 for all lattices, and in Fig. 15 for the 120-cell and 600-cell alone. It can be seen that both maxima and minima of this function exist, along the curves under consideration.
Using these data, we can again integrate Eqs. (43)-(45), and calculate the length of the curve (from the black hole horizon to the centre of the cell face) using Eq. (76). The results of doing this are shown in Fig. 16, for each of the six lattices. In this case the effect of the • (curlH ) term in Eq. (45) is again small, and the curves look very similar to their counterparts when this term is neglected [8]. Once more, we make the lines in Fig. 16 dashed when one end fails our convergence criterion, and we stop plotting the line once both ends fail. Somewhat surprisingly, the curves in this case extend slightly further than was the case with the cell edges. This means the Taylor series approximation is slightly better along these curves, and for the 5-cell it appears to be good enough to evolve the curve all the way until its proper length vanishes (when the horizon makes its way to the centre of the cell face).
In order to visualize the effect of the • (curlH ) term on the evolution of this curve, we plot the difference from its inclusion in Fig. 17. In this case, the difference in the proper length of the curve can be as much as ∼6% of its initial length (in the case of the 8-cell), and ∼3% for the 16-cell. This larger difference is, to some extent, a consequence of the Taylor series approximation lasting longer in this case, meaning that the cumulative effect of integrating Eq. (45) with an extra term is larger. Again, however, the value of δL for the larger lattices is too small to be seen in Fig. 17. We therefore read off its value at t = m, at t = 2m, and at t = 3m (where m is the proper mass of each of the black holes). This information is displayed graphically in Fig. 18.
As before, when presenting the values of δL in the log-log-plot shown in Fig. 18, we do not present points that would correspond to locations along the curve where the Taylor series approximtation has started to break down. There are, however, more points in the plot shown in Fig. 18 than there are in the one shown in Fig. 10. This is again due to the Taylor series approximation lasting longer in the present case than it did along the cell edges. On the other hand, the broad trends in the two plots are largely similar: The value of δL decreases as the number of cells in the lattice is increased, and increases as the time from the initial hypersurface is increased. The largest values of δL in this plot are again at the level of ∼1%, for the smallest lattices at t = 3m.

Curves through vertices
Finally, let us consider the remaining curves depicted in Fig. 1: these are the blue curves that connect the black hole at the centre of the cell to one of the vertices. In the absence of any better terminology, we will call these curves 'diagonals' (although this is probably only really a fitting description for the case of cubic cells). The location of the horizons will be determined in the same way as in the previous subsection. We find their initial position by looking for the minimum in the function • E, along our LRS curves, and then at subsequent times by propagating a null geodesic outwards from this location. In all of Figs. 19, 20, 21, 22, 23, 24 we take the cell centre to correspond to the point χ = χ 1 , and the cell vertex to correspond to the point χ = χ 2 . We remind the reader that these curves are formally infinitely long, so we restrict our integrations to the regions of space exterior to the black hole horizons.
The values of • E along the diagonals, for the six different lattices we are considering, are shown in Figs Just as with the cell edges, there are two different types of behaviour in this case, depending on whether or not the diagonal curve from one cell are contiguous with those of its neighbours, or not. It is still the case that the curves in the 5-cell, the 8-cell and the 120-cell are non-contiguous, while those of the 16-cell, the 24-cell and the 600-cell are contiguous. In the contiguous case one could image holding a mirror up to the right-hand side of the plots to determine the continuing form of the function in question, as the curve extends directly into the next cell. In the non-contiguous case the diagonal curve from one cell joins onto the end of the cell edge of its neighbours. One could therefore join the plots in  including • (curlH ) are shown in Fig. 26. The magnitude of the effect is again most pronounced for the smaller lattices, as in these cases it takes longer for the Taylor series approximation to break down. Again, it is hard to see the effect of • (curlH )   The length of a diagonal curve, normalized by the maximum value of a corresponding curve in an FLRW universe. The six different lattice models are displayed using the same colours as before, and the FLRW solution with the same total mass is displayed as a black dotted line (colour figure online) Fig. 26 The difference between including • (curlH ) in Eq. (45) and neglecting it. This is displayed by showing the effect it has on the length of a diagonal curve, again normalized by the maximum value of a corresponding curve in an FLRW universe. The six different lattice models are displayed as before

Discussion
We have considered the effect of H ab on the evolution of LRS curves in lattice models of the universe. These models treat the matter in the Universe as a collection of pointlike sources, and allow the formulation of cosmology as an initial value problem. They are therefore particularly well suited to the study of relativistic effects in cosmology, including the study of the evolution and effects of the magnetic part of the Weyl tensor.
We find that although the initial data of our models is silent (with H ab = 0), the evolution of the space is not silent. In particular, a curlH ab term appears in the evolu- tion equation for the electric part of the Weyl tensor, which in turn acts as the source for the evolution of LRS curves. This result was identified and studied numerically by Korzyński, Hinder and Bentivegna for the particular case of a universe that contains eight black holes [6]. We extend their study by calculating the leading-order and next-to-leading order terms in a Taylor series approximation that can be used to incorporate the effects of H ab on LRS curves, and by applying it to all possible regular arrangements of black holes in a closed cosmological model. The inclusion of the next-to-leading order term, in our study, allows us to estimate when the series expansion stops converging, and when numerical techniques need to be employed instead.
We find that the effect of H ab is small, while the series expansion remains valid, but grows with time. In particular, we find that while the effect of the H ab on the expansion of LRS curves can be at the level of 1% when the number of masses in the Universe is small (as in [6]), but that this number decreases as the number of masses in the universe is increased (when comparing at the same time, measured in terms of m). While the effect on the expansion rate is small for larger lattices, however, it does appear that the effect of H ab on the expansion is, in some sense, cumulative. That is, the difference in the curve length that results from including the effect of H ab tends to increase over time, as the magnitude of the magnetic Weyl tensor increases over time. When comparing lattices at the same cosmological time (i.e. at the same time, as measured in units of n × m), it therefore appears that the effect of curlH ab on the evolution of the LRS curve increases as the number of masses in the universe is increased. Whether or not similar result holds in cosmological models that expand for all time, rather than re-collapsing, remains to be seen.
It should be noted, however, that while the effects we find are always small (less than 1%, in most cases), the Taylor series approximation used to derive them usually breaks down on time scales that are shorter than the age of the universe. This is especially true in the larger lattices, where it can be seen that increasing the number of masses in the universe results in the series approximation breaking down at earlier cosmological times. This result holds for all of the LRS curves that we studied, but seems to be especially true of the cell edges (which, unfortunately, are probably the best indicators of the scale of the cosmology as a whole). To reliably follow the evolution of these curves any further will require more advanced techniques.
It is interesting to see that the effect of higher-order terms in Einstein's equations can lead to non-negligible effects in the large-scale expansion of space. One may note, for example, that it is not until we get to order t 6 in the Taylor series expansion in Eq. (74) that the influence of curlH ab becomes non-zero at all. In a weak-field expansion of the gravitational field, about a suitably chosen background, this would probably correspond to quite a high order in perturbations. Nevertheless, the terms involved grow rapidly over cosmological time-scales, until they come close in magnitude to the leading-order terms represented by E ab . This is remiscent of the gravitational wave memory effect [22], where the small effects of non-linear gravity accumulate over time until (in the case of gravitational waves from astrophysical sources) they are comparable with the magnitude of the linear, leading-order terms.