Consistent cosmological structure formation on all scales in relativistic extensions of MOND

General relativity manifests very similar equations in different regimes, notably in large scale cosmological perturbation theory, non-linear cosmological structure formation, and in weak field galactic dynamics. The same is not necessarily true in alternative gravity theories, in particular those that possess MONDian behaviour (“relativistic extensions” of MOND). In these theories different regimes are typically studied quite separately, sometimes even with the freedom in the theories chosen differently in different regimes. If we wish to properly and fully test complete cosmologies containing MOND against the ΛCDM paradigm then we need to understand cosmological structure formation on all scales, and do so in a coherent and consistent manner. We propose a method for doing so and apply it to generalised Einstein-Aether theories as a case study. We derive the equations that govern cosmological structure formation on all scales in these theories and show that the same free function (which may contain both Newtonian and MONDian branches) appears in the cosmological background, linear perturbations, and non-linear cosmological structure formation. We show that MONDian behaviour on galactic scales does not necessarily result in MONDian behaviour on cosmological scales, and for MONDian behaviour to arise cosmologically, there will be no modification to the Friedmann equations governing the evolution of the homogeneous cosmological background. We comment on how existing N-body simulations relate to complete and consistent generalised Einstein-Aether cosmologies. The equations derived in this work allow consistent cosmological N-body simulations to be run in these theories whether or not MONDian behaviour manifests on cosmological scales.


Introduction
The current cosmological paradigm, ΛCDM is broadly successful [1], although it is has some ongoing small scale problems and tensions (see e.g. [2][3][4][5]). This paradigm contains two hypothesised forms of matter: cold dark matter and dark energy. The existing observational evidence for these is solely gravitational in nature, so it is natural to consider whether modified gravitational laws (rather than Einstein's General Relativity (GR)) could instead be responsible for these observations, see e.g. [6] for a comprehensive review. A long-running suggestion to account for some of the phenomenology attributed to dark matter is MOdified Newtonian Dynamics (MOND) [7][8][9]. See [10][11][12][13] for reviews of the theory and observational successes and challenges of MOND.

JCAP06(2023)006
list of recent topics and discussions). In this work we do not take a position advocating for or against MOND. Rather, we take the view that until there is non-gravitational evidence for dark matter, it is worth testing and developing the different paradigms in an attempt to falsify each of them. However, to do this we need to understand cosmological structure formation on all scales 1 in a MONDian cosmology and there is as yet no complete and consistent MONDian cosmology that makes clear predictions on all cosmological scales and that can be fully and fairly tested against the ΛCDM paradigm. Partly this is difficult because MOND itself is not a theory of gravity, it is a phenomenological description of how the gravitational laws might behave in certain limits. To be a complete theory, MONDian behaviour must be embedded in a "relativistic extension", several candidates for which have been investigated [24][25][26][27][28][29]. See references in [28] for a more comprehensive list of different theories in which MONDian behaviour can arise in different ways.
An additional difficulty is that studying cosmological structure formation requires looking at several different regimes, none of which are the weak-field galactic dynamics regime for which MOND was originally proposed. MONDian behaviour is often expressed as a modification of the Poisson equation 2 in this regime, but here we run into the issue that the Poisson equation is something that is a little convenient in GR+ΛCDM: a very similar equation arises in conceptually quite different regimes, notably large scale cosmological perturbation theory, non-linear cosmological structure formation, and in weak field galactic (and smaller scale) dynamics. The equations in these regimes are not necessarily so similar in other theories of gravity, and this particularly applies to the case of relativistic extensions of MOND. These theories are usually studied piecemeal in different regimes, with different assumptions and frameworks in each, and sometimes even with the available freedom in the theories being chosen differently in different regimes. This disjointed approach is theoretically problematic as it is unclear if the different choices made in different regimes can be simultaneously realised in the same universe. It is also practically problematic since it means it is unclear how to run cosmological simulations when the different limits are treated so differently, due to the range of scales that is being covered by these simulations. In particular these simulations need to have a cosmological background and large scale perturbations that are consistent with each other, and with the smaller scale non-linear behaviour. Studying cosmological structure formation in GR is more straightforward because of its similarity across the different regimes.
Cosmological N-body simulations up to scales of 512h −1 Mpc have previously been run with some form of MONDian Poisson equation [30][31][32][33][34][35][36], however due to the issues noted above it is somewhat unclear how these simulations relate to specific relativistic extensions of MOND. From the other side, the equations that govern cosmological structure formation on all scales have not been derived for any relativistic extensions of MOND; in particular, equations that can be used to run cosmological N-body simulations for a specific theory. Another use of such equations would be to examine on which scales involved in cosmological structure formation MONDian behaviour can arise, and if MONDian behaviour on these scales necessarily arises from having MONDian behaviour on galactic scales.
In order for the MOND paradigm to develop into sufficient maturity to be fully and consistently compared against ΛCDM, these issues need to be addressed. In this work we 1 Note that throughout, we use "cosmological structure formation on all scales" to mean all scales where a perturbed FLRW metric with weak fields is a reasonable description of the spacetime (a range of scales from super horizon scales to approximately scales of order Mpc), and by "non-linear scales" we mean scales of around 10Mpc and below, where the cosmological density contrast δ becomes greater than 1.
2 Note that in this paper we consider MOND as a modification of gravity not as modified inertia.

JCAP06(2023)006
examine these issues using the post-Friedmann approach [37,38]; this is a weak field post-Newtonian-like expansion designed to work in a cosmological setting, so it is ideal for such a study. In addition to GR+ΛCDM, it has previously been applied to Hu-Sawick [39] f (R) gravity [40] and used to create a model-independent approach to modified gravity that applies on all cosmological scales [41,42]. We create a method for deriving the equations that govern cosmological structure formation (the evolution of inhomogeneities) on all scales that is inspired by the observations in [41]. In particular, the process used to derive the all-scales equations partly draws on the conditions in GR+ΛCDM that cause the simulation situation to be relatively simple, such as the absence of a regime where neither perturbation theory or the Newtonian limit apply (as discussed in [41]). We apply this method to generalised Einstein-Aether (GEA) theories [43][44][45][46][47][48][49][50] as an illustrative example of how a relativistic extension of MOND can be examined more holistically over the full range of cosmological scales. This process gives a single set of equations for running consistent cosmological simulations, including an expansion history, large scale behaviour, and small scale behaviour that are consistent with each other. We examine how MONDian behaviour can arise in the resulting structure formation and if and how this relates to MONDian behaviour on galactic scales, as well as what this means for relating existing cosmological MOND N-body simulations to GEA cosmologies. We also examine how one can check the assumptions that underly our derivation. This paper is laid out as follows: in section 2 we elaborate on the theoretical context and briefly recap some details of GEA theories and the post-Friedmann approach. In section 3 we apply our approach to GEA theories to construct equations that describe cosmological structure formation on all scales in these theories. We discuss some features of these equations in section 4, notably if and how MONDian behaviour can arise cosmologically. We conclude in section 5.

