Towards galaxy cluster models in Aether-Scalar-Tensor theory: isothermal spheres and curiosities

The Aether-Scalar-Tensor (AeST) theory is an extension of General Relativity (GR) which can support Modified Newtonian Dynamics (MOND) behaviour in its static weak-field limit, and cosmological evolution resembling ΛCDM. We consider static spherically symmetric weak-field solutions in this theory and show that the resulting equations can be reduced to a single equation for the gravitational potential. The reduced equation has apparent isolated singularities at the zeros of the derivative of the potential and we show how these are removed by evolving, instead, the canonical momentum of the corresponding Hamiltonian system that we find. We construct solutions in three cases: (i) in vacuum outside a bounded spherical object, (ii) within an extended prescribed source, and (iii) for an isothermal gas in hydrostatic equilibrium, serving as a simplified model for galaxy clusters. We show that the oscillatory regime that follows the Newtonian and MOND regimes, obtained in previous works in the vacuum case, also persists for isothermal spheres, and we show that the gas density profiles in AeST can become more compressed than their Newtonian or MOND counterparts. We construct the Radial Acceleration Relation (RAR) in AeST for isothermal spheres and find that it can display a peak, an enhancement with respect to the MOND RAR, at an acceleration range determined by the value of the AeST weak-field mass parameter, the mass of the system and the boundary value of the gravitational potential. For lower accelerations, the AeST RAR drops below the MOND expectation, as if there is a negative mass density. Similar observational features of the galaxy cluster RAR have been reported. This illustrates the potential of AeST to address the shortcomings of MOND in galaxy clusters, but a full quantitative comparison with observations will require going beyond the isothermal case.


