Aether Scalar Tensor (AeST) theory: Quasistatic spherical solutions and their phenomenology

There have been many efforts in the last three decades to embed the empirical MOND program into a robust theoretical framework. While many such theories can explain the profile of galactic rotation curves, they usually cannot explain the evolution 15 the primordial fluctuations and the formation of large-scale-structures in the Universe. The Aether Scalar Tensor (AeST) theory seems to have overcome this difficulty, thereby providing the first compelling example of an extension of general relativity able to successfully challenge the particle dark matter hypothesis. Here we study the phenomenology of this theory in the quasistatic weak-field regime and specifically for the idealised case of spherical isolated sources. We find the existence of three distinct gravitational regimes, that is, Newtonian, MOND and a third regime characterised by the presence of oscillations in the gravitational potential which do not exist in the traditional MOND paradigm. We identify the transition scales between these three regimes and discuss their dependence on the boundary conditions and other parameters in the theory. Aided by analytical and numerical solutions, we explore the dependence of these solutions on the theory parameters. Our results could help in searching for interesting observable phenomena at low redshift pertaining to galaxy dynamics as well as lensing observations, however, this may warrant proper N-body simulations that go beyond the idealised case of spherical isolated sources.


INTRODUCTION
The presence of an invisible matter -dark matter -permeating the Universe throughout its cosmic evolution is the most popular explanation to the origin and number of large-scalestructures in the Universe as well as their internal dynamics.Yet dark matter particles remain elusive despite the vast array of experimental programmes searching for such particles with great sensitivity, see e.g.Agrawal et al. (2021); Aalbers et al. (2022); Aprile et al. (2023).
While the particle dark matter hypothesis is still far from being excluded, the task of coming up with alternative theories which fit galactic and cosmological observations remains important and becomes increasingly so, especially in the absence of direct evidence for particle dark matter.One of such alternatives, namely Modified Newtonian Dynamics (MOND), was formulated in Milgrom (1983) and Bekenstein & Milgrom (1984) to explain the anomalous rotation curves of galaxies at large radii.MOND postulates that Newton's second law must be modified for accelerations a with magnitude smaller than a0 ∼ 1.2 × 10 −10 m/s 2 to explain the rotation curves of galaxies (Begeman et al. 1991;Milgrom 1988a).This is neatly captured by changing Newton's sec-E-mail: peter.verwayen@sydney.edu.auond law to µ(y) a = −∇Φ, with y ≡ | a|/a0 and Φ being the gravitational potential determined from the standard Poisson equation sourced by baryons only.This works if µ(y) y for y < 1 (this is the so-called MOND regime) and µ(y) = 1 for y > 1 (this is the Newtonian regime).
MOND faces a fundamental theoretical problem: it is not a relativistic theory but rather an empirical model.As such, MOND cannot be used to compute the formation and distribution of large-scale-structures in the Universe.A number of relativistic theories, typically extensions of General Relativity (GR), have been proposed to palliate this issue, and all have been constructed to lead to MOND behaviour in the non-relativistic limit (see Famaey & McGaugh (2012b) for a review).
One of these theories, the Tensor-Vector-Scalar (TeVeS) theory (Sanders 1997;Bekenstein 2004), introduces a unittimelike vector and a scalar field in addition to the usual gravitational metric tensor.These fields are combined together to define a second metric tensor which is used to determine the geodesics of ordinary standard model matter.With this bimetric disformal (owing to the vector field) structure, TeVeS leads to the equality of dynamical and lensing mass and thus can generate the right amount of gravitational lensing produced by baryon-only galaxies as if dark matter was present.Perturbations of the vector field also play a major role cosmologically as their time-evolution can lead to matter power spectra in line with observations (Skordis et al. 2006;Dodelson & Liguori 2006).However, despite these important improvements over other theories, TeVeS fails to fit the angular distribution of the CMB (Skordis et al. 2006).Moreover, TeVeS leads to tensor mode gravitational wave speed cGW which is different than the speed of light cEM and thus has been ruled out by the simultaneous observation of the GW170817 and GRB170817A events (Abbott et al. 2017).
The recent Aether Scalar Tensor (AeST) proposal of Skordis & Zlosnik (2021) retains the unit-timelike vector field and the scalar from TeVeS, but only has one metric tensor and thus, no bi-metric structure.Unit-timelike vector fields have been dubbed "aether fields" in other instances, see Jacobson & Mattingly (2001), hence, the naming of this theory1 .The success of AeST theory rests on maintaining cGW = cEM in all situations and fitting the CMB and matter power spectrum data quite convincingly, while retaining a MOND limit in galaxies and the correct gravitational lensing.
Cosmologically, the scalar in AeST theory evolves as in shift-symmetric k-essense (Scherrer 2004) which results in its cosmological energy density being similar to dust, i.e. ∝ (1+z) 3 plus small decaying corrections.This k-essence-like behaviour of AeST theory leads to spontaneous breaking of time diffeomorphisms as in the Ghost condensate (GC) theory (Arkani-Hamed et al. 2004a) which results in the metric potential Ψ (see Skordis & Zlosnik (2021, 2022) for details) acquiring a mass term µ.The quasi-static long-distance behaviour of AeST theory thus departs from MOND, but is still different than the GC due to the presence of the noncanonical MOND term and the vector field.Our aim in this paper is to precisely understand the particular behaviour of AeST theory in spherically symmetric static situations, due to the new features introduced by µ.
We present the AeST theory in Section 2, and derive the gravitational equations for quasistatic distributions of matter.In Section 3, we study the properties of density profiles assuming spherical symmetry.This work is critical in view of creating N-body simulations of large-scale-structure formation in the late Universe to compare AeST theory with non-linear observables, such as, those related to galaxy clusters (including the bullet cluster).It can also provide a guide on how extensions of GR at both early and late times may be distinguied from models based on the dark matter hypothesis.