Theoretical context
In this section we discuss some of the issues around MOND N-body simulations, relating them to relativistic extensions of MOND, and why we need equations for these theories that apply on all scales. We then briefly recap GEA theories (mostly following [43,45]), and the post-Friedmann approach (mostly following [37,41]), and in the latter we explain the method that will be used to derive the equations for structure formation that apply on all cosmological scales.

(MOND) N-body simulations and the need for coherent and complete cosmologies on all scales in relativistic extensions of MOND
The equations governing cosmological structure formation in ΛCDM+GR are very similar in the linear perturbation limit and in the non-linear Newtonian limit, and there aren't any scales where the leading order dynamics are not well described by one of these two limits (see e.g. the discussion in [41] and references therein). As a result, it is relatively straightforward (from a physics, if not a computational perspective) to run N-body simulations that cover a wide range of scales, from super horizon scales all the way down to scales where the density fluctuations are large. The situation is typically more complicated in modified gravity theories, where the different assumptions and properties of the different regimes means that the equations can differ more between the different regimes. This issue is particularly pronounced for relativistic extensions of MOND for two reasons. Firstly because different limits of the theories are  Schematic showing the different conceptual regimes referred to throughout the paper, the scales they cover relative to each other, and how this relates to MOND and cosmological simulations. The Newtonian scale denotes the scale below which the Newtonian limit of a relativistic theory is a good approximation to the full equations. Note that there is no intermediate regime in ΛCDM: in this case the perturbation and Newtonian regimes overlap. Similarly, the perturbative quasi-static regime typically only exists in cosmologies with no intermediate regime.
typically studied independently, which sometimes even includes choosing the undetermined freedom in the theories differently in different regimes. 3 The second is that the cosmological regimes have different assumptions to the galactic limit in which MONDian behaviour is shown to arise in these relativistic theories. To the authors' knowledge, a single consistent description covering all of the scales for cosmological structure formation has not previously been derived for a relativistic extension of MOND. Instead, behaviour is typically extrapolated from galactic scales into the non-linear regime of cosmological structure formation (as discussed in the next paragraph); such extrapolations involve moving between regimes in which particular limits of the theory are studied and the freedom is chosen in a particular way. To aid the reader, a schematic illustration of the different scales referred to in this work is show in figure 1. This figure shows the different cosmological and non-cosmological regimes, as well as the range of regimes spanned by cosmological N-body simulations. In this work for simplicity we assume that there is a lower limit to the range of scales on which the use of the FLRW metric, and thus the derivation later in this paper, is valid. 4 As such, the backgrounds of the two Newtonian limits differ, and both the additional time dependence of the FLRW background and that it is not a vacuum solution (but Minkowski is) could create differences for other theories of gravity.
From the simulation side, since MOND is often taken to be a modification to the Newtonian Poisson equation, one common route is to modify N-body simulations to include the MONDian Poisson equation. This is valid on galactic scales where a simple weak field (in a Minkowski background) derivation of MONDian behaviour can be carried out on a par- 3 For the Generalised Einstein Aether theories considered later in this work, this freedom is the function F (K) and its form for different ranges of K. 4 This means that we always work within a cosmological context and do not consider for example whether one can always move to FLRW co-ordinates even on smaller (e.g. galactic) scales as carried out for example in [51,52]). If such transformations are possible in the theories we consider here then this work can be extended to include galactic behaviour and constrain these theories further from the point of view of requiring consistent behaviour on all scales.

JCAP06(2023)006
ticular relativistic extension of MOND. Here, we focus on cosmological-type simulations on larger scales, rather than simulations on smaller scales focused around simulating galaxy formation and evolution. Previous works [30][31][32][33][34][35][36] have been run with box sizes varying from 32h −1 Mpc to 512h −1 Mpc, all of which involve scales that will be evolving according to linear and non-linear fluctuation behaviour in a cosmological background.
With "cosmological" box sizes such as these, it is less clear cut whether MONDian behaviour necessarily arises even in theories that have MONDian behaviour on galactic scales, because cosmological structure formation is a different conceptual regime. Again, just because the GR equations are similar in the different regimes, doesn't mean that this will be the case for other theories. This question has not previously been addressed. 5 This issue is compounded by the aforementioned problems about the different cosmological regimes that are spanned by an N-body simulation. For example, an (effective) modified G in the Poisson equation on linear cosmological scales is a generic prediction of modified gravity theories (see e.g. [6]), so such an effect may be required on the larger (linear) scales in MONDian cosmological simulations. Such a modification in linear theory needs to be understood in terms of how it relates to the branches of the MOND Poisson equation in non-linear structure formation, whether it arises in one or both of these branches, and therefore how it should arise in a cosmological N-body simulation.
For concreteness, when we refer to galactic MONDian behaviour in this work we are referring to a Poisson equation given by where ρ is the matter density and Φ is the Newtonian potential. Since MOND was proposed in a galactic context, when applying it cosmologically one has to make a choice about whether the a 0 function is allowed to vary over cosmic history, and whether the quantities that appear in the MOND Poisson equation (such as ρ, Φ and ∇) are the physical or comoving quantities. 6 We make a small generalisation of MONDian behaviour to include an extra function (γ(a)) of the scale factor a to account for these different possibilities and choices Behaviour matching equation (2.4) will be the behaviour we refer to as cosmological MOND behaviour when examining a specific relativistic extension of MOND later in the paper. We note that this simple form does not quite encompass all of the possible choices, but it covers the main ones and is reasonably general whilst maintaining a fairly simple form. This form is similar to that adopted for cosmological MOND N-body simulations (e.g. [32,36]). 5 Although we note that careful thought is given in [30][31][32][33][34][35][36] to the fact that the underlying covariant cosmology is not known, requiring assumptions such as the validity of the Friedmann equations for the background expansion, assuming that MONDian effects only apply to the peculiar acceleration, and the validity of the initial conditions that are typically used in ΛCDM N-body simulations. Since the method in this paper creates a consistent cosmology on all scales and times, these questions will be effectively answered by the same framework, although these are not the primary questions we are investigating. 6 Note that the subscript p on ∇ in equation (2.1) denotes that this is a derivative with respect to a physical co-ordinate; ∇ throughout the rest of this paper is a derivative with respect to a comoving co-ordinate.