Introduction
Galactic and extra-galactic phenomena cannot be accounted for by the dynamics of baryonic matter and radiation (including neutrinos) evolving under gravity described by general relativity (GR).Velocities of stars [1,2] and gas [3] in the outskirts of disk galaxies are found to asymptote, rather than decline.Assuming GR, or Newtonian gravity for weak fields, relaxed galaxy clusters need more than the baryonic mass to account for the velocity dispersion of galaxies [4] in them, for the thermodynamic state of the gas inferred through X-ray emission [5] and the thermal Sunyaev-Zel'dovich effect [6,7], and for the observed gravitational lensing [8].Cosmologically, the observed clustering of galaxies [9,10] and the anisotropies of the Cosmic Microwave Background Radiation are also in stark conflict with a model of the Universe based on GR containing only baryons and radiation [11][12][13].
The common approach to resolving all these discrepancies is to posit the existence of a dark matter (DM) component -a matter component with feeble or no interactions with baryonic matter other than through gravity.Within the DM paradigm, the most successful (and simplest) model is that of Cold Dark Matter (CDM), modelled as a collection of particles evolving with the Boltzmann equation.Cosmologically, this leads to a dust (pressureless matter) component which outweighs baryons approximately 5 : 1 and forms diverse halos and subhalos around galaxies and galaxy clusters.DM provides the additional gravitational force necessary in galaxies, and the additional mass to account for cluster observations.With the inclusion of a cosmological constant Λ, the resulting ΛCDM model is in broad agreement with all cosmological observations at the scale of ∼ M pc or larger.
Direct and indirect [14] searches for the dark matter particle have, however, been negative so far, leaving open the possibility of an extension of GR being responsible for the discrepancies observed.Moreover, there are regularities in the internal dynamics of galaxies permitting to describe them without a DM component but by assuming that the observed acceleration a obs is given by a obs = √ a 0 √ a N when gravitational accelerations are smaller than a threshold a 0 ∼ 1.2 × 10 −10 m/s 2 where a N is Newtonian acceleration sourced only by the observable baryons.This is the Modified Newtonian Dynamics (MOND) [15][16][17] proposal.Immediate consequences are: the independence of the rotational velocity on the orbital radius, and the baryonic Tully-Fisher relation which has broad observational support [18].
At accelerations larger than a 0 Newtonian gravity should be restored.This is generally achieved through an interpolation function M such that M (x) a obs = a N where x = a obs /a 0 .Thus, Newtonian behaviour ensues if M (x) → 1 when x ≫ 1 while MOND behaviour emerges if M (x) → x when x ≪ 1.The modified Poisson equation where G N is Newton's constant and ρ b the mass density of baryons, readily accommodates MOND behaviour as an extension of Newtonian gravity by letting ⃗ a obs = − ⃗ ∇Φ 1 .It has a variational formulation, the AQUAL (Aquadratic Lagrangian) non-relativistic gravity with Lagrangian L = J (y) + 4πGΦρ b where y ≡ |∇Φ| 2 /a 2 0 = x 2 and dJ /dy = M(x), hence, the deep MOND limit requires J → y 3/2 ∝ |∇Φ| 3 [19].Alternatively, MOND can be seen as a modification of inertia [20], though this will not be pursued here.
MOND implies an algebraic relation between the acceleration expected from the baryons only, a N , and the observed acceleration a obs .In the context of rotation curves of disk galaxies this relation is known as the Radial Acceleration Relation (RAR) and it is observationally supported by the rotation curves of a large sample of disk galaxies [21]; see [22] and [23] for recent studies.The observational support for the RAR has recently been extended by about two orders of magnitude in acceleration using weak lensing [24] in place of baryonic tracers, consistent with the deep MOND scaling of the acceleration a obs ∼ √ a N to a previously unexplored range.It has recently been shown that the RAR is a unique fundamental relation in galaxies and that all correlations present in galaxies are manifestations of the RAR [25].
It is commonly accepted that MOND fails in accounting for the profiles of galaxy clusters [26].The amount of gravitational force needed to support the observed profiles is more than MOND predicts.Radial Acceleration Relations for galaxy clusters have also been constructed using weak and strong lensing [27], galaxy kinematics [28] and a combination of X-ray and thermal Sunyaev-Zel'dovich observations [29].In all cases, the galaxy cluster RAR is above the MOND RAR for disk galaxies.For small enough accelerations, at the outskirts of galaxy clusters, observations suggest that the galaxy cluster RAR even goes below the MOND RAR of disk galaxies.In the context of MOND, this suggests an extra dark component with an apparent negative density in the outskirts.
Extensions of GR, limiting to MOND behaviour in the quasistatic, weak-field regime of galaxies, are necessary if one is to further test this line of investigation with gravitational lensing, and with cosmology.GR extensions along these lines have been constructed throughout the years [19,[30][31][32][33][34][35][36]; see [37] for a review.Until recently the most studied candidate has been the Tensor-Vector-Scalar (TeVeS) [32], which, however, is no longer observationally viable on account of the CMB anisotropies [38] and the absence of a significant delay between gravitational waves and the electromagnetic counterpart [39,40] which TeVeS predicts.
Using similar field content as in TeVeS, namely a scalar ϕ and unit-timelike vector field A µ , the Aether-Scalar-Tensor (AeST) theory has recently been proposed [41].It incorporates MOND-like behaviour in its quasi-static limit on small scales, has a ΛCDM on the largest scales, i.e., effective cosmic dust behaviour compatible with observations of the CMB anisotropies and large-scale structure, and has gravitational tensor-mode waves which propagate with the speed of light.Further studies of the AeST theory have been performed in [42][43][44][45][46][47][48][49][50][51].
The phenomenology of the spherically symmetric quasi-static weak-field AeST equations has been studied in [49] in the case of uniform density and Hernquist profile sources as well as in vacuum outside such source profiles.In cases which depart from from spherical symmetry, for instance, disk galaxy solutions, the curl of the vector field can also play a role and the weak-field AeST equations which include this have also been considered [50].Here we extend the previous studies by investigating the weak-field limit of AeST for spherically symmetric systems, having in mind the application to models of galaxy clusters.As such, we consider a simplified model of the galaxy cluster, the (self-gravitating) hydrostatic isothermal gas sphere.
The article is organized as follows.In Section 2.1, we present the AeST action and we derive the weak-field equations in Section 2.3.We then reduce the resulting two-field system to a single equation for the gravitational potential in Section 2.3.3.The reduced equation is a modified Helmholtz equation which includes a mass term in addition to the modified Laplacian term of MOND.We consider the problem of the apparent singularities due to the oscillatory behaviour that the new mass term induces in Section 2.3.5.In Section 3, we find that the same system can be described in Hamiltonian form in a way which is manifestly singularity free.The system can therefore be solved and transformed back to the original variables.We illustrate the procedure for the vacuum case in Section 3.2.1, the case of a prescribed source in Section 3.3, and finally the case of the hydrostatic isothermal source in Section 3.4.In Section 3.5 we construct the RAR in AeST for the case of a hydrostatic isothermal source and uncover novel features which are not present in MOND.

Field content and the theory Lagrangian density
The field content of AeST is the spacetime metric g µν (tensor field), a scalar field ϕ and a unit time-like vector field A µ , i.e., one for which g µν A µ A ν = −1.The time-like vector field A µ allows the construction of kinetic terms not otherwise present in a local, purely metric theory.A projection of the scalar kinetic term onto the preferred time-like direction A µ gives the scalar and a projection onto directions perpendicular to A µ , facilitated by the projector g µν + A µ A ν , gives the scalar Another (projected) kinetic term is the acceleration of A µ which is its gradient projected onto the vector A µ itself The total Lagrangian density is given by the sum of the AeST and matter Lagrangian densities L AeST and L matter respectively, that is L = L AeST + L matter .The AeST Lagrangian density is the Einstein-Hilbert Lagrangian density of g µν , to which matter couples minimally; the Maxwell-like kinetic term F µν F µν associated with A µ ; a coupling between the projected vector gradient J µ and ∇ µ ϕ; and an unspecified function F (Y, Q) of the scalars Y and Q. Explicitly we have where G is the bare gravitational constant, K B is a coupling constant taking values between 0 < K B < 2 [41,42], and λ is a Lagrange multiplier field that enforces the unit time-like constraint g µν A µ A ν = −1.

Forms of the free function F
The form of F determines the dynamics of AeST in a cosmological setting and the effective profiles of the gravitational potential in the quasistatic weak-field limit.In order to recover ΛCDM behaviour, F has features akin to shift-symmetric K-essence [52] which is the low energy limit of ghost condensation [53,54].This means that on a homogeneous-isotropic spacetime, where Y = 0, we may define the function K(Q) ≡ − 1 2 F(0, Q) which is required to have a Taylor expansion of the form . .around a non-zero constant Q 0 , while K 2 is another constant and where (. ..) denote higher order terms in this expansion.Assuming a synchronous coordinate system where g 00 = −1, this ensures that cosmologically Q evolves with redshift z as Q = Q 0 +I 0 (1+z) 3 +. . .where (. ..) denote higher powers of 1+z which are relevant in the early Universe, and I 0 is an initial condition which sets the initial displacement of Q away from its minimum at Q 0 .These considerations on K (and thus on F), result in AeST providing cosmological energy density scaling as (1 + z) 3 plus corrections scaling with higher powers of 1 + z, that is, approximately that of dust.The evolution of Q initially displaced away from Q 0 is to tend towards it, and the cosmological dust density Ω AeST is set by the initial displacement of Q from Q 0 [41].
On quasistatic backgrounds in the late Universe, Q may be assumed to be sufficiently close to its minimum up to small fluctuations.This introduces another reduction of F, the function J (Y) ≡ F(Y, Q 0 )/(2 − K B ), which is subsequently found to appear in the modified Poisson equation and must be appropriately chosen to lead to MOND and to Newton under appropriate conditions.The former (MOND) can happen provided as Y → 0, (that is, when √ Y ≪ a 0 ) which captures the AQUAL requirement |Y| 3/2 (see below).Newtonian behaviour results if as Y → ∞ (that is, when √ Y ≫ a 0 ), where (. ..) denote subleading terms to Y.Here λ s > 0 is a free parameter defined to be part of F; see [42].The effect of λ s is to completely screen the MOND contributions in the large gradient limit as λ s → ∞.
With the above considerations we assume that F consists of two pieces K and J , and that the former has the Taylor expansion discussed above.Then, (2.7) The exact form of J (Y) determines the MOND interpolation function M(x) and as we show below, in spherically symmetric situations, or when curls are subdominant, the simple choice of that we referred to as "totally screened" interpolation function, corresponds to the MOND interpolation function where x ≡ |∇Φ| /a 0 , in the modified Poisson equation The function M in (2.9) has the expected limits M(x) → x for x ≪ 1 (deep MOND regime) and the Newtonian limit M(x) → 1 when x ≫ 1.The exact relation between the function J in the Lagrangian density and the MOND interpolation function M is made explicit in Section 2.3.2 and Section 2.3.3.We refer to Appendix A for an example of a choice of J (Y) for the case of the MOND interpolation function M(x) = x/ (x + 1).Here, we keep to the totally screened choice of (2.8) when discussing our numerical results for the remainder of the article.A generalization of this is found in Appendix A.

2.3
The weak-field, quasi-static limit of AeST 2.3.1 Weak-field limit and the quasistatic Lagrangian density We are interested in the weak-field quasistatic limit of AeST applied to the case of spherical symmetry in order to eventually build simple models of galaxy clusters.We assume that the spacetime is approximately flat, plus small metric fluctuations parametrised in the Newtonian gauge by the potentials Ψ and Φ such that where γ ij is the flat Euclidean metric in arbitrary spatial coordinates.
The unit time-like constraint g µν A µ A ν = −1 implies that and the spatial part of A µ denoted ⃗ A, is decomposed into the gradient of a scalar α and a curl ⃗ α ⊥ so that with ⃗ ∇ • ⃗ h = 0.In the following we set ⃗ h = 0 as it is not relevant to the spherically symmetric situations we study, and is expected to be subdominant even in cases corresponding to disk galaxies; see for example [50].We assume that the quasistatic systems we consider are formed in the late universe where the cosmological contribution of AeST has settled very close to the minimum of the function K, that is Q → Q 0 , up to small fluctuations.Then we may expand ϕ about this minimum so that ϕ = Q 0 t + φ. (2.14) Finally, we impose that all variables are time independent so that Ψ = Φ = α = φ = 0. Defining the gauge invariant variable χ ≡ φ + Q 0 α (see [42]) that combines the scalar perturbation φ and the gradient mode α, the (spatial) kinetic scalar Y from (2.2) equals Hence J (Y) = 2Y 3/2 /3a 0 in (2.8) becomes J (Y) = 2| ⃗ ∇χ| 3 /3a 0 .Let us now consider the Lagrangian density (2.4) and expand the various terms.The Ricci scalar itself only involves the metric potentials Ψ and Φ, and after partial integrations which remove second spatial derivatives we find From the vector field kinetic term comes only F µν F µν = −2| ⃗ ∇Ψ| 2 and from the vector-scalar coupling To calculate the contribution from the matter Lagrangian we set T 00 = ρ b to be the baryon density, while Inserting all these terms into (2.4) leads to where we have multiplied through by 16π G.The variation of Φ gives the condition Φ = Ψ, which is substituted back to give where we have defined the mass scale and the rescaled gravitational constant G as

The weak-field system of equations
The resulting equations, obtained by variations of (2.17) with respect to Φ and χ, respectively, are where we have defined the potential Φ through and where J Y ≡ dJ /dY.The significance of the Φ combination is that as long as the term proportional to µ is small, Φ is the usual potential obtained by solving the Poisson equation sourced only by baryons.Starting from (2.20) and (2.21), we obtain solutions for Φ and χ by casting Φ = Φ + χ and this last equality is used to determine Φ from which particle accelerations ⃗ a (and trajectories) are calculated via ⃗ a = − ⃗ ∇Φ.We refer to the set (2.20) and (2.21) as the two-component system of equations.
Let us now return to the MOND and Newtonian limits.For this it is sufficient to set µ → 0. If J = λ s Y (Newtonian limit) according to (2.6), then (2.21) results in ⃗ ∇ 2 Φ = λ s ⃗ ∇ 2 χ and so up to curl modes which are subdominant and unimportant in spherical symmetry, χ = Φ/λ s , so that (2.20) turns into a Poisson equation for Φ, i.e. ⃗ ∇ 2 Φ = 4πG N ρ with a rescaled Newton's constant corresponding to its measured value as With the above definition of G N , it can be easily checked that J = 2λs 3(1+λs)a 0 |Y| 3/2 according to (2.5) corresponds to a deep MOND limit.
Notice that as λ s → ∞ in our short discussion above, we have that χ → 0 so that Φ → Φ and G N → G.In this sense, the effects of χ are totally screened and so we refer to λ s as the screening parameter.In practice, J is not exactly equal to λ s Y, but contains subdominant terms in Y (such as √ Y, 1/Y or log(Y)).Therefore, screening is not perfect, in the sense that χ ′ is not exactly zero but contains terms scaling as r n with n > −2, so that Φ ′ always dominates as r → 0; it is in the limit r → 0 that screening becomes exact.Thus the λ s → ∞ limit is to be taken that J does not contain a term proportional to Y in this limit.We refer the reader to Appendix A for a more detailed discussion.

Reduction to one equation
It is possible to reduce the two-component system of equations (2.20) and (2.21) to one PDE of the gravitational potential Φ provided that curl components resulting after integrating (2.21) are small.Defining where β = β (|∇χ|/a 0 ), the system (2.20) and (2.21) becomes Integrating the last equation, (2.26), leads to where ⃗ ∇ × ⃗ h is a curl mode.Thus provided that boundary conditions allow ⃗ h to be small, we can eliminate χ in favour of Φ.Then we have that and, after defining We are after y, so isolating it gives and substituting back into (2.27)(with ⃗ h = 0) gives where I(x) ≡ β−1 (x)/x.Taking the divergence of (2.31) results in which allows us to eliminate ⃗ ∇ 2 χ from (2.25) to get the final single-field equation In getting (2.33) we have defined the M = M(x) via the rescaled mass parameter as and the inverse screening parameter as We recognise (2.33) to be the AQUAL equation (1.1) of [19], with the additional term µ 2 Φ, turning it into a modified Helmholtz equation.The new term µ 2 Φ becomes important at very large distances away from sources and introduces novel oscillatory features in the solutions which manifest after the MOND limit occurs [49].

The interpolating function
Consider the simple function (2.8), which implies that λ s = ∞ (and hence β 0 = 0).We first compute from which we find x(y) by forming β(y) = yβ = y + y 2 = x from (2.29), leading to This in turn gives . from which we get the MOND interpolating function using (2.34) One can easily check that it has the correct limits: when x ≫ 1 we have that M(x) → 1 which is the large gradient limit giving Newtonian behaviour, and when x ≪ 1 we have that M (x) → x leading to the low gradient limit which contains MOND when µ 2 → 0.
In Appendix A we compute a more general interpolation function and show that it coincides with (2.39) as λ s → ∞ implying full consistency.

Solving the modified Helmholtz equation in spherical symmetry
We are interested in solving (2.33) in the case of spherical symmetry, for a variety of source profiles ρ b (r), including the isothermal case.In spherical symmetry, (2.33) becomes where ρ b = ρ b (r) and x = |Φ ′ |/a 0 , x being the argument of M and primes denoting derivatives w.r.t.r.
In the large gradient limit where M → 1, (2.40) turns into the Helmholtz equation which has exact oscillatory solutions for any source ρ b (r) given by the homogeneous solution Φ(r) = C 1 cos (μr) /r + C 2 sin (μr) /2μr.For more general M, the oscillations persist but an analytic solution is no longer possible.Numerical solutions to the two-component system, which as we have shown above are equivalent to (2.40), have been studied thoroughly in [49]; see also [50].Assuming a source of mass M and that an interpolation function M is chosen with the correct small gradient limit in order to lead to MOND when µ is negligible (which happens at a scale then oscillations begin [49] at a scale r C given by where rM = λ s r M /(1 + λ s ) and ∆ is a parameter related to the boundary condition which sets Φ on the source surface R * .The above equation holds if the size of the source R * is smaller than r M and so, the oscillations of Φ (and also χ) are occurring in vacuum.We show below that oscillations persist even if the source is extended, such as the case of an isothermal source, however, the scale r C is modified and we have been unable to find an analytic estimate as in (2.42) in that case.
In practice, oscillations in (2.33) for a general M cause problems because the derivative from the divergence acts on |Φ ′ |.Since Φ oscillates it is unavoidable that |Φ ′ | = 0, introducing a singularity when solving for Φ ′′ in (2.40), which we shall see is only apparent.
In the next section we demonstrate that one can find a reduced Lagrangian from which (2.40) is derived by a variational principle.Then, one can find the Hamiltonian and recast (2.40) into a pair of (first order) Hamilton's equations for Φ and its canonical momentum P Φ .The Hamiltonian system does not exhibit these apparent singularities and integration can proceed smoothly.
When an equivalent Hamiltonian system is not known, then it is better to use the twocomponent system of equations rather than (2.40), which in spherical symmetry are 3 Reduced Hamiltonian systems

The master Lagrangian
Consider the Lagrangian where S = S(Φ, r) is an unspecified function and P = P(x) is defined by with x ≡ |Φ ′ |/a 0 as before and M(x) a given interpolation function.The Euler-Lagrange equations resulting from (3.1) lead directly to (2.40), with the density given by In order to resolve the issues with the apparent singularities appearing in (2.40) we choose to work in phase space coordinatised by the potential Φ and its canonical momentum P Φ .We find that which, defining z ≡ |P Φ |/a 0 , may be inverted to give Then the corresponding Hamiltonian H ≡ Φ ′ P Φ − L is From Hamilton's equations we then find the equations of motion as where x = x(z) in (3.8) is evaluated from (3.5).Note that (3.8) is consistent with (3.4) and becomes trivial when substituting x = |Φ ′ |/a 0 .One can now solve the first-order system (3.8) and (3.9) or form the second-order ODE for P Φ as Notice that no apparent singularities occur when the Hamiltonian formulation is used.Let us also note that for our simple choice of interpolation function (2.39) we have that Hence in that case we find 3.2 Vacuum case: S = 0

Vacuum equations and dimensionless variables
Outside a source, we may impose the vacuum condition for which S = 0. Then the system (3.8) and (3.9) becomes and equivalently, the second order equation (3.12) becomes To solve the above equation we impose boundary conditions translated from the original problem to this one using the definition of P Φ (3.4) and its evolution equation (3.14) as In the µ = 0 case, we can readily solve the above system to get that P Φ is a constant from (3.14), so that the force Φ ′ becomes a sum of a Newtonian force 1/r 2 and a MOND force 1/r.
If µ ̸ = 0 we cannot find analytic solutions and a numerical implementation must be employed.We leave the details for Appendix B.2 and here we present our results.
It is easier to work with dimensionless variables adapted to our problem.For the case of the point source with mass M , a good choice is the MOND radius r M defined by (2.41) so that the dimensionless radius is r = r/r M .Then r = 1 is the turning radius from Newtonian behaviour to MOND so that radii larger that r = 1 are in the low gradient regime.Likewise we define the dimensionless mass parameter μ2 ≡ μ2 r and equivalently Note that this means that the case of a galaxy and a galaxy cluster are cases with different μ as their MOND radii are different.For the point source the parameter μ2 is proportional to the mass of the system.

Solutions
As an example we set μ2 = 10 −2 .We consider the gravitational potential in vacuum near a point source and give it Newtonian boundary conditions at r0 = 10 −3 (where it is expected to be Newtonian) so that Φ = −1/r 0 = −10 3 and d Φ/dr = 1/r 2 0 = 10 6 .These boundary conditions for Φ and Φ′ are then translated into boundary conditions for PΦ and P ′ Φ using (3.16) and (3.17), or equivalently, (3.13) and (3.14).The numerical solution to (3.20) of PΦ and its translation to Φ and Φ′ is shown in Figure 1.In the Newtonian regime, Φ ′ ∝ 1/r 2 and since Φ ′ = P Φ / a 0 r 2 this implies that P Φ is approximately constant.Also, Φ ∝ −1/r and since P ′ Φ = −µ 2 r 2 Φ this implies that P ′ Φ ∝ r.In the MOND regime P Φ is also approximately constant but the growth of Φ, which is only logarithmic, does not suppress the r 2 growth of P ′ Φ .These two regimes are not visible in the upper left graph of Figure 1 due to the linear axis but are visible in the upper right graph of Figure 1 where the axes are logarithmic.
The plots of Φ and Φ ′ after their translation from P Φ and P ′ Φ using (3.19) and (3.18) are in the centre and lower left and right parts of Figure 1, respectively.They are here compared with the numerical solutions of the two-component system given by (2.43)The graph of Φ ′ has a vertical tangent at its zeros.This is due to the vertical tangent of the function Sign PΦ | PΦ | at PΦ = 02 , as can be seen in (3.18).
The phase portrait of Φ and PΦ is shown on the left in Figure 2, where Φ has been scaled to r 2 Φ.
The numerical solution of Φ shows that, after a 1/r 2 (Newtonian) decay, and a 1/r (MOND) decay as the MOND radius r = 1 is approached, there is a new regime with oscillatory decay in Φ which begins at a scale r C determined by µ and initial conditions; see (2.42).This regime and its dependence on boundary conditions, interpolation function and AeST parameters µ and λ s has been extensively studied in [49].A logarithmic plot of the absolute value |Φ ′ |/a 0 versus r reveals the trends of the decay, which numerical evidence suggests is a power-law decay of the envelope.This is shown on the right in Figure 2.  where the dimensionless density3 ρb = 4π G 3 M/a 3 0 ρ b , while (3.8) remains the same as it is independent from S. The second-order equation (3.12) for P Φ is then Curiously, in the case of a uniform density sphere of density ρ 0 and radius R * , since ρ ′ b = 0 both inside and outside, the evolution is the same as the vacuum case, only the inside and outside solutions must be matched using (3.21) so as to have Φ continuous.This requires that there is a discontinuity in P ′ Φ so that and hence Although the properties of this case have been extensively treated in [49], it is illustrative to consider, as the final self-gravitating isothermal case shares similar qualitative features.The case of two uniform spheres with different densities and radii, but same total mass, are shown in Figure 3 for Newtonian gravity, MOND and AeST.As more mass is enclosed the acceleration increases until all mass is enclosed.This happens abruptly for a uniform sphere, here at r = 1 and r = 10 −2 , respectively.If the uniform sphere is of low enough density (left figure) there is a low gradient regime throughout, coinciding with the prediction from AQUAL (since in this case R * ≪ µ −1 and so the mass parameter does not play a role).This is evident by the different slopes of the Newtonian (black) and MOND (blue) accelerations after all mass is enclosed.For high-density spheres the MOND regime is only found near r ∼ 1 where the slopes of the Newtonian and MOND accelerations are found to differ.The acceleration in AeST follows that of MOND until a critical radius r C dependent on the gravitational potential, studied quantitatively in [49], after which the acceleration becomes oscillatory.That the transition is dependent on the matter distribution is illustrated in Figure 3 where the first node of the oscillation is found to occur earlier for the high-density case (right figure).

Hydrostatic equilibrium for isothermal profiles
In the case of a gas in hydrostatic equilibrium with temperature T , the pressure gradient ⃗ ∇p balances the gravitational force −ρ ⃗ ∇Φ.If the gas is ideal then p = nk B T where n is the number density and k B is Boltzmann's constant.Expressing the number density in terms of the mass density as ρ = m gas n, and further letting m gas = µ gas m p where m p is the proton mass and µ gas is the mean mass per gas particle in units of the proton mass4 , (not to be confused with AeST parameter µ) allows us to express the pressure w.r.t. the mass density ρ as Then a first order differential equation for ρ is obtained For isothermal profiles that we are here considering, the temperature is constant, so that with the boundary conditions specified at a radius r 0 so that ρ(r 0 ) = ρ 0 and Φ(r 0 ) = Φ 0 .

Hamiltonian formulation
With the above expression of ρ in terms of Φ, we specify appropriate S for the case of isothermal profiles: Once again, (3.8) remains the same as for the previous cases while (3.9) in lieu of (3.30) becomes while (3.12) becomes In the above equation, Φ is expressed in terms of P Φ by solving (3.31), resulting in the Lambert function (product logarithm) W so that where We now have the equation to solve (3.32) and obtain P Φ and the relations (3.6) and (3.33) to translate into Φ ′ and Φ, respectively.

Isothermal profile solutions
We introduce again dimensionless variables, albeit with some differences.First and foremost, we can no longer use the definition of the MOND radius r M from (2.41), as it only holds in vacuum outside some source.This condition is not satisfied here, as we are dealing with a continuous density distribution which extends to infinity (even though in practice, as we show below, the density effectively cuts off at a finite size) and in the Newtonian case it is well known that the total enclosed mass diverges.As such, we need a new scale to use for creating dimensionless variables.Our scale, denoted by r I , should incorporate the constant a 0 so that it is relevant to the problem we are solving.It should also know about the profile (and so c 2 /a 0 , where c is the speed of light cannot be a good choice).A good choice is the parameter ξ −1 which has units of c 2 (same as Φ) and involves the temperature T of the profile.Let us first show that this is a good choice.Starting from (3.32), we define the dimensionless canonical momentum by PΦ = AP Φ for a scale A and consider the μ = 0 case.Then (3.31) leads to so that (3.32) turns into So the choice A = ξ −1 r I and is forced upon us.This choice, also prompts the definition of the dimensionless central density ρ0 ≡ ρ 0 /ρ I with ρ I ≡ a 2 0 ξ/(4πG N ), the dimensionless potential Φ ≡ ξΦ and the dimensionless mass parameter μ ≡ μr I as in the previous sections, so that the dimensionless form of (3.32) is where Φ is calculated from and where Interestingly, r I can also be seen as the distance over which the density decreases by one e-fold in the MOND regime.For the density to decrease one e-fold ξ∆Φ = 1 implying ∆Φ = ξ −1 .
In the MOND regime ∆Φ/∆r ∼ a 0 and so ∆r = ∆Φ/a 0 = 1/ (ξa 0 ) ≡ r I .As we shall see below, it turns out to lead to length scales characteristic of galaxy clusters.As a reference case we shall take T = 5.5 keV and µ gas = 0.6.This provides a characteristic distance scale r I ∼ 0.24 Mpc.This characteristic density ρ M ∼ 1.9•10 −23 kg•m −3 ∼ 2108×ρ crit where ρ crit is the critical cosmic energy density ρ crit = 3H 2 0 /(8πG N ) when assuming H 0 ∼ 70 km s −1 Mpc −1 .Before considering the full solution to (3.38), let briefly us touch upon the case µ = 0 given by (3.36), which in dimensionless form translates to where we have also set μ = 0 into (3.31) to obtain d PΦ /dr = ρ0 r2 e Φ0 − Φ.This equation has no free parameters!Neglecting the first term on the RHS in (3.41) (the term containing | PΦ |) amounts to neglecting the MOND regime, that is, it corresponds to the Newtonian isothermal sphere.The general analytical solution is not known, however, the particular solution PΦ = 2r corresponds to the singular isothermal sphere [55].Neglecting the Newtonian part and focussing on the MOND contribution coming from the term containing the | PΦ | leads to The general solution in this case contains implicit functions, however, a particular solution is found to be (this is (D.11)) which corresponds to a cored density profile controlled by the parameter r0 , given by and potential as The details of the derivation of the pure MOND isothermal profile are found in Appendix D. Interestingly, P Φ is nothing but the mass enclosed within radius r.That is, M = 4π r 0 drr 2 ρ(r) leads to G N M = r 0 dr dP Φ dr = P Φ (r) − P Φ (0).Hence, taking r → ∞, the total mass enclosed is Hence we find that is the MOND radius in disguise!Interestingly though, the MOND radius r M was defined in vacuum outside a source of mass M , while here, we have a continuous source and M is only found after integrating to r = ∞.Nevertheless, the relevant scale remains the same.We now describe the full numerical solution found by integrating (3.38).The numerical procedure is described in Appendix B and for ensuring robustness of results, we compare to the two-component system and demonstrate excellent agreement; see Figure 10    In order to impose boundary conditions, but also to compute Φ which is needed in (3.38), we use the algebraic relations between Φ and d PΦ /dr as well as between d Φ/dr and PΦ .These can be found from the results of the previous sections and are (B.7) and (B.8) respectively.Having Φ and d Φ/dr, we can then reconstruct d Φ/dr and d χ/dr through the interpolation function and further integrate them to obtain Φ and χ, and remember that Φ = Φ + χ.Here the "hats" denote dimensionless quantities by rescaling with ξ as in the case of Φ.
We choose boundary conditions such that Φ(r in ) = 0 and χ(r in ) = χin so that Φin = Φ+ χ = χin .Rather than χin it is sometimes desirable to refer to the final boundary condition as r → ∞, that is, χ∞ , as we use in the figures.The boundary condition for d Φ/dr as well as the connection to the boundary conditions for PΦ and d PΦ /dr are found in appendix B. We set ρ0 = 10 3 which ensures a dense enough state that the high acceleration (Newtonian regime) is valid within the inner part of the isothermal sphere.
The potentials Φ(r), χ(r) and Φ(r) for two cases of μ and two different boundary conditions χ∞ are shown in Figure 4. Just like the vacuum and prescribed source case of the previous sections, beyond a certain radius the potentials oscillate with r.However, here it may happen that the oscillations are around a non-zero asymptotic value, and this depends on the boundary condition.On the left panel of Figure 5 we display the potential gradient |Φ ′ |/a 0 which corresponds to the acceleration felt by a gas particle and on the right panel of Figure 5 we show the gas density profile.It is seen that AeST for shifted potentials χ ∞ < 0, provides a stronger acceleration than MOND.A greater shift in the potential produces a greater deviation from MOND throughout the extended source.However, beyond a certain radius, dependent on μ and χ ∞ , the acceleration goes below the MOND expectation, and then oscillates.It is seen that the gas profiles are significantly more compressed in AeST compared to either Newtonian gravity or MOND.At large radii and small gas densities the oscillations from the potential Φ are imprinted onto the gas density, effectively trapping the gas at the minima and rarefying it at the maxima of the potential.

Effective phantom mass
It is instructive to consider the would-be mass contributed by MOND and AeST, respectively.In the MOND context, this is referred to as the phantom mass.The total apparent mass in the system is defined as the effective source of Φ, if it obeyed the Poisson equation.That is, one determines Φ for a given model (e.g.MOND or AeST), then constructs the r-dependent effective total mass as In order to be able to compare AeST to MOND, we first run an AeST model and compute the corresponding gas density ρ b (r) and potential Φ (AeST) (r).Given this ρ b (r), we use it as a prescribed source in the µ = 0 (i.e.MOND) equations of the previous section and determine the equivalent MOND potential Φ (MOND) (r).From these two potentials, we determine M (AeST) total and M (MOND) total according to (3.48).The enclosed MOND χ-specific part, the phantom mass, is then computed from while the purely AeST part, that is, the contribution from µ 2 , excluding the pure MOND contribution and the gas, is computed as The resulting enclosed mass profiles are shown in Figure 6.M gas (r) M MOND total (r) M AeST total (r) Figure 6.Left panel: The mass profile for the gas M gas (r) for the same models as in Figure 5: AeST (red and green lines), MOND (blue line) and Newtonian gravity (black line).Right panel: The gravitational force of AeST interpreted as an effective mass in Newtonian gravity due to the MOND part of the theory M MOND part (r) (dashed blue line) and the AeST-specific part M AeST part (r) (red dashed line) in addition to the gas mass M gas (r) (black line), giving the total MOND mass M MOND total (r) (gas and phantom mass, full blue line) and the total enclosed AeST mass M AeST total (r) (gas, phantom and µ-mass, full red line), respectively.The main AeST contribution peaks in the outer part.The top axes indicate the physical length in kpc and the right axes the mass in units of the solar mass M ⊙ for the reference case of a T = 5.5 keV gas.
It is seen that the AeST contribution to the apparent mass is more significant in the outskirts, growing steadily until the effects of µ 2 Φ become significant.

The Radial Acceleration Relation of Aether-Scalar-Tensor theory
The Radial Acceleration Relation (RAR) is the relation between the actual acceleration, and the acceleration that would be inferred from just the given matter distribution, in principle visible, in Newtonian gravity.In AeST, there is no universal RAR, as the mass term µ 2 Φ introduces effects that depend on the matter distribution and the boundary condition χ ∞ .
The RAR for the vacuum solutions is shown on the left panel of in Figure 7.A peak goes above the MOND expectation in the case of AeST for potentials shifted by the boundary condition for Φ.The size and position of this peak depend on both the shift in the potential ∆Φ and the dimensionless mass parameter μ.As the square mass parameter μ is proportional to the mass of the source through the dependence of μ2 ≡ µ 2 r 2 M on r M ∝ M it implies that the transition to AeST behaviour is mass-dependent.The peak is made larger by a more negative shift of the potential.Inevitably, beyond a certain range the AeST RAR goes below the MOND RAR.In the context of MOND this has the appearance of a negative mass density.
A similar RAR is present for isothermal gaseous spheres, which can function as a simplified model of a galaxy cluster, as shown on the right panel of Figure 7 for ρ 0 /ρ I = 10 3 and two different choices of μ and two different boundary conditions parametrised by χ∞ .The quantitative features and behaviour with respect to changing μ and χ∞ are similar as in the vacuum case. .The AeST RAR depends on the mass of the system and the shift of the potential ∆Φ (full vs dotted lines).The RAR is shown for Newtonian gravity (trivial, black line with slope 1), MOND (blue line), and various choices of μ and boundary condition for the potential (∆Φ in the vacuum case and χ ∞ in the isothermal case).In both cases the RAR deviates from MOND in the low acceleration regime and may display a peak followed by a sharp descent below it.Smaller values of μ push the peak to smaller accelerations and the peak may disappear altogether, with only the sharp descent remaining, for some boundary conditions.
Interesting behaviour for the isothermal case emerges when the density is so low, that the solution is in the low acceleration regime throughout.In that case, the RAR displays bifurcation such that there are two observed accelerations Φ ′ for a given Newtonian acceleration Φ ′ N .We discuss this in Appendix C.

Conclusion
We studied the weak-field spherically symmetric static solutions in AeST theory, extending previous work done in [49] and [50].We showed that if boundary conditions allow curls to be neglected, such as the case of spherical symmetry, the two-component version of the AeST equations in the weak-field spherically symmetric static limit can be mapped onto a onecomponent equation for the metric potential Φ resembling the standard AQUAL equation plus a mass term which is new and coming from AeST theory.The resulting one-component equation for Φ seems to lead to singularities in the evolution when its gradients | ⃗ ∇Φ| go through zero, which happens in AeST because the mass term µ 2 Φ induces oscillations in Φ as was already reported in [49] and [50].We showed that these are only apparent singularities and nothing physically bad happens at those points.By creating an equivalent Hamiltonian system, and evolving the canonical momentum P Φ rather than Φ, the resulting equation does not have any points where a term may become singular.
We solved the two-component field-based and one-component Hamiltonian system numerically in three cases: in vacuum outside a spherically symmetric source, within a prescribed source and for an isothermal sphere in hydrostatic equilibrium which serves as a simplified model of a galaxy cluster.The vacuum solutions have a Newtonian regime with force scaling as 1/r 2 , an intermediate MOND regime with force scaling as 1/r and an oscillatory phase with an envelope that decays as a power law.For isothermal spheres, the gas profiles in AeST were found to be more compressed than both the Newtonian case and in the case of MOND, signifying a stronger gravitational force than in MOND, or equivalently, more apparent dark matter.In all cases, we found excellent agreement between the two methods.
We constructed the Radial Acceleration Relation (RAR) for the isothermal hydrostatic source.For shifted gravitational potentials we found that there is a peak in the RAR, an enhancement of the acceleration in a range determined by the shift, the weak-field mass parameter of AeST and the mass of the system.The peak in the AeST RAR compared to MOND is followed by a deficit with the MOND expectation, interpreted as a negative phantom mass in MOND.
Our analysis of the isothermal sphere indicates that AeST possesses the qualitative features to address the problem of galaxy clusters in MOND.The (observational) galaxy cluster RAR [27,28,56] displayes a peak above the MOND expectation, with hints of deficit below the MOND expectation (negative phantom mass) in the outskirts.We found such behaviour to happen in AeST, however, it remains to be seen whether this effect can be corroborated with real data after using realistic galaxy cluster models in AeST.A quantitative analysis going beyond the isothermal case, and the presence of multiple components needed for the analysis of a realistic galaxy cluster is left for future work.
In the isothermal profile case, we define in addition the variable σ ≡ − ln (ρ) and evolve the equation This system of ODEs was solved using an explicit Runge-Kutta method of order 5(4) [57] as implemented in SciPy [58].We started the integration at radius r0 with boundary conditions Φ(r in ) = Φin = 0 and χ(r in ) = χin for the potentials, ρin = ρ(r in ) for the density (so that σ(r in ) = − ln (ρ in ) in the isothermal case) and ψ(r in ) = ρin rin /3 − μ2 Φin rin /3 which arises due to the enclosed mass (gas and ) at the position r in .The right condition is that there is zero force at the centre.As there is a coordinate singularity we must evaluate at a position away from the centre, which necessarily encloses a tiny mass (here taken into account).
When a particular asymptotic value χ ∞ is desired, rather than having Φin as boundary condition, the condition χ ≡ χin can first be set to that desired value χin = χ∞ .After one numerical iteration the integration reaches some asymptotic value χn in (indexed by n) and the boundary condition is updated for the next iteration so that ∆ χ(r in ) = χ ∞ − χn in .This is repeated until the asymptotic value of χ equals χ ∞ ; we found 5 iterations usually suffice.

B.2 Reduced Hamiltonian system
In the reduced-one field case, we solve the second-order dimensionless equation for the canonical momentum PΦ given by (3.20) for the vacuum case, (3.22) for the case of prescribed source and (3.38) for the isothermal case, all corresponding to our chosen interpolation function (2.39).Defining Π ≡ d PΦ /dr, The resulting equations in dimensionless form are where in the isothermal case, Φ is obtained from (B.9) below.The equations were solved using an explicit Runge-Kutta method of order 8(5, 3) (DOP853) [59] as implemented in SciPy [58].We specified boundary conditions as in the two-field system above, and translated them into boundary conditions PΦin and Π in using The RAR in AeST for the hydrostatic isothermal gas when densities are low, ρ 0 /ρ I = 1 (in low-acceleration regime throughout).When densities are low there can be two accelerations (bifurcation) Φ ′ in AeST (seen in both red and green lines) corresponding to a given Newtonian acceleration Φ ′ N , when enclosing little mass or when being farther away from the centre.This does not happen in MOND (blue line) which is one-to-one with the Newtonian acceleration.

D.1 Setting up the problem
We would like to solve (3.42) for cases of physical interest.We require that the mass density is positive everywhere, implying also positivity of the mass enclosed (with M (0) = 0. Then (3.35) implies that P ′ Φ ≥ 0, i.e.PΦ cannot decrease, and moreover integrating (3.35) leads to GM (r) = P Φ (r), that is P Φ ≥ 0 with P Φ = 0 only at r = 0.With these considerations, (3.8) implies that Φ ′ ≥ 0 (with strict equality only at r = 0), and thus (3.27) implies that ρ ′ < 0 which means that the density will always decrease as r → ∞.Given that also ρ ≥ 0, then ρ must asymptote to a positive constant as r → ∞, implying also that ρ ′ → 0 as r → ∞.Given that Φ ′ ̸ = 0 except at r = 0, then this implies the stronger condition that ρ → 0 as r → ∞.From (3.35), we then get that P ′ Φ → 0 as r → ∞, hence, PΦ is a monotonically increasing function which asymptotes to a constant as r → ∞.
With the above considerations and since r > 0 we may set t = ln r so that (3.42) becomes an autonomous second order ODE

. 44 )
Details of the setup for the numerical solution of (2.43)-(2.44)are found in Appendix B.1.

2 M
and the dimensionless canonical momentum PΦ = P Φ /G N M and rescaled potential Φ = Φ/ √ G N M a 0 .In these units, (3

Figure 1 .
Figure 1.The momentum frame solutions (first row; left linear, right logarithmic) of PΦ (red line) and P ′ Φ (green line) of (3.20) with clear Newtonian and MOND regimes (black dotted and dashed lines, respectively), translated into the rescaled gravitational potential Φ (second row) and its dimensionless derivative Φ′ (third row) in AeST for μ = 10 −2 .The numerical solution of Φ (red lines, second row) and | Φ′ | and Φ′ (red lines, third row) obtained from the momentum frame of the Hamiltonian system are compared with the solution of the (original) two-component system (2.43) and (2.44) (orange dotted-dashed line) and pure MOND (blue lines of left column).The figures in the centre and lower right column show the oscillations in the region marked by the dashed black lines on the facing figures, but on a linear scales without absolute values.

Figure 2 .
Figure 2. The phase portrait (r 2 Φ, PΦ ) of the AeST Hamiltonian system in vacuum outside a point source with mass M for (non-dimensionalised AeST weak-field mass parameter) μ2 = 10 −2 (left), and the acceleration | Φ′ | (right) in AeST (red line) in vacuum outside the same source vs. radius r/r M ≡ r (where r M = GM/a 0 ), compared to MOND (blue dashed-dotted line), with Newtonian 1/r 2 (dashed black line) and deep MOND 1/r (dashed-dotted black line) trends indicated.

Figure 3 .
Figure3.The acceleration (|Φ ′ |/a 0 ) vs. radius for the case of the uniform density sphere in Newtonian gravity (black line), MOND (blue line) and AeST (for μ = 10 −2 and Φ(r i ) = 0 at r i /r M = 10 −3 , red line) for two different radii R * and densities ρ 0 , low density (left figure) and high density (right figure), respectively, but same total mass.For low enough densities (left figure) there is a MOND regime (black and red lines differ) in the inner part, not present (black and red lines coincide) with higher densities (right figure).The two cases have the same total mass and boundary conditions, yet the transition to oscillatory behaviour occurs earlier for the high-density case.

Figure 5 .
Figure 5. Left panel: Acceleration |Φ ′ |/a 0 vs. radius r/r I for the hydrostatic isothermal gas with initial density ρ 0 /ρ I = 10 3 in Newtonian gravity (black line), MOND (blue line) and the same AeST models as in Figure 4: μ = 10 −1 in red and μ = 10 −2 in green with solid corresponding to χ ∞ = −20 and dotted corresponding to χ ∞ = −10.Right panel: Gas density profiles for the hydrostatic isothermal gas in AeST, compared to MOND and Newtonian gravity for the same models as in the left panel.The density profiles of AeST become more compressed.For the model case {μ = 10 −1 , χ ∞ = −10} one sees the oscillations of the potential imprinted on the density profile.The top axes in both panels indicate the physical length for the reference case of a T = 5.5 keV gas.

Figure 7 .
Figure 7.The AeST Radial Acceleration Relation acceleration for the vacuum case (left) and isothermal sphere (right) with ρ 0 /ρ I = 10 3 .The AeST RAR depends on the mass of the system and the shift of the potential ∆Φ (full vs dotted lines).The RAR is shown for Newtonian gravity (trivial, black line with slope 1), MOND (blue line), and various choices of μ and boundary condition for the potential (∆Φ in the vacuum case and χ ∞ in the isothermal case).In both cases the RAR deviates from MOND in the low acceleration regime and may display a peak followed by a sharp descent below it.Smaller values of μ push the peak to smaller accelerations and the peak may disappear altogether, with only the sharp descent remaining, for some boundary conditions.

2 rVacuum or prescribed source 2 r
Π − μ2 r Sign(P Φ ) | PΦ | + PΦ + r2 dρ b dr Π − ρ0 e Φ0 − Φ + μ2 r Sign(P Φ ) | PΦ | + PΦ Isothermal (B.6) Figure 11.The RAR in AeST for the hydrostatic isothermal gas when densities are low, ρ 0 /ρ I = 1 (in low-acceleration regime throughout).When densities are low there can be two accelerations (bifurcation) Φ ′ in AeST (seen in both red and green lines) corresponding to a given Newtonian acceleration Φ ′ N , when enclosing little mass or when being farther away from the centre.This does not happen in MOND (blue line) which is one-to-one with the Newtonian acceleration.