THEORETICAL MODEL
We now summarise the main features of AeST theory as proposed in Skordis & Zlosnik (2021), the theory parameters, and the transition scales between the three regimes (Newton, MOND and "µ-domination") which feature in this theory.More details regarding the underlying setup and derivation of the weak-field static equations may be found in Skordis & Zlosnik (2021, 2022).

Static weak-field equations
The fields of AeST theory are the metric gµν , a unit-timelike vector field Aµ, such that A µ Aµ = −1, and a scalar field φ.Perturbing the field equations around Minkowski spacetime and ignoring time derivatives reduces the field dependence to just two potentials: Φ and χ, see Skordis & Zlosnik (2021, 2022).The gravitational potential Φ results from the metric perturbation.Matter fields couple minimally to this potential and follow its geodesics such that the Einstein equivalence principle holds.Specifically, gravitational accelerations are determined from ∇Φ.The potential χ is a gauge-invariant combination of a scalar mode contained2 in the spatial part of the vector field Aµ, and the perturbation of the scalar field φ.In addition to Φ and χ, we find it useful to define the potential Φ through treating Φ and χ as the fundamental fields and Φ as a derived variable.With these variables, the weak-field equations of AeST thaory which extend the standard Poisson equation for Φ and determine spatial dependence of the gravitational potentials of the theory Φ, χ and Φ, can be written in the form while a third, non-independent equation, is In the equations above, GN is Newton's constant, ρ b is the baryonic density, J = J (Y) is a function of the variable Y = (g µν + A µ A ν )∇µφ∇ν φ, µ is the mass parameter and β0 is the inverse screening parameter, the screening parameter being λs ≡ 1/β0.The latter determines how much χ contributes to Φ through (1) in the large gradient limit of (2), ( 3) and (4); see below.The mass parameter µ is determined through a specific combination of parameters which enter the action of the theory, see Skordis & Zlosnik (2021, 2022) for details, however, for our purposes we treat it as a free parameter in this work.It can be shown that linear stability of the theory on Minkowski space imposes λs > 0 and µ 2 > 0, however, in practice µ −1 1 Mpc to ensure that MOND solutions are attainable at galactic scales.Indeed, precise observations of the extent of flat galactic rotation curves can be used to put strong constraints on µ −1 (Mistele et al. 2023).The variable Y denotes the "spatial gradient" combination Y = (g µν + A µ A ν )∇µφ∇ν φ and in the static weak-field case one finds that (5) The function J then controls the behaviour of static weak field configurations through its dependence on Y, or more specifically on √ Y = | ∇χ|.When √ Y a0 the function J is constructed so as to recover GR as close as possible.This can happen if which also serves as the definition of the screening parameter λs, a parameter which must be part of the function J .As λs → ∞ (equivalently β0 → 0), then χ → 0 so that the contribution of χ to Φ through (1) is totally screened3 .In the opposite limit, √ Y a0, the function J is constructed to lead to MOND behaviour, which is possible if Three comments are in order: • In deriving (2)-( 4) one also finds that the two standard weak-field metric potentials Φ and Ψ are equal.This ensures that the lensing mass is equal to the dynamical mass.This means that for any source configuration where the MOND limit of AeST theory successfully replaces dark matter, the gravitational lensing signal for the same source will also be as if dark matter is present.
• The case χ = 0 corresponds to having Φ = Φ and the above system of equations reduces to a single equation for Φ which is identical to the non-relativistic weak-field limit of the Ghost condensate model (Arkani-Hamed et al. 2004b).
• Because of the term µ 2 Φ, the system of equations ( 2)-( 4) is not precisely that of TeVeS theory but tends to it as µ → 0. However, dust-like solutions in cosmology require µ to be non-zero, with µ → ∞ corresponding to the pure Higgs phase (Arkani-Hamed et al. 2004b) where the solutions are exactly dust-like throughout the entire history of the Universe.Thus AeST theory will have solutions which depart from the MOND paradigm at a large enough distance away from the matter source.This results in a third regime, different from both GR and MOND as we study thoroughly below.