JCAP06(2023)006
To summarise, the following are still open questions: what does a complete picture of cosmological structure formation look like in relativistic extensions of MOND, how do we run consistent cosmological N-body simulations for these theories, how do said simulations relate to existing cosmological MOND simulations, and does MONDian behaviour necessarily arise on cosmological scales if it arises on galactic scales (and if it does appear, then in which scales, times and locations does it appear)? These issues apply to analytic examinations of structure formation in MOND as well (e.g [53,54]): if a particular relativistic extension of MOND contains MOND on galactic scales but not cosmological scales, then the general lessons drawn from studying cosmological structure formation in MOND will not apply to these theories; conversely cosmological structure formation potentially being found to not be compatible with these studies does not necessarily therefore rule out MOND phenomenology on galactic scales.
To solve these issues in a particular theory, we need to derive equations that describe the evolution of cosmological inhomogeneities on all scales. In this work, we propose a method to do so and illustrate it by applying it to GEA theories as a case study.

Generalised Einstein Aether theories
GEA theories are theories with an additional timelike vector field, A µ , with an action given by where λ is the Lagrange multiplier ensuring that the Aether field A µ has a timelike unit norm and the function F (K) is free and can be chosen to get different phenomenology. The scalar K is given by The c i are constants, and are constrained to be c 1 < 0, c 2 ≤ 0 and c 1 + c 2 + c 3 ≤ 0. For later convenience we also define α = c 1 + 3c 2 + c 3 . There is also a possible 'c 4 ' contribution K µν αβ = c 4 A µ A ν g αβ , see [49] for a discussion. For simplicity, we drop this term and concentrate on theories that have been shown to deliver MOND on galactic scales [43,45]. The method presented later can accommodate non-zero c 4 if required. For these theories, the Einstein equations are given by µν being the stress energy tensor of matter. The vector field equation is given by GEA theories have been studied extensively in terms of their cosmological background and linear perturbations [43,[45][46][47][48][49][50], and also in terms of their gravitational wave propagation [55,56]. These theories have also been looked at in the weak field limit [43,45] that describes the Solar System and galaxies. Notably, the free function F (K) is usually chosen differently in these different regimes in order to get MONDian behaviour in galaxies and also get dark-energy-like cosmological backgrounds. This choice is made on the assumption that there is a single (complicated) F (K) that has different behaviour for different ranges of K, and the ranges of values that K takes in the different regimes are disjoint. However, the validity of this assumption is less clear if we include (non-linear) cosmological structure formation in which MONDian behaviour may arise on scales that are described by a FLRW background. This is one of the issues that our approach will resolve. Note that we are not positing that GEA theories themselves are necessarily strong contenders to ΛCDM, merely using them as a reasonably well-studied example to illustrate how our method can be used to derive complete and consistent cosmologies in (and insight into) relativistic extensions of MOND.

Post-Friedmann formalism
The post-Friedmann formalism is a post-Newtonian-like expansion of the Einstein equations (or modified Einstein equations) in powers of the speed of light c, altered compared to a "Solar-System" type expansion in order to apply to a FLRW cosmology [37,38]. The starting point of the post-Friedmann approach is the perturbed FLRW metric in Poisson gauge, 7 which is expanded up to order c −5 , The two scalar potentials have each been split into their leading order (Newtonian) (U N ,V N ) and higher order (U P ,V P ) components. The gauge freedom is chosen such that the vector potential appears in the 0i part of the metric, and this has also been split up into B N i and B P i . The three-vectors B N i and B P i are both divergence-less, B N i,i = 0 and B P i,i = 0. In addition, the tensor perturbation h ij is transverse and trace-free, h i i = h ,i ij = 0. Time derivatives are associated with a factor of 1 c . The matter content (in addition to a possible cosmological constant) is taken to be pressure-less dust, 8 the four-velocity of which is used to construct the energy-momentum tensor, which is also expanded in powers of c. The parameters describing the pressure-less dust fluid are the background densityρ, the density contrast δ, and the peculiar velocity v i . Crucially, this setup doesn't require the density contrast to be small, unlike standard (linear) cosmological perturbation theory, which allows the approach to be used on small scales. To obtain the equations that apply on all cosmological scales we will JCAP06(2023)006 change variables to the "re-summed potentials" [37,41] For full details of the post-Friedmann formalism, see [37].

Cosmological structure formation on all scales
The post-Friedmann formalism was used in [37] to construct equations that apply on all cosmological scales, including both the large and small scale limits and a possible "intermediate" regime where neither linear perturbation theory nor the Newtonian limit 9 is sufficient for calculating cosmological structure formation. These equations are quite complicated, however they can be simplified [41] under the assumption that there is no intermediate regime, i.e. that all scales of interest are described by either linear perturbation theory or the cosmological Newtonian limit. This assumption holds for ΛCDM and at least some modified gravity theories [40,41]. The simplified equations describe the evolution of cosmological homogeneities on every scale from super-horizon scales to the regime of non-linear cosmological structure formation as long as the assumption that there is no intermediate regime holds. For a ΛCDM+GR cosmology, these simple all scales equations are given by equations (3.20a) − (3.20d) in section IV of [41] .
In this work we build on the idea of [41] to present a method for deriving equations that apply on all scales in a specific theory of modified gravity. In [41] the simplified all scales equations are arrived at by expanding to c −5 order, converting to the re-summed potentials, and then removing certain terms from the equations: specifically the terms that are both "structurally non-linear" and beyond leading order in the 1 c expansion. These same equations can be arrived at with a simpler method, namely for each equation derive the two limiting cases, write them in terms of the re-summed potentials, and then combine the two equations into a single equation that contains both limits (and no more). This is the method we will follow here, and it results in a single consistent set of equations that describe the background, and fluctuations on all cosmological scales in these theories. By definition this is equivalent to constructing the full set of master equations up to order c −5 , converting to the re-summed potentials, and removing the terms that associated only with an intermediate regime where neither the linear perturbation theory or Newtonian limit applies, however it is a much easier computational process. This process amounts to "stitching together" the two limits, which works because they are both weak field expansions and thus can both be written in terms of the re-summed potentials.
This method is not universally valid for any cosmology or matter content. By construction, this derivation only applies in a matter or Λ dominated cosmology, i.e. the dominant clustering component in the universe needs to be pressureless dust. Since the method followed here yields the same equations as the method laid out in [41], the remaining requirements for this method to be valid are largely given by the set of criteria in section V.A of [41]. These requirements can be approximately summarised as follows: only scalar metric potentials are important for cosmological structure formation, the two limits (linear perturbation and Newtonian) are valid, and every relevant scale is well-described by one (or both) of the two limits (in terms of figure 1 this means that there is no "intermediate regime" and thus that the Newtonian and perturbative limits overlap). The first of these is not required in the GEA context. We will examine the remaining requirements for the GEA case in section 3.4, here we just note that in any cosmological N-body simulation, the vector potential ω can be used as a quantitative test of the Newtonian limit assumption [40,41,58,59]. As such, this test is one of the tests for this method, and the required equation to perform this test is naturally derived as part of the method; in the GEA case this is equation (3.26).

Derivation of all-scales equations in GEA theories
As described in the previous section, this derivation has two parts: finding the leading order post-Friedmann equations, and then combining them with the linear perturbation equations.

Applying the post-Friedmann approach to GEA theories
To begin, we need to decide how to expand the vector field A µ in the post-Friedmann 1 c expansion. We do this as follows and for A µ it follows that: where β i = δ ij β j . Hence the homogeneous background part of the aether field is given byĀ µ = (1, 0, 0, 0) as usual. Our choice of expansion for the perturbations is justified by looking at the Lagrange multiplier equation (ensuring A µ is timelike) order by order. It is fairly straightforward to expand the perturbation to A 0 : it typically behaves like a scalar field, which are usually expanded in even powers of c and, from the Lagrange multiplier equation, in the weak field limit it is expected to be equal to U N . Both of these considerations are satisfied by the leading order perturbation to A 0 being of order c −2 . Assuming that the spatial perturbation to the vector field is not of lower order in c than the time perturbation, then the Lagrange multiplier equation at c −3 order sets any possible c −3 contribution to the time perturbation to be zero. Then the only question is the spatial perturbation: the Lagrange multiplier at c −5 order requires where β i(n) refers to the part of β i that comes with an attached power of c n . If β i(2) = 0 then one would see that the leading order terms of the spatial part of the vector equation force this term to either be zero or satisfy a strange constraint, 10 so we choose β i(2) = 0; the 10 c1β (2)j ,ij = 0.
It is not surprising for the pure vector part of the spatial perturbation that there is a c −3 factor, as the pure vector metric potential is also c −3 .
In order to be as general as possible with the freedom allowed in GEA theories we do not specify a form of F (K) here, or assign a power of c to M in the derivations. Instead we initially assume that K, F (K) and F K can be sensibly expanded, and therefore in these equations M 2 F (K) and F K terms are taken to represent the leading order part. Terms denoted F (K) orF K mean only the homogeneous part of any such leading order term. In section 3.3 we examine the issue of assigning a power of c to M and what this means for the expansion and the resulting equations, as well as how this relates to possible choices of F (K). We then consider specific F (K) that might lead to MONDian behaviour in section 4.
With these choices we can now calculate the leading order equations. We include calculations of some of the intermediate quantities in appendix D for the interested reader and to ensure that others can reproduce our results if required. The Lagrange multiplier is calculated from the time part of the vector field equation (throughout, bothẋ and ∂ T x denote differentiation with respect to time) and K is given by The homogeneous part of T (GEA) µν (the GEA contributions to the Einstein equations) at leading order (up to c −3 ) is given bȳ and the inhomogeneous parts are given by The homogeneous equations are given by

JCAP06(2023)006
and the full (inhomogeneous) equations with the homogeneous parts subtracted are given by The spatial part of the vector field equation is We can re-write the inhomogeneous equations in terms of the re-summed potentials, split the spatial perturbation to the vector field into its scalar and vector parts as β i → ∂ i ξ S + ξ i , and simplify the equations to get The scalar K re-written in terms of the same variables is given by Equations (3.12)-(3.16) represent the Newtonian limit of the equations governing cosmological structure formation, and apply regardless of the size of the density contrast. They can be compared to both the cosmological perturbation theory equations and galaxy-scale weak-field equations to understand the phenomenology of these theories and how the different regimes relate to each other. We will use these equations with the process laid out earlier to derive equations that govern cosmological structure formation on all scales.

All scales equations
The leading order 1 c equations and linear perturbation equations (the latter are in appendix A) can be combined to give -12 -

JCAP06(2023)006
We have denoted by LS the large scale terms that only apply on scales close to the horizon, the explicit form of which are given in appendix C. These terms can usually be neglected in N-body simulations (but do not need to be), see discussion in [41]. Their inclusion here is still important as it ensures that we have a single set of equations that describe structure formation on all scales, containing the same variables and with the freedom in the theory chosen consistently. In other words, when examining the possibility of cosmological MON-Dian behaviour in section 4, we can ensure that the behaviour in the different cosmological regimes are all consistent with each other and realisable in the same universe, with the same choice of F (K) (including any conditions on K for different branches of F (K) to be realised).

Consistently expanding M 2 and F (K)
We now consider the issue of M 2 in the expansion in powers of c, and how this relates to expanding F (K) and F K . By construction, K is dimensionless. However M 2 is not, and dimensionful constants typically need to be considered carefully in 1 c type expansions (see e.g. [41,60]). For the original linear EA theory M 2 K, it makes no difference to the leading order equations which power of c is assigned to M , however different choices will result in difference outcomes for the more general case with F (K). The following considerations and observations will be used to guide our choice • The power of c attached to the leading order of F K will differ from that of the leading order of F (K) by the power of c attached to the leading order of K.
• It makes little sense for M 2 F (K) and F K to contribute at lower order than the standard GR and matter terms, thus requiring M 2 F (K) to be of order c −2 or smaller and F K to be of order c 0 or smaller. • For any power of c assigned to M , ifK = 0 then the choice of the functional form F (K) is much more constrained in order to deliver sensible results. For example, a power law F (K) = AK n is problematic for the 1 c expansion if n < 0. The same is probably true for a Taylor expansion of F (K) and F K as typically performed in perturbation theory: ifK = 0 then this expansion is problematic.
Given all of these, the choice M 2 → M 2 * c 2 appears the most sensible, at least as long as there is a background for K (i.e. α = 0). With this choice, the leading order part of K purely homogeneous (i.e. the leading orderK is c 2 larger than the leading order of δK), which makes sense in terms of expanding around a background and the dimensionless nature of K. It also makes the outcome more similar to a perturbative expansion, so many of the terms that could exist in principle in the non-linear regime are not present at leading order. The contributions at leading order become for a powerlaw (K n )

JCAP06(2023)006
i.e. the homogeneous parts contribute at exactly the expected order for all power laws for F (K), and the inhomogeneous parts won't contribute at leading order. There are few explicit forms for F (K) beyond power laws in the literature that we can use to test our choices further. One example is [61], which uses F (K) ∼ √ K + √ K ln K. For the choices we have made, this function also gives a sensible expansion as long as α = 0. With this choice, the equations become where the large scale terms are denoted by LS and their explicit form is given in appendix C. Note that not all of the terms that get removed from the leading order equations as a result of applying this choice for M appear in the revised "large scale terms", since some of the removed terms no longer contribute at leading order in either limit, and thus they can be removed entirely. We also note that when examining MOND phenomenology in section 4 the unusual properties of the specific F (K) that is used necessitates keeping some of the sub-leading terms that have been neglected in deriving this second set of equations. More generally, if one is concerned about the properties of a particular F (K) then one can always insert M 2 → M 2 * c 2 into the first set of all-scales equations and check whether any additional terms should be kept.
Taking the vector part of the spatial vector field and G 0i equations (and ignoring the large scale terms) we find which will be used in the next subsection.

Conditions for the derivation to hold
Here we paraphrase section V.A of [41] to formally lay out the conditions that are required for this derivation to be valid, and thus have been implicitly assumed above. Since the quantitative tests require detailed numerical simulations, we do not carry the tests out here, -14 -

JCAP06(2023)006
however we lay out these conditions to give a full picture of how the method works conceptually and how to test the assumptions that underlie the method. These assumptions are often explicitly or implicitly in other calculations carried out in modified gravity cosmologies (for example when using N-body simulations).
As noted earlier, we don't require the first condition in section V.A of [41] that only scalar fluctuations are important for structure formation, so the remaining requirements are equivalent to the two limits (linear perturbation and Newtonian) are valid on the largest and smallest scales, and every relevant scale is well-described by one (or both) of the two limits, i.e. there is no intermediate regime (see figure 1) where neither of the limits apply and there is at least a small region of overlap where both limits apply. As such, these requirements can be written more concretely as 1. A weak field metric is appropriate on all cosmological scales.
2. Check for the existence of a scale k * , which is between the horizon scale and the scale at which the density fluctuations become non-linear, such that on length scales below the length scale corresponding to k * the only terms that contribute significantly in the linear perturbation equations are the same terms that are present in the linearised Newtonian equations.
3. Calculate the metric vector potential ω from GEA N-body simulations using equation (3.25), and check that this is small enough on all non-linear scales for the Newtonian approximation to be valid on all of these scales.
The first of these is assumed to be axiomatic in most cosmologies examined in the literature (including ΛCDM and modified gravity cosmologies) and there is no evidence to the contrary (see e.g. the discussion in [41]), so we do not discuss it further. The second and third conditions on this list are quantitative tests of the validity of this derivation that require numerically solving the equations.
The second test requires using a modified cosmological Boltzmann code that includes the GEA modifications to the linear perturbation equations, and in practice this is equivalent to testing the quasi-static approximation (see discussion in [41]). The quasi-static approximation is a fairly standard approximation in modified gravity, and is typically assumed to hold in GEA theories (e.g. [48]), so we assume the same here and leave an in-depth examination of this to future work.
The hardest condition to test is the third one, namely that the Newtonian limit of the gravitational equations is a good approximation on all scales below the scale at which density fluctuations become non-linear. This condition has been checked explicitly for GR+ΛCDM [58,59] and f (R) [40], and is implicitly assumed to be true whenever N-body simulations are run for a particular cosmology, so it can be seen as a formal check of the working assumptions under which simulations are run, irrespective of the issue we consider in this paper of how to consistently describe structure formation on all scales. This test requires N-body simulations for the cosmology in question, which can now be run for GEA theories using the equations laid out earlier.

Discussion and phenomenology of the equations
Equations (3.21)-(3.24) are the GEA equivalent to the ΛCDM+GR equations (3.20) in section IV of [41]. These gravitational equations should be combined with the matter equations -15 -

JCAP06(2023)006
(3.20c) and (3.20d) from [41] (repeated in appendix B for completeness) in order to provide a complete and unified set of equations for evolving the density fluctuations in the universe on any scale in a generalised Einstein-Aether cosmology. Here we examine these equations and their consequences.
One key point of these equations is that since cosmological structure formation becomes non-linear in the matter dominated era, and by construction the all-scales equations smoothly connect to the perturbation theory equations in the matter dominated era, then the combination of the equations here (describing the matter dominated universe onwards) and the (well-studied) cosmological perturbation theory for GEA theories comprises a complete and coherent description of GEA cosmologies. To the authors' knowledge this is the first complete and coherent cosmology of a relativistic extension of MOND covering all necessary scales and regimes. These equations also contain the source equation for the vector potential in the metric, which can be used to check the validity of the Newtonian limit on all non-linear scales [41,58,59] and can generate novel lensing signals, so it is a useful equation to have derived for GEA cosmologies irrespective of its uses in the context in this work.
We note that in these equations the same F (K) (with the same branches and conditions for the branches to be realised) is present on all cosmological scales, and therefore the chosen function will govern all of the cosmological dynamics, including the background expansion, large (linear) scales and non-linear scales where the density contrast is large. For consistency of the theory and the cosmology, this F (K) cannot be chosen differently in these different regimes (unless the branch conditions are such that a different part of F (K) applies in each regime, which is unlikely numerically). As a result, we can see that introducing MONDian behaviour on non-linear scales is likely to also introduce MONDian behaviour into the cosmological perturbation equations (as long as the numerical values of a 0 and the gravitational potential are such that the perturbations are in the MONDian regime, or at least in transition regime between the Newtonian and MONDian branches). In other words these two regimes are intertwined and one cannot arbitrarily put MONDian behaviour into only one of the two regimes. We will see this in more detail in subsection 4.2.
In the weak field limit that is present in galaxies, deep MONDian behaviour arises using lim |∇Φ| a 0 F (K) = − 2 c 1 K + BK 3/2 [43]. When we apply this F (K) to our cosmological equations, we find that in general MOND behaviour does not arise, due to the presence of the background terms in K, so we examine this choice of F (K) separately for the two cases and α = 0 and α = 0. We start by assuming we are in the deep MOND regime in each case, and then later comment on choosing F (K) for the Newtonian branch and the criteria determining the different branches for each case.

α = 0
Using K as in equation (3.16) and F (K) = − 2 c 1 K + BK 3/2 , we find that the leading order (in 1 c ) G 00 equation is given by The contributions that are usually responsible for delivering MONDian behaviour are higher order c −3 in this case, so they do not contribute at leading order. There are no issues performing the perturbative expansion in this case, and the perturbative equations match the 1 c equations in the overlap regime, so the all scales equations follow from substituting the -16 -

JCAP06(2023)006
appropriate K and F K into the equations in section 3.3. We do not make this substitution here for brevity. It is often taken that M * ∼ H 0 ; in this case there is an order unity change to Newton's constant in the cosmological Poisson equation, with no scale dependence and only simple time dependence. Thus, GEA theories would not manifest any cosmological MONDian behaviour if the background expansion is modified from the GR Friedmann equations. These theories might still have MONDian behaviour on galaxy scales below where the FLRW metric applies. As MONDian behaviour does not arise for this case, we do not look further into the issue of defining different branches or choosing different F (K)s in different branches.

α = 0
We now consider the case α = 0, such that K = −c 1 φ P,k φ ,k P M 2 * a 2 c 2 . This is a somewhat unusual K in thatK = 0, but F K contains a constant (and thus effectively homogeneous and not perturbatively small) part for any K. As such, a Taylor expansion of F (K) 11 (as sometimes used perturbatively) fails. Truncating the G 00 equation at order c −2 gives 8πGρδ c 2 = 0, so the leading order equation (in the c −1 expansion) must be expanded to include additional terms of order 1 c 3 or higher. The K 3 2 part of F (K) does not contribute at order c −2 , but it does contribute at c −3 order through the term − c 1 c 2 a 2 (φ P,i F K ) ,i in equation (3.17). The other beyond-leading-order terms arise at c −4 or higher order, or are of order 1 c 3 but occur with a prefactor α, so disappear for this case, leaving As such, the only terms at this order are exactly the terms that give rise to a MONDian Poisson equation as in the galactic weak field case, therefore in this case it is possible to get MONDian behaviour on cosmological scales. The unusual nature of the choice of F (K) is what makes the derivation of MOND behaviour in this way in a relativistic extension of MOND somewhat strange, in that F K cancels the leading order term from GR, and replaces it with a term that is nominally not leading order, but is actually the dominant term that remains after cancelling the usual leading order terms. An analogous cancellation occurs in the weak field galactic MOND case. Allowing for this case in the general G 00 all-scales equation requires only a very minor adjustment (reverting one of the terms from equation (3.21) to its form in equation (3.17)), This equation replaces equation (3.21) when it is desired to include GEA theories with cosmological MONDian behaviour. One of the advantages of the all-scales equations we have constructed is that the leading order (in the c −1 expansion) parts of these equations cover the quasi-static linear and nonlinear regimes consistently and simultaneously. As such, we can see how this MONDian behaviour manifests on scales that are well below the horizon where the density contrast is small: this is what is usually referred to as the (perturbative) quasi-static regime. If one naively applied the quasi static approximation and linear perturbation theory in this regime, 11 I.e. expanding as F (K) → F (K) + FK (K)δK and FK → FK (K) + FKK (K)δK. then one would have a nonsenical leading order equation and no MONDian behaviour, despite the Newtonian limit of the (non-linear) equations suggesting that such a MONDian term should be present and dominant. The combination of the all-scales equations and the tests to determine their validity gives a way out of this problem. Firstly one can run simulations with the leading order 1 c equations and determine what scale these are valid up to (using the vector potential check described earlier) and what size the density contrast is on those scales. If the 1 c equations are valid up to a scale where the density contrast is small (linear), then the problem is the linearisation in perturbation theory preventing the MONDian term. One can then use the linearised all-scales equations in a Boltzmann code and check that the quasi static approximation is valid to a larger scale than the non-linear scale (from the perspective of perturbation theory). If this second test is passed then the all-scales equations are valid, despite the problem with a naive application of perturbation theory. If either test is failed, then this means the cosmology in question, unlike ΛCDM, has a range of scales where neither perturbation theory nor the Newtonian limit is valid. The authors are not aware of any tools for computing cosmological structure formation in any cosmology that has such a regime. Relatedly, since the all-scales equations do not apply to a radiation-dominated cosmology, it remains unclear if and how cosmological MONDian behaviour might arise in the perturbation equations during the radiation dominated era.
As deep MOND behaviour can manifest for this case, we now consider how a Newtonian branch could behave, and what the conditions for being in different branches would be. 12 The condition to be in the deep MONDian branch of equation ( where this identification is independent of any power of c assigned to M . We have observed from the first equation that for GEA theories manifesting cosmological MOND, γ(a) = a, however we leave γ(a) in the equations for the rest of this subsection to ease comparison with the deep MONDian branch of equation (2.4). Interestingly, this form of γ(a) means that when converted to a physical derivative, the condition to be in the deep MONDian branch, and thus the effective a 0 , does not vary with time. This is similar to the behaviour found in a different relativistic extension of MOND in [62]. We can substitute into the expression for K and define K * = 16c 4 to relate K to the definition of the deep MOND branch As such, the MONDian and Newtonian regimes are defined by how K compares to the (constant) critical value K * , with some explicit time dependence in the relationship. The

JCAP06(2023)006
MONDian branch and MONDian behaviour occurs when 0 < K K * , and the Newtonian branch occurs when K K * . In principle, the function F (K) can be chosen to have a different form in the Newtonian branch, however despite this the Newtonian branch is significantly constrained, as the term that modifies the Poisson equation is proportional toF K ∇ 2 φ P , which will not contribute for most F (K) if α = 0. As such, requiring cosmological MOND behaviour in GEA theories means that the Newtonian branch in the quasi-static limit of perturbation theory is the GR Poisson equation, i.e. there is no modification to gravity and G eff = 1 in this branch, including in linear perturbation theory.
There is one final technical nuance to consider, which is whether the emergence of MONDian behaviour, the identification with a 0 and the conditions for the different branches are consistent with how we have derived our equations and the power of c assigned to M . Assuming no power of c assigned to M , the leading order term that we have neglected is (4.8) The As an additional consistency check, we can also estimate the effect of the first term in the galactic weak-field regime by noting that it is akin to adding a term ∇Φ · ∇Φ on the left hand side of 2.1. It can be shown that this modifies the solution Φ in the deep MOND branch at a = 1 from ∂ r Φ = √ GM a 0 /r to ∂ r Φ = e −a 0 r √ GM a 0 /r so its effect is subdominant for r a 0 . Typical values for a 0 correspond to length scales of the order of the Hubble radius and so the first term is expected to be subdominant for the subhorizon scales of interest.
Since ∇φ P is free to vary substantially numerically just as in a usual MOND weak field expansion, the conditions for the branches of the MOND Poisson equation (|∇φ| γ(a)a 0 and |∇φ| γ(a)a 0 ) can both be realised without violating any assumptions in the derivation presented here (including the power of c assigned to M ). Note that in the discussion in this subsection, we do not comment on the scales and times in the universe where the different branch conditions would be realised numerically. Rather we are just showing what the branch conditions are, what the equations are in the limiting case of each branch, and that the branch conditions and equations are consistent across the different cosmological regimes (so e.g. the condition to be in the deep MOND regime is not different depending on whether the cosmological density contrast is non-linear or perturbative). We leave a full numerical solution of this model to future work, here we just note that the condition for the different branches is something that needs to be considered within a cosmological Boltzmann code if one is seeking to calculate predictions for a theory with MOND on cosmological scales.

The G 0i equation
If we apply α = 0 and the same ansatz for F (K) (in the deep MOND limit 0 < K K * ) to the scalar part of the G 0i equation, this equation behaves similarly to the Poisson equation: the leading order GR terms are again cancelled by the leading order parts of F K . The scalar part of G 0i is given by where any terms that are not pure scalars are considered to have had their pure vector (divergenceless) parts subtracted. The leading order F (K) term is Since φ P = ψ P at leading order for any F (K), this means all of the non-matter terms vanish. In other words, the MONDian function that is designed to cancel the leading order non-matter part in the GR Poisson equation, also cancels the leading order non-matter part in the GR G 0i equation.
As for the Poisson equation we can go beyond leading order. The subleading terms that create MONDian behaviour in the Poisson equation will also be the terms of lowest remaining order here (these terms are of order c −4 , whereas the large scale B2 terms in equation (C.6) are of order c −5 ), so the scalar G 0i equation in the deep MOND limit is given by To the authors' knowledge, such an examination of the G 0i Einstein equation in the deep MOND regime has not previously been carried out. It is interesting that a similar cancellation occurs in both the Poisson equation and the G 0i equation for the same choice of F (K) and α = 0. We leave to future work an investigation of this potential coincidence, its possible physical significance, and whether it occurs in other relativistic extensions of MOND. As for the G 00 equation, this case can be included in the all-scales G 0i equation by simply returning a few of the terms that were dropped earlier, where again the explicit form of the large scale term LS B2 is given in appendix C. This equation replaces equation (3.22) when it is desired to include GEA theories with cosmological MONDian behaviour.

N-body simulations
The all-scales equations derived above can be implemented in cosmological N-body simulations to give simulations that represent GEA cosmologies, with a consistent background expansion, and where the inhomogeneities are evolved correctly no matter whether they are outside the horizon, around the horizon scale, well inside the horizon, or in the non-linear regime. In practice, for most cosmological simulations, particularly those with a smaller box size, the "large scale" terms can be neglected in the equations and this should make little difference to the output of the simulations. Initial conditions for these simulations can be -20 -

JCAP06(2023)006
calculated from a Boltzmann solver that implements GEA perturbation theory. Due to the smooth connection during the matter dominated era, these initial conditions and simulations will form a self-consistent and coherent complete cosmological picture for GEA theories. In particular, if models with α = 0 (and F (K) as described above for MONDian behaviour) were run, these simulations would correspond to the first fully consistent N-body simulations for a relativistic extension of MOND where MONDian behaviour arises cosmologically. We can also use these equations to examine how GEA cosmologies relate to cosmological MOND N-body simulations that have already been run [30][31][32][33][34][35][36]. These works themselves acknowledge the complications of running consistent cosmological MOND simulations, and in particular include discussions about different choices and options regarding the background expansion history and the initial conditions. The point we wish to focus on here is a further issue, which was alluded to in [36]: what are the actual equations governing the inhomogeneities in a consistent cosmology for a relativistic extension of MOND, and does this correspond to an implementation of MOND as carried out in N-body simulations that have been run so far.
As discussed earlier, when MOND is phrased as a modification to the Poisson equation, it is perhaps intuitive when one is used to GR to think that this can be applied to the Poisson equations that arise in different regimes, however these are not necessarily so similar in a modified gravity theory, and particularly a relativistic extension of MOND. This is borne out explicitly by the equations that we have derived above: for α = 0, although these models may still contain MOND on galactic scales, a cosmological N-body simulation should not contain MONDian behaviour, even if the free function F (K) in these theories is the same cosmologically as on galactic scales. So, despite having MOND on galactic scales these models do not relate to cosmological MOND N-body simulations.
For models with α = 0 however, there is a connection between MONDian behaviour on galactic and cosmological scales. In this case, we have shown that a Poisson equation with both MONDian and Newtonian branches exists, that the condition to be in either of these branches matches the usual MOND condition, and that both branches can be realised in principle. In particular, the Newtonian branch Poisson equation matches that in GR (for the vast majority of F (K), i.e. unless F (K) in the Newtonian branch is chosen very specifically to avoid this). In addition, we have shown that in these models the background Friedmann equations are not modified from those in GR. Combined, these results mean that existing MOND N-body codes can be used to represent structure formation in these theories, for a particular choice of matter, as long as the simulation parameters for the matter (for the background and the particle content) are chosen consistently, as they would be in GR. It may be however that the initial conditions used in previous works do not match those required to interpret the output of these simulations as part of a consistent GEA cosmology; determining this is beyond the scope of the work here. As noted above, the equations in this paper can be implemented in a Boltzmann code and used to generate initial conditions that would be consistent with interpreting the output of the N-body codes as representative of structure formation in these GEA models containing cosmological MOND.

Implications for GEA theories' ability to replace cold dark matter
It appears that GEA theories that have cosmological MONDian behaviour are quite limited, in that the Newtonian branch of the quasi-static (linear and non-linear) Poisson equation and the Friedmann equations are all essentially restricted to have their GR form. Given the range of evidence for cold dark matter and the different scales and environments in which cold dark matter phenomenology appears, these theories are probably too restricted to make -21 -

JCAP06(2023)006
for good alternatives to ΛCDM without themselves adding additional matter species to the universe, although we do not investigate this issue in detail here.
However, we also see that it is not required for GEA theories to have cosmological MOND just because it arises on galactic scales. It may be that the GEA theories with only galactic MOND can keep the successes of the MOND paradigm without some of the problems (such as galaxy clusters), but determining this will require running cosmological simulations with the equations derived here and no cosmological MOND behaviour. The same situation may be true for other relativistic extensions of MOND. It may be that attempts to build a MOND+sterile neutrino paradigm (see e.g. [33,63]) can benefit from the finding that MOND on cosmological scales is not required by having MOND on galactic scales; we leave a detailed examination of this prospect to future work.

Conclusion
To fairly and fully compare MONDian cosmologies against ΛCDM, we need to be able to examine cosmological structure formation on all scales in relativistic extensions of MOND, in a coherent and consistent fashion. In this paper we have laid out a method for deriving equations that govern the evolution of cosmological inhomogeneities on all scales in these theories. For any given theory, these equations allow consistent cosmological simulations to be run, and to examine when and how MONDian behaviour arises cosmologically. We have illustrated this method using the concrete example of GEA theories, resulting in the set of equations (3.21)-(3.26) that govern the cosmological background and evolution of cosmological inhomogeneities on all scales in a unified manner with a single F (K).
These equations show that MONDian behaviour does not necessarily arise cosmologically, even if the free function F (K) in these theories is the same cosmologically as on galactic scales. Specifically, cosmological MONDian behaviour requires the model parameter α to be zero, and thus the background expansion history of the universe and the Newtonian branch of the cosmological Poisson equation should both be governed by the same equations as in GR. Although structure formation as carried out in existing MOND N-body simulations can correspond to these GEA theories, in practice these theories are quite limited, but it seems difficult to get cosmological MONDian behaviour without these limitations. It may be that the GEA theories with only galactic MOND can keep the successes of MOND without some of the problems (such as galaxy clusters), but determining this will require running cosmological simulations with the equations derived here. These results demonstrate the strengths of our method for deriving a single set of equations that govern cosmological structure formation on all scales.
We have also laid out the assumptions and checks required to validate our method for a given cosmology and commented on these for our case study; notably we derived the source term for the vector potential in GEA theories, equation (3.26), which can be extracted from N-body simulations and should be suitably small on all scales. If larger than in GR (but small enough to not jeopardise the 1 c expansion), this vector potential can also generate novel lensing phenomenology and thus provide a smoking gun for modified gravity behaviour on non-linear scales [64].
The process laid out in this manuscript can be applied to any relativistic extension of MOND, thus paving the way for consistent cosmological simulations of these theories covering many decades in scale in order to fairly and fully compare them to the ΛCDM paradigm. It may also be the case that simply by deriving and analysing the equations for other relativistic extensions, as we have done here, the connection between some of these theories and galactic -22 -

JCAP06(2023)006
MONDian behaviour becomes clearer, potentially highlighting shortcomings in these theories without the need to run simulations. However, we caution against inferring what conclusions might be reached by such studies of a specific relativistic extension of MOND from the single study carried out here. It is clear from this work that the cosmological phenomenology of such theories is very rich, and trying to understand each of these theories consistently across all scales, including in which regimes MONDian behaviour can manifest, is not straightforward.
In particular, one recent theory that it would be interesting to apply our method to is the AeST (Aether Scalar Tensor) theory [28,29]. This theory satisfies the generic constraint from the CMB that the additional fields in these theories must behave similarly to the behaviour of cold dark matter in linear perturbation theory in the conditions that exist in the early universe [65,66], but the extra fields in this theory can still behave substantially differently to cold dark matter in other regimes, potentially giving rise to MONDian behaviour (on galactic or cosmological scales) or cosmological structure formation that proceeds differently to GR+CDM structure formation. In this sense, whilst the theory can be argued to possess dark matter for the purposes of CMB calculations, it doesn't possess particle dark matter in the sense of the "CDM" in ΛCDM, or in the sense of a dark matter particle of the kind that would be detected in direct, indirect, or collider searches for dark matter. The method laid out in this paper provides a way to examine cosmological structure formation in this theory and determine in which regimes the behaviour is similar or different to GR with cold dark matter. Interestingly, in AeST theory the role of the cosmological dark matter overdensity involves the time derivative of a scalar field perturbation, which may require careful thought when defining a Newtonian limit.

A Linear perturbation equations
In this appendix we present the linear perturbation equations written in terms of the resummed potentials. Note that we have split things into their background and perturbative parts, i.e. F K →F K + δF K and M 2 F → M 2F + M 2 δF .
Space-Space components of the Einstein equations.
Time-Space components of the Einstein equations.

B Matter equations
These are the "all scales" matter equations from [41] dv i dt = −Hv i + φ P,i a (B.1)