Newtonian, MOND and µ-domination limits and corresponding transition scales
We define the dimensionless variable and the new function which corresponds to the usual MOND interpolation function and is designed to provide a continuous transition between the Newtonian and MOND regime as follows: f → λs : large gradient limit, x 1 x 1+β 0 : small gradient limit, x 1 .
(10) with β0 ≡ 1/λs (see also the equivalent limits in terms of J in equations ( 6) and ( 7)).We note that the large gradient limit corresponds to Newtonian behaviour to lowest order in the potentials and GRlike behaviour in the full relativistic strong-field case (up to parametrically small post-Newtonian corrections).
2.2.1 Large gradient limit: strong-field/Newtonian regime Consider first the large gradient limit which leads to the strong-field (and Newtonian) regime.Setting f = λs in equation (3) and integrating twice, results in χ = β0 Φ + χ0, where we have neglected a possible curl which is not relevant in spherically symmetric situations (see also the discussion in 3.3).Substituting this relation in (1) gives Φ = (Φ − χ0)/(1 + β0) and with these conditions (2) takes the form of the inhomogeneous Helmholtz equation Equation ( 11) also results from taking the non-relativistic weak-field limit of the Ghost condensate model (Arkani-Hamed et al. 2004b).For scales smaller than ∼ µ −1 / √ 1 + β0, (11) reduces to the standard Poisson equation describing Newtonian gravity.The solution to (11) for scales larger than ∼ µ −1 / √ 1 + β0 is not physically relevant to the models studied in this article, since the systems we study first transition to the small gradient regime before the scale ∼ µ −1 / √ 1 + β0 is reached.
Note that the integration constant χ0 does not affect (11), and so, does not play a role in the large gradient regime.However, it's value will survive until the small gradient regime, and as we discuss below, it will lead to observable effects at very large distances from the source.

Small gradient limit: MOND and µ regimes
Consider now the small gradient limit where f = x/(1 + β0); see (10).Then (4) becomes which reduces to leading to the MOND equation for χ in the limit where µ → 0 and Eq. ( 2) turns into The above equation has an exterior solution scaling as | ∇ Φ| ∼ 1/r 2 while the deep MOND solution to (13) (for µ = 0) is | ∇χ| ∼ 1/r.Hence, for large enough r (but not too large so that we keep within the regime where µ has no influence and can be set to zero) we have | ∇χ| ∼ 1/r implying that the solution for Φ is also MONDlike.
When r becomes even larger, the (1 + β0)µ 2 Φ term in the above equation becomes important, leading to a new regime specific to AeST theory which departs from MOND.This happens because since | ∇χ| ∼ 1/r, the 1st term is expected to scale as 1/r 3 , while the 2nd term grows as ln r and is always bound to become important at some scale rC.We call this the µ-dominated regime.
We now derive the transition scales from Newton to MOND and from MOND to µ-domination in the case of a point source.

Transition from Newtonian to MOND regime
The transition from Newtonian to MOND behaviour in the "classical" MOND paradigm occurs when the Newtonian force is equal to the MOND force, and is given by the familiar MOND radius equation However, in AeST theory, the MOND force is coupled only to the χ component of the field.Therefore, determining the radius at which the total force deviates from the Newtonian force is determined by the χ component only.In the Newtonian limit of AeST theory, χ = β0Φ/(1+β0) with ∇Φ tracking the Newtonian force (i.e.| ∇Φ| = GN M/r 2 ).Hence ∇χ departs from the Newtonian limit and enters the MOND regime when ∇χ is no longer is proportional to | ∇Φ| = GN M/r 2 .This happens when the two limiting situations in (10) become comparable, that is, when λs ∼ x/(1 + β0).Hence, different transition scales may occur depending on whether we consider interior or exterior solutions.In the exterior case, we call this transition rχ where such that rχ rM always.Given Eq.( 16) and the fact that β0 1, we expect rχ rM .To determine when the total AeST force enters the MOND regime, we also need to require that χ is the dominant component of Φ, that is ∇Φ ∇χ.This domination occurs at a distance rM which is defined as the point when such that rM rM always.In practice, β0 is expected to be small so that the two scales rM and rM are approximately equal.
The different transition scales above are depicted in Fig. 7 where we discuss how the solutions change when we vary the parameter λs.

Transition from MOND to µ-domination
To estimate the transition between MOND and µ-domination we consider (13) for a point source of mass M which we define to occur at the scale rC .Estimating derivatives as d/dr ∼ 1/r and considering the Log terms appearing in the MOND solution Φ, i.e. ln(r/rM), as being O(1)4 , the terms on the LHS in ( 13) become comparable at the scale rC when where χ ∼ Φ ∼ √ GNM a0×O(1).This gives a simple estimate of rC ∼ rM/µ 2 1/3 .The above estimate, however, ignores possible effects coming from the boundary condition and these can be important.Let us set the boundary condition for χ at a radius r0, that is, χ(r0).For later convenience, we normalize it with respect to √ GNM a0 and define χout ≡ χ(r0)/ √ GNM a0 as the constant free parameter setting the boundary condition.Consider now the extreme case where χout is so large that it dominates the solution, while ∇χ still retains a (subdominant) 1/r MOND component.Then, Φ ∼ χ ∼ √ GNM a0 χout and ( 18) leads to rC ∼ rM/µ 2 χout 1/3 .The last relation above would appear to diverge when χout → 0, but nothing obviously wrong should happen when setting the boundary condition for χ to zero (and this is verified numerically).Moreover, the point r0 where the boundary condition is set is arbitrary, thus the "zero point" where χ(r0) (or equivalently χout ) is zero is also arbitrary.
Nevertheless, defining χout can be done in an unambiguous way, so that a better estimate for rC which interpolates between the two extreme cases above can be derived.We summarise this here and leave the details for appendix-A.In the small gradient regime, see 2.2.2, the details of the interpolation function are unimportant and we may set f → x/(1+β0) from (10) resulting in (13).First, assuming that µ = 0 the solution is χ = √ GNM a0 χout + ln r rM , which serves to define χout as χout ≡ χ(rM)/ √ GNM a0.To determine when the full µ = 0 solution deviates from the pure MOND solution, we expand the former in terms of the MOND solution plus a small perturbation.We find that when χout takes a specific value that we denote as χ(max) out , the transition scale rC reaches a maximum.In appendix A we show that this specific value is given by We then denote deviations from this extreme value as ∆, so that a general value for χout is determined from thus, the actual value for the boundary condition is fully determined from ∆.
We then determine (appendix A) a better estimate for rC which is where the factor of 1/3 is inserted to create a more conservative estimate.Thus, our first naive estimate above (rC ), simply corresponds to the maximum rC case which is obtained when ∆ = 0, while the 2nd extreme case above is equivalent to ∆ 1.To summarize, rC denotes the scale where deviations from pure MOND solutions can occur due to the onset of the µ regime.As we show below, for distances r > rC the solution to the weak-field equations becomes oscillatory but still decaying (similarly to the Helmholtz equation).Importantly, one may use data to constrain rC, leading to global constraints on µ and specific constraints for ∆ which may be different for individual astrophysical sources.We discuss the dependence of the solution on ∆ in section 4.2.3.

NUMERICAL SOLUTIONS: SETTING UP THE PROBLEM
In order to investigate how AeST theory solutions differ from those of the classical MOND theory, we define a baryonic system representative of an idealized spherical galaxy and calculate the solutions of the field equations, as described in the previous section.In the case of MOND, solutions for the force for spherical systems can be obtained analytically starting from a Newtonian solution.However, due to the mass term µ, analytical solutions are no longer possible in AeST theory and we need to rely on a numerical approach.We describe in this section the setup that we use to obtain the spherically symmetric solutions.

Density profiles
We use two different density profiles for the analysis: a tophat and a Hernquist profile, the latter being a good description of the baryonic component of spherical galaxies, see Hernquist (1990).
The top-hat profile is defined as a sphere of radius a h of constant density embedded in a uniform background: where ρc is the baryonic density at the centre of the profile.For our purposes, we fix the parameters of the profile to ρc = 3.45 × 10 9 M kpc −3 which is ∼ 10 7 larger than the critical density of the Universe, and a h = 2 kpc.This is representative of a galactic system.We discuss how changing the density affects the solutions in Section 4.2.5.The Hernquist density profile is given by (Hernquist 1990), Note that for simplicity, we use the same symbol a h to represent the radius of the top-hat sphere as well as the scale radius for the Hernquist profile.The Newtonian potential associated to this profile is given by, Binney & Tremaine (2008) ΦN 10 2 10 1 10 0 10 1 10 2 10 3  27), ( 28) and ( 29).For each function, we show the result for three different values of λs: 1, 2 and 10.

AeST field equations for spherical isolated sources
In this section, we focus on solving (2) and ( 4) in order to predict the radial dependence of Φ and χ which are critical for understanding the behaviour of gravity for spherical isolated sources in AeST theory.The potential Φ used for particle accelerations is calculated from (1).The spherical version of (2) and ( 4) from which we obtained the numerical solutions is where x = |∇χ|/a0.For setting boundary conditions we employ the analytic solutions of these equations at the (non-zero) radius where we start the numerical integration.Throughout this work, we use a0 = 1.2 × 10 −10 m s −2 (Begeman et al. 1991).We will use different expressions for the function f .Two are given in the literature (Famaey & Binney 2005;Zhao & Famaey 2006): which agree with the limits defined in (10).However, we will also explore the consequences of a different function that we define as which exhibits a sharp Newtonian to MOND transition and returns the correct limits defined in (10).This function has a turning point at x → λs + 1 and is not designed to be fully consistent with observations, but rather as a test function to show the Newtonian to MOND transition radii and viceversa.Fig. 1 shows these three functions for three different values of λs.See McGaugh (2008) for additional discussion on interpolation functions.

Analytic solutions for the top-hat profile
Analytic solutions are required for setting up inner boundary conditions for the numerical solvers as well as for testing that our numerical implementation of the complete solutions is correct, at least for the cases where µ = 0 or r rC , for which these analytic solutions exist.While our numerical solutions are obtained using ( 25) and ( 26), the analytical calculations are easier to obtain using (25) together with the spherical version of (3).
We now derive the analytic solutions for Φ and χ in the case of a top-hat source profile with the simple interpolation function.In appendix B we present solutions for the tophat profile with the sharp interpolation function, while in appendix C we display the solution for the Hernquist profile with the simple interpolation function.
The solution for Φ can be obtained by re-scaling the gravitational constant in the solution of the standard Poisson's equation (e.g.Binney & Tremaine 2008) and takes the following form: We then integrate (3) once to obtain The field ∇×k is divergenceless and was discussed in detail in Bekenstein & Milgrom (1984) where it has been shown to be exactly zero for particular symmetries (including spherical), and to behave at least as r −3 for non-symmetric situations (the effects of this so called curl term on non-linear structure formation with pure MOND were studied in detail by Llinares et al. (2008) and Llinares (2011)).Since we are assuming spherical symmetry, we can ignore ∇×k and invert equation ( 32) to find ∇χ.Applying this procedure with the simple interpolation function ( 27) gives which can be integrated to obtain a solution for χ.Letting M ≡ 4 3 πa 3 h ρc be the total mass of the system, and defining the scale rI ≡ 4a 3 h /r 2 χ , the solution is where χin = χ(0)/ √ GNM a0 and χout are two integration constants corresponding to the inner and outer regions respectively, both normalized to √ GNM a0 for later convenience.These two integration constants are not independent but are related by matching the solution at the boundary r = a h in (34).Lastly, we start our numerical integration, at radius r0 (in the interior) and set the boundary condition as χ0 = χ(r0)/ √ GNM a0, which by appropriate use of (34) may then be related to χin and χout.We study the dependence of the solution on the boundary condition in Section 4.2.3.

RESULTS
We first analyse the dependence of the numerical solutions on the model and density profile parameters.We then present a list of physical effects which make AeST theory different from standard MOND.

General properties of the solutions
Fig. 2 shows the numerical solution for the potentials (left) and their derivatives (right) for a fiducial set of parameters (λs, µ) = (1, 1 Mpc −1 ), the "simple" interpolation function f (x) from ( 27) and the Hernquist density profile defined in (23).The blue solid curve in Fig. 2 corresponds to the Newtonian solution for Φ which we denote as ΦN and the green solid curve is solution for Φ in AeST theory with µ = 0 which is equivalent to MOND, thus denoted as ΦMOND.Both ∇ΦN and ∇ΦMOND agree in the central region of the galaxy and thus give rise to the same force profile there.Farther away from the center, the gradients fall below a0 (see horizontal lines in the right panel) leading to the characteristic logarithmic MOND potential for ΦMOND outside the source and a force that follows a 1/r relation.
The three additional curves in Fig. 2 in both panels are solutions of the field equations ( 25) and ( 26) and the total potential provided by the relation (1).The green dotted curve is the solution for Φ (and ∇ Φ); the yellow dashed curves are the solution for χ (and ∇χ); and finally, the pink dash-dot curve is the solutions for the potential Φ (and ∇Φ).As we move farther away from the center, between the MOND radius rM and rC, the theory tends to the classical MOND behaviour (i.e. a force law which not only has the same dependence with r, but also the same normalization).Farther away when r > rC , the solutions enter the oscillatory regime, where the potential develops additional potential wells and the force can become repulsive.
In the next two sections, we describe how variations of the AeST theory parameters affect these reference solutions and what additional physical effects are associated with them.).The blue, green and red regions delineate the Newtonian, MOND and Oscillatory regions respectively.The yellow and green dashed lines are for the potentials Φ and χ respectively, and the purple dash-dot line is the potential Φ which is responsible for defining particle accelerations through its derivative.We have included the Newtonian (blue) and classical MOND (green) solutions for comparison.The break in the blue curve at ∇Φ = 10 −5 is not physical, but related to the symlog scaling that we use for the vertical axis of the right panel. .

Dependence on the free parameters
We now study the dependence of the solutions on the free parameters of the model and on the central density (and mass) of the source, assuming a top-hat density profile (22).

Interpolation function
Fig. 3 shows solutions for the gravitational potential (left column) and force (right column) that were obtained with the three interpolation functions f (x) defined by ( 27), ( 28) and ( 29).The first three rows correspond to a different potential; from top to bottom, these are Φ, χ and Φ.We fix the boundary condition χ(r0) at r0 to be the same for all three interpolation functions, hence, the inner Newtonian regime is identical in all cases displayed in the first three rows of Fig. 3.For the same reason, the MOND regime for the forces (right panels of Fig. 3) is also identical in all solutions, however, the potentials (left panels of Fig. 3) differ.The reason for the potentials reaching different values as the system evolves towards the MOND regime, is the sharpness or smoothness of the interpolation function f (x).Smooth transitions from the Newtonian to the MOND regime lead to higher values of the potential before it settles into the MOND track (see also Fig. 1), after which, its evolution becomes the same for all interpolation functions.When r rC, the solutions enter the oscillatory regime, and the value of Φ = Φ + χ then plays a role due to the mass term µ (see ( 25) and ( 26)).At that stage, since the value of the potential between interpolation functions is already different, the result is a change in the oscillation pattern, although the overall oscillation envelope stays the roughly the same; see Fig. 3.
The last row of Fig. 3 shows the result when the boundary conditions are calculated individually for each interpolation function according to the prescription described in section 4.2.3, as well as in Appendix A and B. With the boundary conditions appropriately calculated, the µ-dominated region shows identical behaviour across all the interpolation functions shown.This is expected, as all interpolation functions have the same limiting form at small gradients.However, given that they now have different boundary conditions at r0, the inner Newtonian (and partially the transition to the MOND regime) is slightly different across each function.We do not show the exponential interpolation function since we do not have an analytical solution for it in order to specify the χ boundary condition according to our prescription of section 2.2.4 and Appendix A.
For the remainder of this work will mainly focus on the "Simple" interpolation function ( 27) so as to simplify the plots and standardise the results.

Mass parameter µ
Setting the mass parameter µ to zero gives a standard MOND-like behaviour in which the force follows the Newtonian solution in the inner regions at r < rM and transitions into a MOND regime when r > rM.Letting µ take positive values initiates spatial oscillations in the solutions for r rC .Choosing λs = 1 and top-hat profile (leading to rM = 12.5 kpc), Fig. 4 shows solutions with three different values of µ = {0, 1, 10} Mpc −1 , corresponding to rC: {∞, 156, 33.6} kpc; see (21).All solutions for the force give the same result up to ∼ rC.After this point, the solutions de- ).On the left we show the potentials and on the right their derivatives.From top to bottom, the panels contain Φ, χ and Φ.The horizontal dash-dot lines in the right panels denote Milgrom's constant a 0 .For each panel, the three lines each show an interpolation function and we also show the Newtonian solution in blue for comparison.The break in the blue curve at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the right panels.viate from each other with the oscillations in the µ-dominated region starting at smaller r for larger µ.Observations of galaxies do not show evidence for the presence of repulsive gravitational forces affecting their internal dynamics, but only evidence of the MOND regime.Requiring rC to be larger than the virial radius of the Milky Way (∼ 200 kpc) gives an estimate that µ −1 1 Mpc.Param-eter estimation using several different astrophysical objects and proper modelling of the observable quantities may provide more accurate bounds on this parameter but this goes beyond the scope of our current work.See however (Mistele et al. 2023) where the first attempt in doing so is considered.Fig. 5 shows Newtonian and AeST forces for the region r rC assuming µ = 1 Mpc −1 and λs = 1.Thanks to the non-linearity of the AeST theory equations, the oscillations are far from sinusoidal, acquiring rather a very steep slope when crossing zero.We generally find that the mean spatial frequency of the oscillations decreases with µ, however, given the non-linear nature of the equations, we do not have an analytical estimate for the exact relation.Furthermore, the wavelength is not constant through the domain, but decreases slightly with r.Further detail of the zero crossing is shown in the bottom panel of Fig. 5 indicating that the 2nd derivative might become singular at the crossing.This, however, is an artifact, and an alternative way of solving the system of equations (2) and (4) using a Hamiltonian approach (Durakovic & Skordis 2023) shows that nothing bad happens there.

Boundary conditions
The Poisson and classical MOND equations depend on the potential only through its derivatives, hence, adding a constant to the solution is not physically relevant and does not affect the gravitational force.The situation is different in AeST theory due to the presence of the mass term µ 2 Φ in (25) and ( 26) which makes the value of the potentials physically relevant at r > rC.In other words, the boundary conditions are important in AeST theory.We now investigate their impact on the solution determining the depth and radial profile of the gravitational potential Φ and its derivative.
As discussed in Sec.2.2.4 the boundary condition ∆ directly affects where the solution departs from pure MOND and en- .Sensitivity of the solutions to changes on the inner boundary condition employed to obtain numerical solutions.We show solutions for metric perturbation Φ for the top-hat profile and (λs, µ) = (1, 1Mpc −1 ).Left and right columns are the field and its radial derivative (i.e. the force that defines particle trajectories).The horizontal dash-dot lines in the right panels are the MOND constant a 0 .The two different rows correspond to mild (upper row) and large (bottom row) offsets on the fiducial boundary condition which we define as the value of the Newtonian solution in the center.The blue line is the Newtonian solution, shown for comparison.The break in the blue curve at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the right panels.
ters the µ-regime.This is estimated according to (21).The effect of varying ∆ is illustrated in Fig. 6 and as the numerics confirm, ∆ = 0 constitutes the fiducial zero point corresponding to maximum rC.
To investigate this, we choose a value χ0 at r = r0, such that, χout = χ(max) out + 1 − 1 2 ln 2 through the use of (34), where χ(max) out is defined in ( 19).This corresponds to ∆ = 0 in (20) as shown by the magenta dash-dot curve in Fig. 6.The expectation from ( 21) is that this also corresponds to the maximum rC when all the other parameters, e.g.µ and λs are kept fixed.Indeed, this is what is observed in the upper row of Fig. 6.Varying ∆ to negative (positive) values, results in χout being smaller (larger).However, in both cases rC, decreases from the value reached through ∆ = 0, as expected through (21) where ∆ enters with an absolute value.
There are two further related effects that are worth describing.Firstly, when ∆ is positive, the force ∇Φ passes through zero at smaller distances than the fiducial ∆ = 0 case.Hence, one observes a smaller oscillation phase.Conversely, for negative ∆ the opposite happens, that is, ∇Φ passes through zero at larger distances than the fiducial ∆ = 0 case.Hence, one observes a larger oscillation phase.This is seen in both upper and lower panels of Fig. 6.Secondly, when ∆ deviates significantly from zero, the force reaches an extensive plateau which is either negative (for ∆ > 0) or positive (for ∆ < 0) before entering the oscillatory regime.This effect is more pronounced in the lower panel of Fig. 6. .Sensitivity of the solutions to changes on the model parameter λs.We show solutions for metric perturbation Φ for the tophat profile and µ = 1Mpc −1 .In this particular figure we used the Sharp interpolation function to highlight the transition between the Newtonian and MOND limits.Top left and right plots are the field and its radial derivative (i.e. the force that defines particle trajectories).
The horizontal dash-dot lines are the MOND constant a 0 .The blue line is the Newtonian solution, shown for comparison.The break in the blue curve at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the top right panel.See section 4.2.4 for the explanation of the bottom panels.
In addition, λs controls the effective gravitational strength which couples Φ to the density ρ, see (25) and also (30) and (31).Fig. 7 shows the affect of varying λs on the solution of the potential Φ whilst keeping other parameters fixed.We have used the Sharp interpolation function (29) to show these effects clearly.The zoom-in portion of the figure (bottom left) shows the transition from Newtonian to MOND behaviour in the exterior of the source.We see the two scales rM and rχ at work, and their dependence on λs.The scale rχ signifies when the χ component, and the total field Φ = χ + Φ, changes from a Newtonian to MOND behaviour.The scale The transition to the full MOND force is essentially completed when the point rM is reached.We note that both Newtonian and MOND magnitudes of the force are independent of λsonly the transition is affected by this parameter.Towards larger radii, when r ∼ rC, the solution approaches the µ-dominated regime as discussed above.There is a mild dependence on λs as to where that happens, since rC ∝ λ 1/3 s /(1 + λs) 1/3 , and this is clearly seen in the zoom-in region (bottom right).The main effect of this shift, is the starting point of the oscillations, which sets the overall oscillation phase.

Central Density and total mass
In our discussion so far, we have assumed a fiducial central density representative of the bulge of the Milky Way or an average spherical galaxy (Widrow et al. 2008).In Fig. 8 we show the effect of changing this central density, on the total potential Φ (left panel), on the total force ∇Φ (middle panel) and on the total force when the mass is kept constant (by adjusting the source size a h ) (right panel) We observe that larger central densities (resulting in larger total mass), make the inner potentials deeper and increase the overall force throughout the Newtonian and MOND regimes, but also increase the magnitude of the oscillations in the µdominated regime.The opposite happens with decreasing the central density, and this overall behaviour is what we expect based on the transition scale dependence on the total mass M .Indeed, in the right panel, we see the effect of keeping M constant, which is that the exterior solution remains the same, and the only differences occur in the interior.

CONCLUSIONS
In this paper, we predict the profile of the gravitational potential and associated force expected in spherical galaxies in the AeST theory.The later is the first extension of GR that successfully fits the CMB angular and matter power spectra without a dark matter component.The non-relativistic limit of AeST theory differs from the classical MOND theory in the fact that the field equations (equivalent to the Poisson's equation in the standard gravity case) include a mass term.This new ingredient leads to a different gravitational phenomenology which we investigate here for two different spherical density distributions: a top-hat and a Hernquist profile.
We identified three characteristic regimes in the solutions, independent of the density profile: a Newtonian regime, a MOND regime and finally, an oscillatory regime, where the mass term dominates and the fields (as well as the forces and their associated dynamical mass) develop spatial oscillations.Focusing on the case where the MOND regime appears for intermediate scales and the oscillatory regime at larger scales, we find that the transition from the Newtonian to the MOND regime depends on the usual acceleration parameter a0, and the screening parameter λs, while the transition from MOND to the oscillatory regime depends on the mass parameter µ, the inner boundary condition and the total mass of the gravitating object through (21).This means that the oscillations do not appear at a fixed radius, but have a more complex dependence which is different for every object and depends on their mass distribution.On the other hand, we find no strong dependence of the solutions with the free function that regulates the speed of the transition between the Newtonian and MOND regimes.
There should be remarkable consequences from the distinct behaviour of the gravitational potential and forces observed here.Accurate predictions would require N-body simulations but we do expect the matter density and distribution to be impacted and different from that in the ΛCDM model.Given the oscillations at large scales corresponding to the µ dominated regime, one could speculate that this will translate into the existence of ring-like structures far away from the galactic centre, which may eventually resemble observed structures such as galactic rings (Buta & Combes 1996).Furthermore, the enhancement (suppression) of the AeST gravitational potential in the inner galactic regions for certain values of the inner boundary condition (see Fig. 6) and interpolation function (see Fig. 3) could be misinterpreted as a higher (lower) dark matter density in the ΛCDM framework.This could af-fect lensing analyses, galaxy cluster profiles, dark matter indirect (and potentially also direct) detection predictions, and presumably the whole galactic (and possibly early Universe) evolution.
The main conclusion of our work is that AeST theory is a potential alternative to particle dark matter but N-body simulations are needed to make sure that the non-linear formation of large-scale-structures in this theory and the gravitational structure of galaxies are consistent with observations.

Figure 1 .
Figure 1.Interpolation functions explored in this work.Different colours correspond to the different functions defined by (27), (28) and (29).For each function, we show the result for three different values of λs: 1, 2 and 10.

Figure 2 .
Figure2.Numerical solution for the potentials (left) and their derivatives (right) for the Hernquist density profile and the fiducial model parameters with (λs, µ) = (1, 1 Mpc −1 ).The blue, green and red regions delineate the Newtonian, MOND and Oscillatory regions respectively.The yellow and green dashed lines are for the potentials Φ and χ respectively, and the purple dash-dot line is the potential Φ which is responsible for defining particle accelerations through its derivative.We have included the Newtonian (blue) and classical MOND (green) solutions for comparison.The break in the blue curve at ∇Φ = 10 −5 is not physical, but related to the symlog scaling that we use for the vertical axis of the right panel.

Figure 3 .
Figure3.Sensitivity of the solutions to changes on the interpolation function f for the top-hat profile with (λs, µ) = (1, 1Mpc −1 ).On the left we show the potentials and on the right their derivatives.From top to bottom, the panels contain Φ, χ and Φ.The horizontal dash-dot lines in the right panels denote Milgrom's constant a 0 .For each panel, the three lines each show an interpolation function and we also show the Newtonian solution in blue for comparison.The break in the blue curve at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the right panels.

Figure 4 .Figure 5 .
Figure 4. Sensitivity of the solutions to changes on the mass parameter µ on the potential Φ (Left) and its derivative (Right) for the top-hat profile and λs = 1.The horizontal dash-dot lines in the right panels denote Milgrom's constant a 0 .The blue solid line is the Newtonian solution, shown for comparison.The other three lines correspond to three different values of µ = {0, 1, 10}Mpc −1 .The break in the curves at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the right panels.
Figure6.Sensitivity of the solutions to changes on the inner boundary condition employed to obtain numerical solutions.We show solutions for metric perturbation Φ for the top-hat profile and (λs, µ) = (1, 1Mpc −1 ).Left and right columns are the field and its radial derivative (i.e. the force that defines particle trajectories).The horizontal dash-dot lines in the right panels are the MOND constant a 0 .The two different rows correspond to mild (upper row) and large (bottom row) offsets on the fiducial boundary condition which we define as the value of the Newtonian solution in the center.The blue line is the Newtonian solution, shown for comparison.The break in the blue curve at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the right panels.
Figure7.Sensitivity of the solutions to changes on the model parameter λs.We show solutions for metric perturbation Φ for the tophat profile and µ = 1Mpc −1 .In this particular figure we used the Sharp interpolation function to highlight the transition between the Newtonian and MOND limits.Top left and right plots are the field and its radial derivative (i.e. the force that defines particle trajectories).The horizontal dash-dot lines are the MOND constant a 0 .The blue line is the Newtonian solution, shown for comparison.The break in the blue curve at ∇Φ = 10 −5 is related to the symlog scaling that we use for the vertical axis of the top right panel.See section 4.2.4 for the explanation of the bottom panels.

Figure 8 .
Figure8.Sensitivity of the solutions to changes on the normalization of the density profile.The left panel shows the total potential Φ and the middle and right panels shows its radial derivative.While in the left and middle panel the profile radius a h is kept constant (so that the total mass is different), in the right panel the mass is kept the same for all curves by adjusting a h .The continuous cyan curve is the fiducial model used in previous figures with ρc = 3.45 × 10 9 M kpc −3 .The other curves show perturbations around this value.The inner boundary condition for the potential is the same as used in previous figures.