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

Amel Durakovic, Constantinos Skordis

<figure><img src="image_1.png" /></figure>

► To cite this version:

Amel Durakovic, Constantinos Skordis. Towards galaxy cluster models in Aether-Scalar-Tensor theory: isothermal spheres and curiosities. Journal of Cosmology and Astroparticle Physics, 2024, 04, pp.040. 10.1088/1475-7516/2024/04/040. [hal-04356166](https://hal.science/hal-04356166v1)

HAL Id: hal-04356166

https://hal.science/hal-04356166v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

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

To cite this article: A. Durakovic and C. Skordis JCAP04(2024)040

View the [article online](#) for updates and enhancements.

You may also like

- The Effect of Adiabatic Compression on Dark Matter Halos and the Radial Acceleration Relation
  Pengfei Li, Stacy S. McGaugh, Federico Lelli et al.

- A Strong Falsification of the Universal Radial Acceleration Relation in Galaxies
  Man Ho Chan

- Dwarf Galaxies with Central Cores in Modified Newtonian Dynamics Gravity
  J. Sánchez Almeida
# Towards galaxy cluster models in Aether-Scalar-Tensor theory: isothermal spheres and curiosities

A. Durakovic<sup>1, a,b,*</sup> and C. Skordis<sup>1, a</sup>

<sup>a</sup>CEICO — FZU, Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182 00 Prague 8, Czechia

<sup>b</sup>Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 11 Rue de l'Université, 67000 Strasbourg, France

E-mail: amel@fzu.cz, skordis@fzu.cz

**ABSTRACT:** 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.

**KEYWORDS:** modified gravity, galaxy clusters

ARXIV EPRINT: 2312.00889

<figure><img src="image_3.png" /></figure>

<footnote>*Corresponding author.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>The Aether-Scalar-Tensor theory</td><td>4</td></tr><tr><td>2.1</td><td>Field content and the theory Lagrangian density</td><td>4</td></tr><tr><td>2.2</td><td>Forms of the free function F</td><td>4</td></tr><tr><td>2.3</td><td>The weak-field, quasi-static limit of AeST</td><td>6</td></tr><tr><td>3</td><td>Reduced Hamiltonian systems</td><td>10</td></tr><tr><td>3.1</td><td>The master Lagrangian</td><td>10</td></tr><tr><td>3.2</td><td>Vacuum case: S = 0</td><td>12</td></tr><tr><td>3.3</td><td>Prescribed source profile: S = ρ(r)Φ</td><td>13</td></tr><tr><td>3.4</td><td>Hydrostatic isothermal gas source</td><td>16</td></tr><tr><td>3.5</td><td>The Radial Acceleration Relation of Aether-Scalar-Tensor theory</td><td>21</td></tr><tr><td>4</td><td>Conclusion</td><td>23</td></tr><tr><td>A</td><td>Interpolation functions and screening</td><td>24</td></tr><tr><td>A.1</td><td>Generalization of our totally screened interpolation function (2.9)</td><td>24</td></tr><tr><td>A.2</td><td>Musings on the “simple” MOND interpolation function</td><td>25</td></tr><tr><td>B</td><td>Setup of numerical solutions</td><td>27</td></tr><tr><td>B.1</td><td>Two-component system</td><td>27</td></tr><tr><td>B.2</td><td>Reduced Hamiltonian system</td><td>28</td></tr><tr><td>C</td><td>Low-density RAR</td><td>29</td></tr><tr><td>D</td><td>Pure MOND isothermal sphere</td><td>29</td></tr><tr><td>D.1</td><td>Setting up the problem</td><td>29</td></tr><tr><td>D.2</td><td>The case v₀ = 0</td><td>31</td></tr><tr><td>D.3</td><td>The case v₀ > 0</td><td>32</td></tr></table>

## 1 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–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 $\Lambda$, the resulting $\Lambda$CDM model is in broad agreement with all cosmological observations at the scale of $\sim$ Mpc 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 observed discrepancies. Moreover, there are regularities in the internal dynamics of galaxies permitting their description without a DM component but by assuming that the observed acceleration $a_{\text{obs}}$ is given by $a_{\text{obs}} = \sqrt{a_0}\sqrt{a_N}$ when gravitational accelerations are smaller than a threshold $a_0 \sim 1.2 \times 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–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_{\text{obs}} = a_N$ where $x = a_{\text{obs}}/a_0$. Thus, Newtonian behaviour ensues if $M(x) \to 1$ when $x \gg 1$ while MOND behaviour emerges if $M(x) \to x$ when $x \ll 1$. The modified Poisson equation

$$ \vec{\nabla} \cdot (M \vec{\nabla} \Phi) = 4\pi G_N \rho_b \qquad (1.1) $$

where $G_N$ is Newton’s constant and $\rho_b$ the mass density of baryons, readily accommodates MOND behaviour as an extension of Newtonian gravity by letting $\vec{a}_{\text{obs}} = -\vec{\nabla}\Phi$.<fnref n="1" /> It has a variational formulation, the AQUAL (Aquadratic Lagrangian) non-relativistic gravity with Lagrangian $\mathcal{L} = \mathcal{J}(y) + 4\pi G_N \Phi \rho_b$ where $y = |\nabla\Phi|^2/a_0^2 = x^2$ and $d\mathcal{J}/dy = M(x)$, hence, the deep MOND limit requires $\mathcal{J} \to y^{3/2} \propto |\nabla\Phi|^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_{\text{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, 25] in place of baryonic tracers, consistent with the deep MOND scaling of the acceleration $a_{\text{obs}} \sim \sqrt{a_N}$ to a previously

<footnote n="1">¹Additional non-gravitational forces may be present such that $\vec{a}_{\text{obs}} = \vec{a}_{\text{grav}} + \vec{a}_{\text{non-grav}}$. In that case only $\vec{a}_{\text{grav}} \equiv \vec{\nabla}\Phi$ is affected through (1.1).</footnote>
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 [26].

It is commonly accepted that MOND fails in accounting for the profiles of galaxy clusters [27]. 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 [28], galaxy kinematics [29] and a combination of X-ray and thermal Sunyaev-Zel'dovich observations [30]. 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, 31–37]; see [38] for a review. Until recently the most studied candidate has been the Tensor-Vector-Scalar (TeVeS) [33], which, however, is no longer observationally viable on account of the CMB anisotropies [39] and the absence of a significant delay between gravitational waves and the electromagnetic counterpart [40, 41] which TeVeS predicts.

Using similar field content as in TeVeS, namely a scalar field $\phi$ and unit-timelike vector field $A_{\mu}$, the Aether-Scalar-Tensor (AeST) theory has recently been proposed [42]. It incorporates MOND-like behaviour in its quasi-static limit on small scales, has a $\Lambda$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 [43–52].

The phenomenology of the spherically symmetric quasi-static weak-field AeST equations has been studied in [50] in the case of uniform density and Hernquist profile sources as well as in vacuum outside such source profiles. In cases which depart 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 [51]. 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.
## 2 The Aether-Scalar-Tensor theory

### 2.1 Field content and the theory Lagrangian density

The field content of AeST is the spacetime metric $g_{\mu\nu}$ (tensor field), a scalar field $\phi$ and a unit time-like vector field $A^\mu$, i.e., one for which $g_{\mu\nu}A^\mu A^\nu = -1$. The time-like vector field $A^\mu$ 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^\mu$ gives the scalar

$$Q \equiv A^{\mu} \nabla_{\mu} \phi, \qquad (2.1)$$

and a projection onto directions perpendicular to $A^\mu$, facilitated by the projector $g^{\mu\nu} + A^\mu A^\nu$, gives the scalar

$$\mathcal{Y} \equiv (g^{\mu\nu} + A^\mu A^\nu) \nabla_\mu \phi \nabla_\nu \phi. \qquad (2.2)$$

Another (projected) kinetic term is the acceleration of $A^\mu$ which is its gradient projected onto the vector $A^\mu$ itself

$$J_{\mu} = A^{\alpha} \nabla_{\alpha} A_{\mu}. \qquad (2.3)$$

The total Lagrangian density is given by the sum of the AeST and matter Lagrangian densities $\mathcal{L}_{\text{AeST}}$ and $\mathcal{L}_{\text{matter}}$ respectively, that is $\mathcal{L} = \mathcal{L}_{\text{AeST}} + \mathcal{L}_{\text{matter}}$. The AeST Lagrangian density is the Einstein-Hilbert Lagrangian density of $g_{\mu\nu}$, to which matter couples minimally; the Maxwell-like kinetic term $F_{\mu\nu}F^{\mu\nu}$ associated with $A^\mu$; a coupling between the projected vector gradient $J_\mu$ and $\nabla_\mu\phi$; and an unspecified function $\mathcal{F}(\mathcal{Y}, \mathcal{Q})$ of the scalars $\mathcal{Y}$ and $\mathcal{Q}$. Explicitly we have

$$\mathcal{L}_{\text{AeST}} = \frac{\sqrt{-g}}{16\pi\tilde{G}} \left\{ R - 2\Lambda - \frac{K_B}{2} F_{\mu\nu} F^{\mu\nu} + (2 - K_B)(2J^\mu \nabla_\mu \phi - \mathcal{Y} - \mathcal{F}) - \lambda (A^\mu A_\mu + 1) \right\} \qquad (2.4)$$

where $\tilde{G}$ is the bare gravitational constant, $\Lambda$ is the cosmological constant, $K_B$ is a coupling constant taking values between $0 < K_B < 2$ [42, 43], and $\lambda$ is a Lagrange multiplier field that enforces the unit time-like constraint $g_{\mu\nu}A^\mu A^\nu = -1$.

### 2.2 Forms of the free function $\mathcal{F}$

The form of $\mathcal{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 $\Lambda$CDM behaviour, $\mathcal{F}$ has features akin to shift-symmetric K-essence [53] which is the low energy limit of ghost condensation [54, 55]. This means that on a homogeneous-isotropic space, where $\mathcal{Y} = 0$, we may define the function $\mathcal{K}(\mathcal{Q}) \equiv -\frac{1}{2}\mathcal{F}(0, \mathcal{Q})$ which is required to have a Taylor expansion of the form $\mathcal{K}(\mathcal{Q}) = \mathcal{K}_2(\mathcal{Q} - \mathcal{Q}_0)^2 + \dots$ around a non-zero constant $\mathcal{Q}_0$, while $\mathcal{K}_2$ is another constant and where $(\dots)$ denote higher order terms in this expansion. Assuming a synchronous coordinate system where $g_{00} = -1$, this ensures that cosmologically $\mathcal{Q}$ evolves with redshift $z$ as $\mathcal{Q} = \mathcal{Q}_0 + I_0(1 + z)^3 + \dots$ where $(\dots)$ 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 $\mathcal{Q}$ away from its minimum at $\mathcal{Q}_0$. These considerations on $\mathcal{K}$ (and
thus on $\mathcal{F}$), result in AeST providing cosmological energy density scaling as $(1+z)^3$ with corrections scaling with higher powers of $1+z$, that is, approximately that of dust. The evolution of $\mathcal{Q}$ initially displaced away from $\mathcal{Q}_0$ is to tend towards it, and the cosmological dust density $\Omega_{\text{AeST}}$ is set by the initial displacement of $\mathcal{Q}$ from $\mathcal{Q}_0$ [42].

On quasistatic backgrounds in the late Universe, $\mathcal{Q}$ may be assumed to be sufficiently close to its minimum up to small fluctuations. This introduces another reduction of $\mathcal{F}$, the function $\mathcal{J}(\mathcal{Y}) = \mathcal{F}(\mathcal{Y}, \mathcal{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 Newtonian gravity under appropriate conditions. The former (MOND) can happen provided

$$ \mathcal{J}(\mathcal{Y}) = \frac{2\lambda_s}{3(1+\lambda_s)a_0} |\mathcal{Y}|^{3/2} + \mathcal{O}(\mathcal{Y}^2) \qquad (2.5) $$

as $\mathcal{Y} \to 0$, (that is, when $\sqrt{\mathcal{Y}} \ll a_0$) which captures the AQUAL requirement $|\mathcal{Y}|^{3/2}$ (see below). Newtonian behaviour results if

$$ \mathcal{J}(\mathcal{Y}) = \lambda_s \mathcal{Y} + \dots \qquad (2.6) $$

as $\mathcal{Y} \to \infty$ (that is, when $\sqrt{\mathcal{Y}} \gg a_0$), where $(\dots)$ denote subleading terms to $\mathcal{Y}$. Here $\lambda_s > 0$ is a free parameter defined to be part of $\mathcal{F}$; see [43]. The effect of $\lambda_s$ is to completely screen the MOND contributions in the large gradient limit as $\lambda_s \to \infty$.

With the above considerations we assume that $\mathcal{F}$ consists of two pieces $\mathcal{K}$ and $\mathcal{J}$, and that the former has the Taylor expansion discussed above. Then,

$$ \mathcal{F}(\mathcal{Y}, \mathcal{Q}) = -2\mathcal{K}_2(\mathcal{Q} - \mathcal{Q}_0) + \mathcal{J}(\mathcal{Y}). \qquad (2.7) $$

The exact form of $\mathcal{J}(\mathcal{Y})$ determines the MOND interpolation function $\mathcal{M}(x)$ and as we show below, in spherically symmetric situations, or when curls are subdominant, the simple choice of

$$ \mathcal{J}(\mathcal{Y}) = \frac{2\mathcal{Y}^{3/2}}{3a_0} \qquad (2.8) $$

that we referred to as the “totally screened” interpolation function, corresponds to the MOND interpolation function

$$ \mathcal{M}(x) = \frac{\sqrt{1+4x}-1}{\sqrt{1+4x}+1} \qquad (2.9) $$

where $x = |\nabla\Phi|/a_0$, in the modified Poisson equation

$$ \nabla \cdot [\mathcal{M}(x)\nabla\Phi] + \dots = 4\pi G\rho_b. \qquad (2.10) $$

The function $\mathcal{M}$ in (2.9) has the expected limits $\mathcal{M}(x) \to x$ for $x \ll 1$ (deep MOND regime) and the Newtonian limit $\mathcal{M}(x) \to 1$ when $x \gg 1$. The exact relation between the function $\mathcal{J}$ in the Lagrangian density and the MOND interpolation function $\mathcal{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 $\mathcal{J}(\mathcal{Y})$ for the case of the MOND interpolation function $\mathcal{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 this work. 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 $\Psi$ and $\Phi$ such that

$$ds^2 = -(1+2\Psi)dt^2 + (1-2\Phi)\gamma_{ij}dx^i dx^j \qquad (2.11)$$

where $\gamma_{ij}$ is the flat Euclidean metric in arbitrary spatial coordinates.

The unit time-like constraint $g_{\mu\nu}A^\mu A^\nu = -1$ implies that

$$A_0 = -\sqrt{1+2\Psi}\sqrt{1+(1-2\Phi)^{-1}|\vec{A}|^2} = -1 - \Psi + \frac{\Psi^2}{2} - \frac{|\vec{A}|^2}{2} + \dots \qquad (2.12)$$

and the spatial part of $A_\mu$ denoted $\vec{A}$, is decomposed into the gradient of a scalar $\alpha$ and a curl $\vec{\alpha}_\perp$ so that

$$\vec{A} = \vec{\nabla}\alpha + \vec{\alpha}_{\perp} \qquad (2.13)$$

with $\vec{\nabla} \cdot \vec{\alpha}_\perp = 0$. In the following we set $\vec{\alpha}_\perp = 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 [51]. 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 $\mathcal{K}$, that is $\mathcal{Q} \to \mathcal{Q}_0$, up to small fluctuations. Then we may expand $\phi$ about this minimum so that

$$\phi = \mathcal{Q}_0 t + \varphi. \qquad (2.14)$$

Finally, we impose that all variables are time independent so that $\dot{\Psi} = \dot{\Phi} = \dot{\alpha} = \dot{\varphi} = 0$.

Defining the gauge invariant variable $\chi \equiv \varphi + \mathcal{Q}_0\alpha$ (see [43]) that combines the scalar perturbation $\varphi$ and the gradient mode $\alpha$, the (spatial) kinetic scalar $\mathcal{Y}$ from (2.2) equals

$$\mathcal{Y} = (g^{\mu\nu} + A^\mu A^\nu) \nabla_\mu (\mathcal{Q}_0 t + \varphi) \nabla_\nu (\mathcal{Q}_0 t + \varphi) = |\vec{\nabla}\chi|^2. \qquad (2.15)$$

Hence $\mathcal{J}(\mathcal{Y}) = 2\mathcal{Y}^{3/2}/3a_0$ in (2.8) becomes $\mathcal{J}(\mathcal{Y}) = 2|\vec{\nabla}\chi|^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 $\Psi$ and $\Phi$, and after partial integrations which remove second spatial derivatives we find $\sqrt{-g}R = 2|\vec{\nabla}\Phi|^2 - 4\vec{\nabla}\Phi \cdot \vec{\nabla}\Psi$. From the vector field kinetic term comes only $F_{\mu\nu}F^{\mu\nu} = -2|\vec{\nabla}\Psi|^2$ and from the vector-scalar coupling $J^\mu\nabla_\mu\phi = \vec{\nabla}\Psi \cdot \vec{\nabla}\chi$. The free function $\mathcal{F}$ evaluates to $\mathcal{F} = -2\mathcal{K}_2\mathcal{Q}_0^2\Psi^2 + \mathcal{J}$. To calculate the contribution from the matter Lagrangian we set $T_{00} = \rho_b$ to be the baryon density, while $T_{0i} = T_{ij} = 0$, and use $-2\delta(\sqrt{-g}\mathcal{L}_{\text{matter}}) = \sqrt{-g}T_{\mu\nu}\delta g^{\mu\nu} = 2\rho_b\Psi$.

Inserting all these terms into (2.4) leads to

$$\begin{aligned} \mathcal{L}_{\text{AeST}}^{\text{Quasistatic}} = 2|\vec{\nabla}\Phi|^2 & - 4\vec{\nabla}\Phi \cdot \vec{\nabla}\Psi + K_B|\vec{\nabla}\Psi|^2 + 2(2-K_B)\vec{\nabla}\Psi \cdot \vec{\nabla}\chi \\ & - (2-K_B)(|\vec{\nabla}\chi|^2 + \mathcal{J}) + 2\mathcal{K}_2\mathcal{Q}_0^2\Psi^2 - 16\pi\tilde{G}\rho_b\Psi \end{aligned} \qquad (2.16)$$
where we have multiplied through by $16\pi\tilde{G}$. The variation of $\Phi$ gives the condition $\Phi = \Psi$, which is substituted back to give

$$-\frac{1}{2-K_B} \mathcal{L}_{\text{AeST}}^{\text{Quasistatic}} = |\vec{\nabla}\Phi|^2 - 2\vec{\nabla}\Phi \cdot \vec{\nabla}\chi + |\vec{\nabla}\chi|^2 + \mathcal{J} - \mu^2\Phi^2 + 8\pi G \rho_b \Phi \quad (2.17)$$

where we have defined the mass scale

$$\mu^2 = \frac{2K_2 Q_0^2}{2 - K_B} \quad (2.18)$$

and the rescaled gravitational constant $G$ as

$$G = \frac{2\tilde{G}}{2 - K_B} \quad (2.19)$$

### 2.3.2 The weak-field system of equations

The resulting equations, obtained by variations of (2.17) with respect to $\Phi$ and $\chi$, respectively, are

$$\vec{\nabla}^2 \tilde{\Phi} + \mu^2 \Phi = 4\pi G \rho_b \quad (2.20)$$

$$\vec{\nabla}^2 \tilde{\Phi} = \vec{\nabla} \cdot (\mathcal{J}_Y \vec{\nabla} \chi) \quad (2.21)$$

where we have defined the potential $\tilde{\Phi}$ through

$$\tilde{\Phi} = \Phi - \chi \quad (2.22)$$

and where $\mathcal{J}_Y \equiv d\mathcal{J}/d\mathcal{Y}$. The significance of the $\tilde{\Phi}$ combination is that as long as the term proportional to $\mu$ is small, $\tilde{\Phi}$ 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 $\tilde{\Phi}$ and $\chi$ by casting $\Phi = \tilde{\Phi} + \chi$ and this last equality is used to determine $\Phi$ from which particle accelerations $\vec{a}$ (and trajectories) are calculated via $\vec{a} = -\vec{\nabla}\Phi$. 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 $\mu \to 0$. If $\mathcal{J} = \lambda_s \mathcal{Y}$ (Newtonian limit) according to (2.6), then (2.21) results in $\vec{\nabla}^2 \tilde{\Phi} = \lambda_s \vec{\nabla}^2 \chi$ and so up to curl modes which are subdominant and unimportant in spherical symmetry, $\chi = \tilde{\Phi}/\lambda_s$, so that (2.20) turns into a Poisson equation for $\Phi$, i.e. $\vec{\nabla}^2 \Phi = 4\pi G_N \rho$ with a rescaled Newton's constant corresponding to its measured value as

$$G_N \equiv \frac{1 + \lambda_s}{\lambda_s} G = \frac{1 + \lambda_s}{\lambda_s} \frac{2\tilde{G}}{2 - K_B} \quad (2.23)$$

With the above definition of $G_N$, it can easily be checked that $\mathcal{J} = \frac{2\lambda_s}{3(1+\lambda_s)a_0} |\mathcal{Y}|^{3/2}$ according to (2.5) corresponds to a deep MOND limit.

Notice that, as $\lambda_s \to \infty$ in our short discussion above, we have that $\chi \to 0$ so that $\tilde{\Phi} \to \Phi$ and $G_N \to G$. In this sense, the effects of $\chi$ are totally screened and so we refer to $\lambda_s$ as the screening parameter. In practice, $\mathcal{J}$ is not exactly equal to $\lambda_s \mathcal{Y}$, but contains subdominant terms in $\mathcal{Y}$ (such as $\sqrt{\mathcal{Y}}$, $1/\mathcal{Y}$ or $\log(\mathcal{Y})$). Therefore, screening is not perfect, in the sense that $\chi'$ is not exactly zero but contains terms scaling as $r^n$ with $n > -2$, so that $\Phi'$ always dominates as $r \to 0$; it is in the limit $r \to 0$ that screening becomes exact. We refer the reader to appendix A for a more detailed discussion.
### 2.3.3 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 $\Phi$ provided that curl components resulting after integrating (2.21) are small.

Defining

$$ \beta = 1 + \mathcal{J}_y \qquad (2.24) $$

where $\beta = \beta(|\nabla\chi|/a_0)$, the system (2.20) and (2.21) becomes

$$ \vec{\nabla}^2 \Phi - \vec{\nabla}^2 \chi + \mu^2 \Phi = 4\pi G \rho_b \qquad (2.25) $$

$$ \vec{\nabla} \cdot (\beta \vec{\nabla} \chi) = \vec{\nabla}^2 \Phi. \qquad (2.26) $$

Integrating the last equation, (2.26), leads to

$$ \beta \vec{\nabla} \chi = \vec{\nabla} \Phi + \vec{\nabla} \times \vec{h} \qquad (2.27) $$

where $\vec{\nabla} \times \vec{h}$ is a curl mode. Thus provided that boundary conditions allow $\vec{h}$ to be small, we can eliminate $\chi$ in favour of $\Phi$. Then we have that

$$ \beta |\vec{\nabla} \chi| = |\vec{\nabla} \Phi| \qquad (2.28) $$

and, after defining $y \equiv |\vec{\nabla}\chi|/a_0$ and $x \equiv |\vec{\nabla}\Phi|/a_0$, we define a new function $\tilde{\beta}$ through

$$ \tilde{\beta}(y) \equiv \beta(y) \quad y = x. \qquad (2.29) $$

We are after $y$, so isolating it gives

$$ y = \tilde{\beta}^{-1}(x) \qquad (2.30) $$

and substituting back into (2.27) (with $\vec{h} = 0$) gives

$$ \vec{\nabla} \chi = \tilde{\beta}^{-1}(x) \frac{\vec{\nabla} \Phi}{|\vec{\nabla} \Phi|} \equiv \mathcal{I}(x) \vec{\nabla} \Phi \qquad (2.31) $$

where $\mathcal{I}(x) \equiv \tilde{\beta}^{-1}(x)/x$. Taking the divergence of (2.31) results in

$$ \vec{\nabla}^2 \chi = \vec{\nabla} \cdot (\mathcal{I} \vec{\nabla} \Phi) \qquad (2.32) $$

which allows us to eliminate $\vec{\nabla}^2 \chi$ from (2.25) to get the final single-field equation

$$ \vec{\nabla} \cdot (\mathcal{M} \vec{\nabla} \Phi) + \tilde{\mu}^2 \Phi = 4\pi G_N \rho_b. \qquad (2.33) $$

In getting (2.33) we have defined the $\mathcal{M} = \mathcal{M}(x)$ via

$$ \mathcal{M} \equiv (1 + \beta_0) (1 - \mathcal{I}), \qquad (2.34) $$

the rescaled mass parameter as

$$ \tilde{\mu}^2 \equiv (1 + \beta_0) \mu^2 \qquad (2.35) $$
and the inverse screening parameter as

$$ \beta_0 \equiv \frac{1}{\lambda_s}. \qquad (2.36) $$

We recognise (2.33) to be the AQUAL equation (1.1) of [19], with the additional term $\mu^2\Phi$, turning the modified Poisson equation into a modified Helmholtz equation. The new term $\mu^2\Phi$ becomes important at very large distances away from sources and introduces novel oscillatory features in the solutions which manifest after the MOND limit occurs [50].

### 2.3.4 The interpolating function

Consider the simple function (2.8), which implies that $\lambda_s = \infty$ (and hence $\beta_0 = 0$). We first compute

$$ \beta(y) = 1 + \mathcal{J}' = 1 + y, \qquad (2.37) $$

from which we find $x(y)$ by forming $\tilde{\beta}(y) = y\beta = y + y^2 = x$ from (2.29), leading to

$$ y = \tilde{\beta}^{-1}(x) = \frac{-1 + \sqrt{1 + 4x}}{2}. \qquad (2.38) $$

This in turn gives $\mathcal{I}(x) = \tilde{\beta}^{-1}(x)/x = (-1 + \sqrt{1 + 4x})/2x$. from which we get the MOND interpolating function using (2.34)

$$ \mathcal{M}(x) = 1 - \mathcal{I}(x) = \frac{-1 + \sqrt{1 + 4x}}{1 + \sqrt{1 + 4x}}. \qquad (2.39) $$

One can easily check that it has the correct limits: when $x \gg 1$ we have that $\mathcal{M}(x) \to 1$ which is the large gradient limit giving Newtonian behaviour, and when $x \ll 1$ we have that $\mathcal{M}(x) \to x$ leading to the low gradient limit which contains MOND when $\mu^2 \to 0$.

In appendix A we compute a more general interpolation function and show that it coincides with (2.39) as $\lambda_s \to \infty$ implying full consistency.

### 2.3.5 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 $\rho_b(r)$, including the isothermal case. In spherical symmetry, (2.33) becomes

$$ \frac{1}{r^2} \frac{d}{dr} \left( r^2 \mathcal{M} \frac{d\Phi}{dr} \right) + \tilde{\mu}^2 \Phi = 4\pi G_N \rho_b \qquad (2.40) $$

where $\rho_b = \rho_b(r)$ and $x = |\Phi'|/a_0$, $x$ being the argument of $\mathcal{M}$ and primes denoting derivatives with respect to $r$.

In the large gradient limit where $\mathcal{M} \to 1$, (2.40) turns into the Helmholtz equation which has exact oscillatory solutions in vacuum given by the homogeneous solution $\Phi(r) = C_1 \cos(\tilde{\mu}r)/r + C_2 \sin(\tilde{\mu}r)/2\tilde{\mu}r$. For more general $\mathcal{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 [50]; see also [51].
Assuming a source of mass $M$ and that an interpolation function $\mathcal{M}$ is chosen with the correct small gradient limit in order to lead to MOND when $\mu$ is negligible which happens at a scale

$$r_M = \sqrt{\frac{G_N M}{a_0}}, \qquad (2.41)$$

then oscillations begin [50] at a scale $r_C$ given by

$$r_C \approx \frac{1}{3} \left( \frac{18 \hat{r}_M}{\mu^2} \frac{1}{1 + 3|\Delta|} \right)^{1/3} \qquad (2.42)$$

where $\hat{r}_M = \lambda_s r_M / (1 + \lambda_s)$ and $\Delta$ is a parameter related to the boundary condition which sets $\Phi$ on the source surface $R_*$. The above relation holds if the size of the source $R_*$ is smaller than $r_M$ and so, the oscillations of $\Phi$ (and also $\chi$) 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) for that case.

Oscillations in (2.33) for a general $\mathcal{M}$ cause problems because the derivative from the divergence acts on $|\Phi'|$. Since $\Phi$ oscillates it is unavoidable that $|\Phi'| = 0$, introducing a singularity when solving for $\Phi''$ 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 $\Phi$ and its canonical momentum $P_\Phi$. The Hamiltonian system does not exhibit these apparent singularities and integration can proceed smoothly.

When an equivalent Hamiltonian system is not known, we must resort to using the two-component system of equations rather than (2.40), which in spherical symmetry is

$$\frac{1}{r^2} \frac{d}{dr} [r^2(\Phi' - \chi')] + \tilde{\mu}^2 \Phi = 4\pi G_N \rho_b \qquad (2.43)$$

$$\beta \frac{d\chi}{dr} = \Phi'. \qquad (2.44)$$

Details of the setup for the numerical solution of (2.43)–(2.44) are found in appendix B.1.

# 3 Reduced Hamiltonian systems

## 3.1 The master Lagrangian

Consider the Lagrangian

$$L(\Phi, \Phi', r) = r^2 a_0^2 \mathcal{P} - \frac{\tilde{\mu}^2 r^2}{2} \Phi^2 + 4\pi G_N r^2 \mathcal{S}(\Phi, r) \qquad (3.1)$$

where $\mathcal{S} = \mathcal{S}(\Phi, r)$ is a yet unspecified function and $\mathcal{P} = \mathcal{P}(x)$ is defined by

$$\mathcal{P}(x) \equiv \int^x dx' \mathcal{M}(x') x', \qquad (3.2)$$
with $x \equiv |\Phi'|/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

$$ \rho_b = \frac{\partial S}{\partial \Phi}. \qquad (3.3) $$

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 $\Phi$ and its canonical momentum $P_{\Phi}$. We find that

$$ P_{\Phi} \equiv \frac{\partial L}{\partial \Phi'} = r^2 a_0^2 \frac{dP}{d\Phi'} = r^2 M \Phi' \qquad (3.4) $$

which, defining $z \equiv |P_{\Phi}|/a_0$, may be inverted to give

$$ \tilde{z} \equiv \frac{z}{r^2} = xM \qquad \Rightarrow x = x(\tilde{z}) \qquad (3.5) $$

so that

$$ \Phi' = \frac{x}{z} P_{\Phi}. \qquad (3.6) $$

Then the corresponding Hamiltonian $H \equiv \Phi' P_{\Phi} - L$ is

$$ H(\Phi, P_{\Phi}, r) = \frac{x}{z} P_{\Phi}^2 - a_0^2 z x + r^2 a_0^2 \int^{\tilde{z}} d\tilde{z}' x + \frac{\tilde{\mu}^2 r^2}{2} \Phi^2 - 4\pi G_N r^2 S. \qquad (3.7) $$

From Hamilton's equations we then find the equations of motion as

$$ \Phi' = \mathrm{Sign}(P_{\Phi})a_0x \qquad (3.8) $$

$$ P'_{\Phi} = r^2 \left( -\tilde{\mu}^2 \Phi + 4\pi G_N \frac{\partial S}{\partial \Phi} \right) \qquad (3.9) $$

where $x = x(\tilde{z})$ in (3.8) is evaluated from (3.5). Note that (3.8) is consistent with (3.4) and becomes trivial when substituting $x = |\Phi'|/a_0$. One can now solve the first-order system (3.8) and (3.9) or form the second-order ODE for $P_{\Phi}$ as

$$ P_{\Phi}'' = \frac{2}{r} P_{\Phi}' + r^2 \left[ \left( 4\pi G_N \frac{\partial^2 S}{\partial \Phi^2} - \tilde{\mu}^2 \right) \mathrm{Sign}(P_{\Phi}) a_0 x + 4\pi G_N \frac{\partial^2 S}{\partial \Phi \partial r} \right]. \qquad (3.10) $$

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

$$ x = \sqrt{\tilde{z}} + \tilde{z}. \qquad (3.11) $$

Hence in that case we find

$$ P_{\Phi}'' = \frac{2}{r} P_{\Phi}' + r \left[ \left( 4\pi G_N \frac{\partial^2 S}{\partial \Phi^2} - \tilde{\mu}^2 \right) \left( \mathrm{Sign}(P_{\Phi}) \sqrt{a_0 |P_{\Phi}|} + \frac{1}{r} P_{\Phi} \right) + 4\pi G_N r \frac{\partial^2 S}{\partial \Phi \partial r} \right]. \qquad (3.12) $$
## 3.2 Vacuum case: $S = 0$

### 3.2.1 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

$$ \Phi' = \frac{\mathrm{Sign}(P_{\Phi})}{r} \sqrt{a_0 |P_{\Phi}|} + \frac{1}{r^2} P_{\Phi} \quad (3.13) $$

$$ P'_{\Phi} = -r^2 \tilde{\mu}^2 \Phi \quad (3.14) $$

and equivalently, the second order equation (3.12) becomes

$$ P''_{\Phi} = \frac{2}{r} P'_{\Phi} - \tilde{\mu}^2 \left[ r \mathrm{Sign}(P_{\Phi}) \sqrt{a_0 |P_{\Phi}|} + P_{\Phi} \right]. \quad (3.15) $$

To solve the above equation we impose boundary conditions translated from the original problem to this one using the definition of $P_{\Phi}$ (3.4) and its evolution equation (3.14) as

$$ P_{\Phi}(r_0) = r_0^2 M \left( \frac{|\Phi'(r_0)|}{a_0} \right) \Phi'(r_0) \quad (3.16) $$

$$ P'_{\Phi}(r_0) = -r_0^2 \tilde{\mu}^2 \Phi(r_0). \quad (3.17) $$

In the $\mu = 0$ case, we can readily solve the above system to get that $P_{\Phi}$ is a constant from (3.14), so that the force $\Phi'$ becomes a sum of a Newtonian force $1/r^2$ and a MOND force $1/r$. If $\mu \neq 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 $\hat{r} = r/r_M$. Then $\hat{r} = 1$ is the turning radius from Newtonian behaviour to MOND so that radii larger that $\hat{r} = 1$ are in the low gradient regime. Likewise we define the dimensionless mass parameter $\hat{\mu}^2 \equiv \tilde{\mu}^2 r_M^2$ and the dimensionless canonical momentum $\hat{P}_{\Phi} = P_{\Phi}/G_N M$ and rescaled potential $\hat{\Phi} = \Phi/\sqrt{G_N M a_0}$. In these units, (3.13) and (3.14) become

$$ \frac{d\hat{\Phi}}{d\hat{r}} = \frac{\mathrm{Sign}(\hat{P}_{\Phi})}{\hat{r}} \sqrt{|\hat{P}_{\Phi}|} + \frac{1}{\hat{r}^2} \hat{P}_{\Phi} \quad (3.18) $$

$$ \frac{d\hat{P}_{\Phi}}{d\hat{r}} = -\hat{r}^2 \hat{\mu}^2 \hat{\Phi} \quad (3.19) $$

and equivalently

$$ \frac{d^2 \hat{P}_{\Phi}}{d \hat{r}^2} = \frac{2}{\hat{r}} \frac{d \hat{P}_{\Phi}}{d \hat{r}} - \hat{\mu}^2 \left[ \hat{r} \mathrm{Sign}(P_{\Phi}) \sqrt{|\hat{P}_{\Phi}|} + \hat{P}_{\Phi} \right]. \quad (3.20) $$

Note that the case of a galaxy and a galaxy cluster are cases with different $\hat{\mu}$ as their MOND radii are different. For the point source the parameter $\hat{\mu}^2$ is proportional to the mass of the system.
### 3.2.2 Solutions

As an example we set $\hat{\mu}^2 = 10^{-2}$. We consider the gravitational potential in vacuum near a point source and give it Newtonian boundary conditions at $\hat{r}_0 = 10^{-3}$ (where it is expected to be Newtonian) so that $\hat{\Phi} = -1/\hat{r}_0 = -10^3$ and $d\hat{\Phi}/d\hat{r} = 1/\hat{r}_0^2 = 10^6$. These boundary conditions for $\hat{\Phi}$ and $\hat{\Phi}'$ are then translated into boundary conditions for $\hat{P}_{\Phi}$ and $\hat{P}_{\Phi}'$ using (3.16) and (3.17), or equivalently, (3.13) and (3.14). The numerical solution to (3.20) of $\hat{P}_{\Phi}$ and its translation to $\hat{\Phi}$ and $\hat{\Phi}'$ is shown in figure 1. In the Newtonian regime, $\hat{\Phi}' \propto 1/r^2$ and since $\hat{\Phi}' = P_{\Phi}/(a_0 r^2)$ this implies that $P_{\Phi}$ is approximately constant. Also, $\hat{\Phi} \propto -1/r$ and since $\hat{P}_{\Phi}' = -\hat{\mu}^2 \hat{r}^2 \hat{\Phi}$ this implies that $P_{\Phi}' \propto r$. In the MOND regime $P_{\Phi}$ is also approximately constant but the growth of $\hat{\Phi}$, which is only logarithmic, does not suppress the $r^2$ growth of $P_{\Phi}'$. 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 $\hat{\Phi}$ and $\hat{\Phi}'$ after their translation from $P_{\Phi}$ and $P_{\Phi}'$ 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) and (2.44), and with pure MOND. The graph of $\hat{\Phi}'$ has a vertical tangent at its zeros. This is due to the vertical tangent of the function $\text{Sign}(\hat{P}_{\Phi}) \sqrt{|\hat{P}_{\Phi}|}$ at $\hat{P}_{\Phi} = 0,^2$ as can be seen in (3.18).

The phase portrait of $\hat{\Phi}$ and $\hat{P}_{\Phi}$ is shown on the left in figure 2, where $\hat{\Phi}$ has been scaled to $r^2\hat{\Phi}$.

The numerical solution of $\hat{\Phi}$ shows that, after a $1/r^2$ (Newtonian) decay, and a $1/r$ (MOND) decay as the MOND radius $\hat{r} = 1$ is approached, there is a new regime with oscillatory decay in $\hat{\Phi}$ which begins at a scale $r_C$ determined by $\mu$ and initial conditions; see (2.42). This regime and its dependence on boundary conditions, interpolation function and AeST parameters $\mu$ and $\lambda_s$ have been extensively studied in [50]. A logarithmic plot of the absolute value $|\Phi'|/a_0$ versus $\hat{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.

## 3.3 Prescribed source profile: $S = \rho(r)\Phi$

The case of a (baryonic) source profile of $r$-varying density $\rho_b(r)$ is captured by setting $S = \rho_b(r)\Phi$. Using dimensionless variables as for the vacuum case, the master equation (3.9) turns into

$$ \frac{d\hat{P}_{\Phi}}{d\hat{r}} = \hat{r}^2 (-\hat{\mu}^2 \hat{\Phi} + \hat{\rho}_b) \qquad (3.21) $$

where the dimensionless density$^3$ $\hat{\rho}_b = 4\pi\sqrt{G^3 M/a_0^3} \rho_b$, while (3.8) remains the same as it is independent from $S$. The second-order equation (3.12) for $P_{\Phi}$ is then

$$ \frac{d^2 \hat{P}_{\Phi}}{d \hat{r}^2} = \frac{2}{\hat{r}} \frac{d \hat{P}_{\Phi}}{d \hat{r}} - \hat{\mu}^2 \left( \hat{r} \text{Sign}(P_{\Phi}) \sqrt{|\hat{P}_{\Phi}|} + \hat{P}_{\Phi} \right) + \hat{r}^2 \frac{d \hat{\rho}_b}{d \hat{r}}. \qquad (3.22) $$

Curiously, in the case of a uniform density sphere of density $\rho_0$ and radius $R_*$, since $\rho'_b = 0$ both inside and outside, the evolution is the same as the vacuum case, only the

<footnote>$^2$The second derivative $\Phi''$ is therefore ill-defined which is the apparent singularity problem that the Hamiltonian formulation solves.</footnote>

<footnote>$^3$This means that $\int d\hat{r} \hat{r}^2 \hat{\rho} = 1$.</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Figure 1.** The momentum frame solutions (first row; left linear, right logarithmic) of $\hat{P}_{\Phi}$ (red line) and $\hat{P}'_{\Phi}$ (green line) of (3.20) with clear Newtonian and MOND regimes (black dotted and dashed lines, respectively), translated into the rescaled gravitational potential $\hat{\Phi}$ (second row) and its dimensionless derivative $\hat{\Phi}'$ (third row) in AeST for $\hat{\mu} = 10^{-2}$. The numerical solution of $\hat{\Phi}$ (red lines, second row) and $|\hat{\Phi}'|$ and $\hat{\Phi}'$ (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><img src="image_10.png" /><figcaption>**Figure 2.** The phase portrait ($\hat{r}^2\hat{\Phi}$, $\hat{P}_\Phi$) of the AeST Hamiltonian system in vacuum outside a point source with mass $M$ for (non-dimensionalised AeST weak-field mass parameter) $\hat{\mu}^2 = 10^{-2}$ (left), and the acceleration $|\hat{\Phi}'|$ (right) in AeST (red line) in vacuum outside the same source vs. radius $r/r_M \equiv \hat{r}$ (where $r_M = \sqrt{GM/a_0}$), compared to MOND (blue dashed-dotted line), with Newtonian $1/\hat{r}^2$ (dashed black line) and deep MOND $1/\hat{r}$ (dashed-dotted black line) trends indicated.</figcaption></figure>

inside and outside solutions must be matched using (3.21) so as to have $\Phi$ continuous. This requires that there is a discontinuity in $P'_\Phi$ so that

$$ \lim_{r \to R_*^-} \Phi(r) = \lim_{r \to R_*^+} \Phi(r) \Rightarrow \frac{4\pi G \rho_0}{\mu^2} - \frac{1}{\mu^2 R_*^2} \left. \frac{dP_\Phi}{dr} \right|_{\text{in}} = - \frac{1}{\mu^2 R_*^2} \left. \frac{dP_\Phi}{dr} \right|_{\text{out}} \quad (3.23) $$

and hence

$$ \Delta P'_{\Phi} \equiv \left. \frac{dP_{\Phi}}{dr} \right|_{\text{out}} - \left. \frac{dP_{\Phi}}{dr} \right|_{\text{in}} = -4\pi G \rho_0 R_*^2. \quad (3.24) $$

Although the properties of this case have been extensively treated in [50], 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 $\hat{r} = 1$ and $\hat{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_* \ll \mu^{-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 $\hat{r} \sim 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 [50], 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).
<figure><img src="image_11.png" /><figcaption>**Figure 3.** The acceleration ($|\Phi'|/a_0$) vs. radius for the case of the uniform density sphere in Newtonian gravity (black line), MOND (blue line) and AeST (for $\hat{\mu} = 10^{-2}$ and $\Phi(r_i) = 0$ at $r_i/r_M = 10^{-3}$, red line) for two different radii $R_*$ and densities $\rho_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.</figcaption></figure>

## 3.4 Hydrostatic isothermal gas source

### 3.4.1 Hydrostatic equilibrium for isothermal profiles

In the case of a gas in hydrostatic equilibrium with temperature $T$, the pressure gradient $\vec{\nabla}p$ balances the gravitational force $-\rho\vec{\nabla}\Phi$. If the gas is ideal then $p = nk_BT$ where $n$ is the number density and $k_B$ is Boltzmann's constant. Expressing the number density in terms of the mass density as $\rho = m_{gas}n$, and further letting $m_{gas} = \mu_{gas}m_p$ where $m_p$ is the proton mass and $\mu_{gas}$ is the mean mass per gas particle in units of the proton mass,<fnref n="4" /> (not to be confused with AeST parameter $\mu$) allows us to express the pressure with respect to the mass density $\rho$ as

$$p = \frac{\rho}{\mu_{\text{gas}}m_p}k_B T. \qquad (3.25)$$

Then a first order differential equation for $\rho$ is obtained

$$-\rho\Phi' = p' = \frac{k_B}{\mu_{\text{gas}}m_p}(T\rho' + \rho T'). \qquad (3.26)$$

For isothermal profiles that we are here considering, the temperature is constant, so that

$$-\rho\Phi' = \frac{k_B T}{\mu_{\text{gas}}m_p}\rho' \qquad (3.27)$$

leading to the solution

$$\rho = \rho_0 e^{-\xi(\Phi - \Phi_0)} \qquad (3.28)$$

<footnote n="4">⁴Note that the mean mass per gas particle in units of the proton mass is less than one, which may at first seem counterintuitive, given that both hydrogen and helium have mass equal to or greater than that. Recall then that the gas is ionised and the electrons themselves behave as an ideal gas. Typically $\mu_{gas} = 0.6$.</footnote>
where

$$ \xi = \frac{\mu_{\text{gas}}m_p}{k_B T} \qquad (3.29) $$

with the boundary conditions specified at a radius $r_0$ so that $\rho(r_0) = \rho_0$ and $\Phi(r_0) = \Phi_0$.

### 3.4.2 Hamiltonian formulation

With the above expression of $\rho$ in terms of $\Phi$, we specify appropriate $S$ for the case of isothermal profiles:

$$ S(r, \Phi) = -\frac{\rho_0}{\xi} e^{-\xi(\Phi - \Phi_0)}. \qquad (3.30) $$

Once again, (3.8) remains the same as for the previous cases and (3.9) in lieu of (3.30) becomes

$$ P'_{\Phi} = r^2 \left( -\tilde{\mu}^2 \Phi + 4\pi G_N \rho_0 e^{-\xi(\Phi - \Phi_0)} \right) \qquad (3.31) $$

while (3.12) becomes

$$ P''_{\Phi} = \frac{2}{r} P'_{\Phi} - r \left( 4\pi G_N \rho_0 \xi e^{-\xi(\Phi - \Phi_0)} + \tilde{\mu}^2 \right) \left( \text{Sign}(P_{\Phi}) \sqrt{a_0 |P_{\Phi}|} + \frac{1}{r} P_{\Phi} \right). \qquad (3.32) $$

In the above equation, $\Phi$ is expressed in terms of $P_{\Phi}$ by solving (3.31), resulting in the Lambert function (product logarithm) $W$ so that

$$ \Phi = \frac{1}{\xi} W(u) - \frac{P'_{\Phi}}{\mu^2 r^2} \qquad (3.33) $$

where

$$ u = \frac{4\pi G \rho_0 \xi}{\mu^2} e^{\xi \left( \frac{P'_{\Phi}}{\mu^2 r^2} + \Phi_0 \right)}. \qquad (3.34) $$

We now have the equation to solve (3.32) and obtain $P_{\Phi}$ and the relations (3.6) and (3.33) to translate into $\Phi'$ and $\Phi$, respectively.

### 3.4.3 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 forming 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 $\xi^{-1}$ which has units of $c^2$ (same as $\Phi$) 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 $\hat{P}_{\Phi} = AP_{\Phi}$ for a scale $A$ and consider the $\tilde{\mu} = 0$ case. Then (3.31) leads to

$$ 4\pi G \rho_0 e^{-\xi(\Phi - \Phi_0)} = \frac{A}{r_I^3 \hat{r}^2} \frac{d\hat{P}_{\Phi}}{d\hat{r}} \qquad (3.35) $$
so that (3.32) turns into

$$ \frac{d^2 \hat{P}_\Phi}{d\hat{r}^2} = \frac{2}{\hat{r}} \frac{d\hat{P}_\Phi}{d\hat{r}} - \frac{\xi A}{r_1 \hat{r}} \frac{d\hat{P}_\Phi}{d\hat{r}} \left( \text{Sign}(\hat{P}_\Phi) \sqrt{\frac{r_1^2 a_0}{A} |\hat{P}_\Phi|} + \frac{1}{\hat{r}} \hat{P}_\Phi \right). \quad (3.36) $$

So the choice $A = \xi^{-1}r_1$ and

$$ r_1 = \frac{1}{a_0 \xi} = \frac{k_B T}{\mu_{\text{gas}} m_p a_0} \quad (3.37) $$

is forced upon us. This choice, also prompts the definition of the dimensionless central density $\hat{\rho}_0 \equiv \rho_0/\rho_I$ with $\rho_I \equiv a_0^2\xi/(4\pi G_N)$, the dimensionless potential $\hat{\Phi} \equiv \xi\Phi$ and the dimensionless mass parameter $\hat{\mu} \equiv \tilde{\mu}r_1$ as in the previous sections, so that the dimensionless form of (3.32) is

$$ \frac{d^2 \hat{P}_\Phi}{d\hat{r}^2} = \frac{2}{\hat{r}} \frac{d\hat{P}_\Phi}{d\hat{r}} - \hat{r} (\hat{\rho}_0 e^{\hat{\Phi}_0 - \hat{\Phi}} + \hat{\mu}^2) \left( \text{Sign}(\hat{P}_\Phi) \sqrt{|\hat{P}_\Phi|} + \frac{1}{\hat{r}} \hat{P}_\Phi \right) \quad (3.38) $$

where $\Phi$ is calculated from

$$ \hat{\Phi} = W(u) - \frac{1}{\hat{\mu}^2 \hat{r}^2} \frac{d\hat{P}_{\Phi}}{d\hat{r}} \quad (3.39) $$

and where

$$ u = \frac{\hat{\rho}_0}{\hat{\mu}^2} e^{\frac{1}{\hat{\mu}^2 \hat{r}^2} \frac{d\hat{P}_\Phi}{d\hat{r}}} + \hat{\Phi}_0. \quad (3.40) $$

Interestingly, $r_1$ 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 $\xi\Delta\Phi = 1$ implying $\Delta\Phi = \xi^{-1}$. In the MOND regime $\Delta\Phi/\Delta r \sim a_0$ and so $\Delta r = \Delta\Phi/a_0 = 1/(\xi a_0) \equiv r_1$. 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 $\mu_{\text{gas}} = 0.6$. This provides a characteristic distance scale $r_1 \sim 0.24$ Mpc. This characteristic density $\rho_M \sim 1.9 \cdot 10^{-23}$ kg·m$^{-3} \sim 2108 \times \rho_{\text{crit}}$ where $\rho_{\text{crit}}$ is the critical cosmic energy density $\rho_{\text{crit}} = 3H_0^2/(8\pi G_N)$ when assuming $H_0 \sim 70$ km s$^{-1}$Mpc$^{-1}$.

Before considering the full solution to (3.38), let briefly us touch upon the case $\mu = 0$ given by (3.36), which in dimensionless form translates to

$$ \frac{d}{d\hat{r}} \left( \frac{1}{\hat{r}^2} \frac{d\hat{P}_\Phi}{d\hat{r}} \right) = -\frac{1}{\hat{r}^3} \frac{d\hat{P}_\Phi}{d\hat{r}} \left[ \text{Sign}(\hat{P}_\Phi) \sqrt{|\hat{P}_\Phi|} + \frac{1}{\hat{r}} \hat{P}_\Phi \right] \quad (3.41) $$

where we have also set $\tilde{\mu} = 0$ into (3.31) to obtain $d\hat{P}_\Phi/d\hat{r} = \hat{\rho}_0\hat{r}^2e^{\hat{\Phi}_0 - \hat{\Phi}}$. This equation has no free parameters!

Neglecting the first term on the r.h.s. of (3.41) (the term containing $\sqrt{|\hat{P}_\Phi|}$) 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 $\hat{P}_\Phi = 2\hat{r}$ corresponds to the singular isothermal sphere [56]. Neglecting the Newtonian part and focussing on the MOND contribution coming from the term containing the $\sqrt{|\hat{P}_\Phi|}$ leads to

$$ \frac{d}{d\hat{r}} \left( \frac{1}{\hat{r}^2} \frac{d\hat{P}_\Phi}{d\hat{r}} \right) = -\frac{1}{\hat{r}^3} \frac{d\hat{P}_\Phi}{d\hat{r}} \text{Sign}(\hat{P}_\Phi) \sqrt{|\hat{P}_\Phi|}. \quad (3.42) $$
The general solution in this case contains implicit functions, however, a particular solution is found to be (this is (D.11))

$$ \hat{P}_{\Phi} = \frac{81 \hat{r}^3 / \hat{r}_0^3}{4 \left( \sqrt{\hat{r}^3 / \hat{r}_0^3 + 1} \right)^2} \quad \forall \ \hat{r} > 0 \qquad (3.43) $$

which corresponds to a cored density profile controlled by the parameter $\hat{r}_0$, given by

$$ \hat{\rho} = \frac{243 a_0^2 \xi}{16 \pi G_N \hat{r}_0^3 [1 + (\hat{r}/\hat{r}_0)^{3/2}]^3} \qquad (3.44) $$

and potential as

$$ \Phi = \Phi_0 + \frac{1}{\xi} \ln \frac{4\hat{\rho}_0 \hat{r}_0^3}{243} + \frac{3}{\xi} \ln [1 + (\hat{r}/\hat{r}_0)^{3/2}]. \qquad (3.45) $$

The details of the derivation of the pure MOND isothermal profile are found in appendix D.

Interestingly, $P_{\Phi}$ is nothing but the mass enclosed within radius $r$. That is, $M = 4\pi \int_{0}^{r} dr r^2 \rho(r)$ leads to $G_N M = \int_{0}^{r} dr \frac{dP_{\Phi}}{dr} = P_{\Phi}(r) - P_{\Phi}(0)$. Hence, taking $\hat{r} \to \infty$, the total mass enclosed is

$$ M = \frac{81}{4G_N a_0 \xi^2} = \frac{81 r_I^2 a_0}{4G_N}. \qquad (3.46) $$

Hence we find that

$$ r_I = \frac{2}{9} \sqrt{\frac{G_N M}{a_0}} = \frac{2}{9} r_M \qquad (3.47) $$

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 = \infty$. 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 of appendix B.2.

In order to impose boundary conditions, but also to compute $\Phi$ which is needed in (3.38), we use the algebraic relations between $\hat{\Phi}$ and $d\hat{P}_{\Phi}/d\hat{r}$ as well as between $d\hat{\Phi}/d\hat{r}$ and $\hat{P}_{\Phi}$. These can be found from the results of the previous sections and are (B.7) and (B.8) respectively. Having $\hat{\Phi}$ and $d\hat{\Phi}/d\hat{r}$, we can then reconstruct $d\hat{\Phi}/d\hat{r}$ and $d\hat{\chi}/d\hat{r}$ through the interpolating function and further integrate them to obtain $\hat{\Phi}$ and $\hat{\chi}$, recalling that $\hat{\Phi} = \hat{\Phi} + \hat{\chi}$. Here the circumflexes denote dimensionless quantities rescaled using $\xi$ as in the case of $\hat{\Phi}$.

We choose boundary conditions such that $\hat{\Phi}(\hat{r}_{in}) = 0$ and $\hat{\chi}(\hat{r}_{in}) = \hat{\chi}_{in}$ so that $\hat{\Phi}_{in} = \hat{\Phi} + \hat{\chi} = \hat{\chi}_{in}$. Rather than $\hat{\chi}_{in}$ it is sometimes desirable to refer to the final boundary condition as $r \to \infty$, that is, $\hat{\chi}_{\infty}$, as we use in the figures. The boundary condition for $d\hat{\Phi}/d\hat{r}$ as well as the connection to the boundary conditions for $\hat{P}_{\Phi}$ and $d\hat{P}_{\Phi}/d\hat{r}$ are found in appendix B. We set $\hat{\rho}_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.
<figure><img src="image_12.png" /><figcaption>**Figure 4.** Left panel: the AeST potentials $\hat{\chi}(r)$ (blue lines) and $\hat{\Phi}(r)$ (black lines) for the case of the hydrostatic isothermal gas with central density $\hat{\rho}_0 = \rho_0/\rho_I = 10^3$ for $\hat{\mu} = 10^{-1}$ and two choices of the asymptotic value of $\hat{\chi}(r): \hat{\chi}_{\infty} = -20$ (dotted lines) and $\hat{\chi}_{\infty} = -10$ (full lines). Interestingly, the asymptotic values of the two fields are negatives of one another $\hat{\chi}_{\infty} = -\hat{\Phi}_{\infty}$ so that $\hat{\Phi} \sim 0$. Right panel: the gravitational potential $\hat{\Phi}(r) \equiv \hat{\chi}(r) + \hat{\Phi}(r)$ for two different choices of the AeST mass parameter, $\hat{\mu} = 10^{-2}$ (green line) and $\hat{\mu} = 10^{-1}$ (red line), and two different choices of the asymptotic values $\hat{\chi}_{\infty} = -20$ (solid) and $\hat{\chi}_{\infty} = -10$ (dotted). The gravitational potential $\hat{\Phi}(r)$ decays while oscillating about zero. The top axes indicate the physical length for the reference case of a $T = 5.5$ keV ideal gas.</figcaption></figure>

The potentials $\hat{\Phi}(r)$, $\hat{\chi}(r)$ and $\hat{\Phi}(r)$ for two cases of $\hat{\mu}$ and two different boundary conditions $\hat{\chi}_{\infty}$ 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 $|\Phi'|/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 $\chi_{\infty} < 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 $\hat{\mu}$ and $\chi_{\infty}$, 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 $\Phi$ are imprinted onto the gas density, effectively trapping the gas at the minima and rarefying it at the maxima of the potential.

### 3.4.4 Effective phantom mass

It is instructive to consider, from a Newtonian perspective, 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 $\Phi$, if it obeyed the Poisson equation. That is, one determines $\Phi$ for a given model (e.g., MOND or
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 5.** Left panel: acceleration $|\Phi'|/a_0$ vs. radius $r/r_I$ for the hydrostatic isothermal gas with initial density $\rho_0/\rho_I = 10^3$ in Newtonian gravity (black line), MOND (blue line) and the same AeST parameters as in figure 4: $\hat{\mu} = 10^{-1}$ in red and $\hat{\mu} = 10^{-2}$ in green with solid corresponding to $\chi_\infty = -20$ and dotted corresponding to $\chi_\infty = -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 $\{\hat{\mu} = 10^{-1}, \chi_\infty = -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.</figcaption></figure>

AeST), then constructs the r-dependent effective total mass as

$$M_{\text{total}}(r) \equiv -r^2\Phi'/G. \qquad (3.48)$$

In order to be able to compare AeST to MOND, we first run an AeST model and compute the corresponding gas density $\rho_b(r)$ and potential $\Phi^{(\text{AeST})}(r)$. Given this $\rho_b(r)$, we use it as a prescribed source in the $\mu = 0$ (i.e. MOND) equations of the previous section and determine the equivalent MOND potential $\Phi^{(\text{MOND})}(r)$. From these two potentials, we determine $M_{\text{total}}^{(\text{AeST})}$ and $M_{\text{total}}^{(\text{MOND})}$ according to (3.48). The enclosed MOND $\chi$-specific part, the phantom mass, is then computed from

$$M_{\text{part}}^{\text{MOND}}(r) \equiv M_{\text{total}}^{\text{MOND}}(r) - M_{\text{gas}}(r) \qquad (3.49)$$

while the purely AeST part, that is, the contribution from $\mu^2$, excluding the pure MOND contribution and the gas, is computed as

$$M_{\text{part}}^{\text{AeST}}(r) \equiv M_{\text{total}}^{\text{AeST}}(r) - M_{\text{total}}^{\text{MOND}}(r). \qquad (3.50)$$

The resulting enclosed mass profiles are shown in figure 6.

It is seen that the AeST contribution to the apparent mass is more significant in the outskirts, growing steadily until the effects of $\mu^2\Phi$ become significant.

## 3.5 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
<figure><img src="image_15.png" /><figcaption>**Figure 6.** Left panel: the mass profile for the gas $M_{\text{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_{\text{part}}^{\text{MOND}}(r)$ (dashed blue line) and the AeST-specific part $M_{\text{part}}^{\text{AeST}}(r)$ (red dashed line) in addition to the gas mass $M_{\text{gas}}(r)$ (black line), giving the total MOND mass $M_{\text{total}}^{\text{MOND}}(r)$ (gas and phantom mass, full blue line) and the total enclosed AeST mass $M_{\text{total}}^{\text{AeST}}(r)$ (gas, phantom and $\mu$-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_{\odot}$ for the reference case of a $T = 5.5$ keV gas.</figcaption></figure>

visible, in Newtonian gravity. In AeST, there is no *universal* RAR, as the mass term $\mu^2\Phi$ introduces effects that depend on the matter distribution and the boundary condition $\chi_\infty$.

The RAR for the vacuum solutions is shown on the left panel of figure 7. A peak goes above the MOND expectation in the case of AeST for potentials shifted by the boundary condition for $\Phi$. The size and position of this peak depend on both the shift in the potential $\Delta\Phi$ and the dimensionless mass parameter $\hat{\mu}$. As the square mass parameter $\hat{\mu}$ is proportional to the mass of the source through the dependence of $\hat{\mu}^2 \equiv \mu^2 r_M^2$ on $r_M \propto 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 $\rho_0/\rho_1 = 10^3$ and two different choices of $\hat{\mu}$ and two different boundary conditions parametrised by $\hat{\chi}_\infty$. The quantitative features and behaviour with respect to changing $\hat{\mu}$ and $\hat{\chi}_\infty$ are similar as in the vacuum case.

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 $\Phi'$ for a given Newtonian acceleration $\Phi'_{\text{N}}$. We discuss this in appendix C.
<figure><img src="image_16.png" /></figure>

**Figure 7.** The AeST Radial Acceleration Relation acceleration for the vacuum case (left) and isothermal sphere (right) with $\rho_0/\rho_I = 10^3$. The AeST RAR depends on the mass of the system and the shift of the potential $\Delta\Phi$ (full vs dotted lines). The RAR is shown for Newtonian gravity (trivial, black line with slope 1), MOND (blue line), and various choices of $\hat{\mu}$ and boundary condition for the potential ($\Delta\Phi$ in the vacuum case and $\chi_\infty$ 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 $\hat{\mu}$ push the peak to smaller accelerations and the peak may disappear altogether, with only the sharp descent remaining, for some boundary conditions.

## 4 Conclusion

We studied the weak-field spherically symmetric static solutions in AeST theory, extending previous work done in [50] and [51]. 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 one-component equation for the metric potential $\Phi$ resembling the standard AQUAL equation plus a mass term which is new and coming from AeST theory.

The resulting one-component equation for $\Phi$ seems to lead to singularities in the evolution when its gradients $|\vec{\nabla}\Phi|$ go through zero, which happens in AeST because the mass term $\mu^2\Phi$ induces oscillations in $\Phi$ as was already reported in [50] and [51]. We showed that these are only apparent singularities that are otherwise inconsequential. By finding an equivalent Hamiltonian system, and evolving the canonical momentum $P_\Phi$ rather than $\Phi$, the resulting equation does not have any points where a term may become singular.

We solved both the two-component system and the corresponding 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 respect to 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 [28, 29, 57] 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.

## Acknowledgments

We thank Benoît Famaey, Tom Złośnik, Tobias Mistele and William Barker for useful discussions. AD, and CS in part, were supported by the European Structural and Investment Funds and the Czech Ministry of Education, Youth and Sports: project CoGraDS — CZ.02.1.01/0.0/0.0/15003/0000437. AD was from 2023/03/01 supported by the European Regional Development Fund and the Czech Ministry of Education, Youth and Sports: project MSCA Fellowship CZ FZU I — CZ.02.01.01/00/22_010/0002906.

## A Interpolation functions and screening

In this appendix we present two further interpolation functions: one for which $\lambda_s \to \infty$ results in the totally screened function (2.9) and one corresponding to the “simple” MOND interpolation function.

### A.1 Generalization of our totally screened interpolation function (2.9)

Consider the function

$$ \mathcal{J} = \lambda_s \left\{ \mathcal{Y} - 2a_0(1 + \lambda_s)\sqrt{\mathcal{Y}} + 2(1 + \lambda_s)^2 a_0^2 \ln \left[ 1 + \frac{\sqrt{\mathcal{Y}}}{(1 + \lambda_s)a_0} \right] \right\} \quad (A.1) $$

which upon further inspection has the correct Newtonian and MOND limits given by (2.6) and (2.5) respectively. Taking derivatives we find that

$$ \mathcal{J}_Y = \frac{\lambda_s \sqrt{\mathcal{Y}}}{a_0(1 + \lambda_s) + \sqrt{\mathcal{Y}}}. \quad (A.2) $$

Following the procedure in section 2.3.3 this leads to the interpolation function

$$ \mathcal{M}(x) = \frac{1 + \lambda_s + (1 + 2\lambda_s)x - \sqrt{(\lambda_s + 1)^2 + 2(\lambda_s + 1)(2\lambda_s + 1)x + x^2}}{2\lambda_s x} \quad (A.3) $$
<figure><img src="image_17.png" /><figcaption>**Figure 8.** The interpolating function $M_{\lambda_s}(x)$ vs. $x \equiv |\vec{\nabla}\Phi|/a_0$ in with linear axes (left) and logarithmic axes (right) for several choices of $\lambda_s$.</figcaption></figure>

which has the right MOND ($M \to x$) and Newtonian ($M \to 1$) limits. Now taking the $\lambda_s \to \infty$ limit we find that

$$M(x) \to \frac{1 + 2x - \sqrt{1 + 4x}}{2x} = \frac{-1 + \sqrt{1 + 4x}}{1 + \sqrt{1 + 4x}} \qquad (A.4)$$

which is (2.9). In the same limit, we have that $J_y \to \sqrt{y}/a_0$ (which when integrated would correspond to (2.8)), although $J$ diverges.

In order to quantify the effects of $\lambda_s$ on the interpolation function, we plot it for different choices of $\lambda_s$, including the case $\lambda_s = \infty$ in figure 8. We observe that as $\lambda_s \to \infty$ the interpolation function has a shallower transition from MOND to Newton, while as $\lambda_s \to 0$ the transition becomes somewhat sharper but not quite instantaneous. This behaviour may cause confusion about the exact meaning of screening.

To elucidate the effect of screening, consider the ratio of the scalar field gradient $\vec{\nabla}\chi$ to the gravitational potential gradient $\vec{\nabla}\Phi$. We have that $|\vec{\nabla}\chi|/|\vec{\nabla}\Phi| = I_{\lambda_s}(x) = 1 - M_{\lambda_s}(x)$ and we plot this in figure 9 for several choices of $\lambda_s$. We observe that the $\lambda_s \to \infty$ case is the case where the scalar field contributes fractionally the least to the total gravitational potential. In other words, as one goes to higher gravitational gradients relevant in the strong-field regime, screened cases are the ones where the scalar contribution is mostly suppressed. The amount of suppression increases for higher gradients (smaller scales). On the contrary, if $\lambda_s$ is small, then this ratio tends to a constant value and therefore the scalar always contributes, even as $\vec{\nabla}\Phi \to \infty$ — it is not screened.

## A.2 Musings on the “simple” MOND interpolation function

The simple MOND interpolating function is defined as

$$M(x) = \frac{x}{1+x} \qquad (A.5)$$

Suppose that we would like to find the equivalent $J(y)$ leading to the above function, which is also screened in the deep Newtonian regime (i.e. $\beta_0 = 0$). We then have that
<figure><img src="image_18.png" /><figcaption>**Figure 9.** The ratio of the scalar field gradient to the gravitational potential gradient $|\nabla\chi|/|\nabla\Phi| = I_{\lambda_s}$ vs. $x \equiv |\nabla\Phi|/a_0$.</figcaption></figure>

$I = 1 - x/(1+x) = 1/(1+x)$ so that $y(x) = xI(x) = x/(1+x)$. This already signifies a possible issue because $y$ is bounded; we have that $0 \le y < 1$, where the limit $y \to 1$ is reached as $x \to \infty$. Continuing nevertheless, $y(x)$ is a linear equation with one solution: $x(y) = y/(1-y)$, so that $\beta(y) = x(y)/y = 1/(1-y)$ leading to $\mathcal{J}_Y(y) = y/(1-y)$. Rescaling by setting $y = \sqrt{Y}/a_0$ leads to

$$ \mathcal{J}_Y = \frac{\sqrt{Y}}{a_0 - \sqrt{Y}} \quad \Rightarrow \quad \mathcal{J} = -Y - 2a_0\sqrt{Y} - 2a_0^2 \ln \left( 1 - \frac{\sqrt{Y}}{a_0} \right) \quad (A.6) $$

where we have discarded the irrelevant constant of integration. Despite appearances, $\mathcal{J} \ge 0$ for $0 \le \sqrt{Y} < a_0$, and so it is expected to obey the correct stability requirements in the MOND limit found in [43]. This function has the correct (MOND) limit as $y \to 0$ but is bounded by $Y < a_0^2$ and so does not fall under our ansatz that the Newtonian limit is reached as $Y \to \infty$. This may not be a problem, however, expanding at the limiting value by setting $Y = a_0^2(1-\epsilon)$, leads to $2\ln 2 - 3 - 2\ln \epsilon + 3\epsilon/2 + \dots$. The $\epsilon$ term is welcome, but the logarithmic singularity $\ln \epsilon$ means that a weak-field expansion in the limit $Y \to a_0^2$ is ill-defined, and hence, the Newtonian limit (and stability requirements) are questionable and warrant further investigation.

One could try to fix the above issue with $\ln \epsilon$ by adding the term $2Y + 2a_0^2 \ln(1 - Y/a_0^2)$, which still retains a MOND limit and removes the $\ln \epsilon$ divergence. Unfortunately, this results to $\mathcal{J} \to 2\ln 2 - 1 - \epsilon/2 + \dots$, i.e. $\epsilon$ has the wrong sign which signifies a ghost in the weak-field expansion on a Minkowski background, in the large gradient limit (Newtonian limit). Moreover, the resulting interpolation function $M(x)$ is no longer of the form $x/(1+x)$.

A final consideration is to have a non-zero $\beta_0$ in (2.34). In that case we arrive at

$$ y = x - \frac{1}{1 + \beta_0} \frac{x^2}{1 + x} \quad (A.7) $$

which, unlike the case $\beta_0 = 0$, has two solutions. We choose

$$ \tilde{\beta}(y) = x(y) = \frac{1 + \lambda_s}{2} \left[ -1 + y + \sqrt{(1 - y)^2 + \frac{4y}{1 + \lambda_s}} \right] \quad (A.8) $$
so that from (2.29) and (2.24) we find

$$ \mathcal{J}_y(y) = -\frac{1-\lambda_s}{2} + \frac{1+\lambda_s}{2y} \left[ -1 + \sqrt{(1-y)^2 + \frac{4y}{1+\lambda_s}} \right] \quad (A.9) $$

where $y = \sqrt{|\mathcal{Y}|}/a_0$. This may be further integrated to get $\mathcal{J}(\mathcal{Y})$ analytically but we do not present its exact form here. In the limit $\beta_0 \to 0$, we recover the same functional form as (A.6), provided that $y < 1$. The function (A.9), however, does not have any problem in crossing the $y=1$ line if $\beta_0 \neq 0$, and indeed $y$ can be taken to $\infty$. Moreover, we see that both $\mathcal{J}_y$ and its integral $\mathcal{J}(\mathcal{Y})$ have the correct limits (Newton and MOND) according to (2.6) and (2.5). Thus, (A.9) represents a one-parameter family of functions, depending on the parameter $\beta_0$ (or equivalently $\lambda_s$), from which the same interpolation function (A.5) emerges.

# B Setup of numerical solutions

## B.1 Two-component system

In this section we present the numerical procedure for solving the two-component system of equations and our three cases corresponding to the density $\rho$: (i) $\rho = 0$ (vacuum), (ii) $\rho(r)$ being a prescribed functional form, and (iii) $\rho$ being the solution of the hydrostatic condition (3.27). All our results were produced using our simple choice of function $\mathcal{J}(\mathcal{Y})$ given by (2.8) and we have used dimensionless variables, as defined in the relevant sections for vacuum 3.2.1, prescribed source 3.3, and isothermal profile 3.4, respectively. We evolve the dimensionless potentials $\hat{\Phi}$ and $\hat{\chi}$ according to (2.20) and the integral of (2.21), while $\hat{\Phi} = \hat{\Phi} + \hat{\chi}$ is a derived variable.

Defining the variables $\psi \equiv d\hat{\Phi}/d\hat{r}$, we numerically integrated the following system

$$ \frac{d\hat{\Phi}}{d\hat{r}} = \psi \quad (B.1) $$

$$ \frac{d\psi}{d\hat{r}} = \hat{\rho} - \frac{2}{\hat{r}}\psi - \hat{\mu}^2 \hat{\Phi} \quad (B.2) $$

$$ \frac{d\hat{\chi}}{d\hat{r}} = \frac{\psi}{\sqrt{|\psi|}} \quad (B.3) $$

In the isothermal profile case, we define in addition the variable $\sigma \equiv -\ln(\hat{\rho})$ and evolve the equation

$$ \sigma' = \left( 1 + \frac{1}{\sqrt{|\psi|}} \right) \psi. \quad (B.4) $$

This system of ODEs was solved using an explicit Runge-Kutta method of order 5(4) [58] as implemented in SciPy [59]. We started the integration at radius $\hat{r}_{in}$ with boundary conditions $\hat{\Phi}(\hat{r}_{in}) = \hat{\Phi}_{in} = 0$ and $\hat{\chi}(\hat{r}_{in}) = \hat{\chi}_{in}$ for the potentials, $\hat{\rho}_{in} = \hat{\rho}(\hat{r}_{in})$ for the density (so that $\sigma(\hat{r}_{in}) = -\ln(\hat{\rho}_{in})$ in the isothermal case) and $\psi(\hat{r}_{in}) = \hat{\rho}_{in}\hat{r}_{in}/3 - \hat{\mu}^2\hat{\Phi}_{in}\hat{r}_{in}/3$ which arises due to the enclosed mass (gas and condensate) at the position $\hat{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 $\hat{\chi}_{\infty}$ is desired, rather than having $\hat{\Phi}_{in}$ as boundary condition, the condition $\hat{\chi} \equiv \hat{\chi}_{in}$ can first be set to that desired value $\hat{\chi}_{in} = \hat{\chi}_{\infty}$. After one numerical iteration the integration reaches some asymptotic value $\hat{\chi}_{in}^n$ (indexed by $n$) and the boundary condition is updated for the next iteration so that $\Delta\hat{\chi}(\hat{r}_{in}) = \hat{\chi}_{\infty} - \hat{\chi}_{in}^n$. This is repeated until the asymptotic value of $\hat{\chi}$ equals $\hat{\chi}_{\infty}$; we found for the cases studied that 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 $\hat{P}_{\Phi}$ 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 $\Pi \equiv d\hat{P}_{\Phi}/d\hat{r}$, The resulting equations in dimensionless form are

$$ \frac{d\hat{P}_{\Phi}}{d\hat{r}} = \Pi \qquad (B.5) $$

$$ \frac{d\Pi}{d\hat{r}} = \begin{cases} \frac{2}{\hat{r}}\Pi - \hat{\mu}^2 \left[ \hat{r}\text{Sign}(P_{\Phi})\sqrt{|\hat{P}_{\Phi}|} + \hat{P}_{\Phi} \right] + \hat{r}^2 \frac{d\hat{\rho}_b}{d\hat{r}} & \text{Vacuum or prescribed source} \\ \frac{2}{\hat{r}}\Pi - (\hat{\rho}_0 e^{\hat{\Phi}_0 - \hat{\Phi}} + \hat{\mu}^2) \left[ \hat{r}\text{Sign}(P_{\Phi})\sqrt{|\hat{P}_{\Phi}|} + \hat{P}_{\Phi} \right] & \text{Isothermal} \end{cases} \qquad (B.6) $$

where in the isothermal case, $\hat{\Phi}$ is obtained from (B.9) below. The equations were solved using an explicit Runge-Kutta method of order 8(5,3) (DOP853) [60] as implemented in SciPy [59]. We specified boundary conditions as in the two-field system above, and translated them into boundary conditions $\hat{P}_{\Phi_{in}}$ and $\Pi_{in}$ using

$$ \frac{d\hat{\Phi}}{d\hat{r}} = \frac{1}{\hat{r}} \left[ \text{Sign}(\hat{P}_{\Phi})\sqrt{|\hat{P}_{\Phi}|} + \frac{\hat{P}_{\Phi}}{\hat{r}} \right] \Leftrightarrow \sqrt{|\hat{P}_{\Phi}|} = \frac{\hat{r}}{2} \left( -1 + \sqrt{1 + 4 \left| \frac{d\hat{\Phi}}{d\hat{r}} \right|} \right) \qquad (B.7) $$

which is found after combining (3.8) and (3.11) for our choice of interpolation function and $\text{Sign}(\hat{P}_{\Phi}) = \text{Sign}(d\hat{\Phi}/d\hat{r})$. Likewise the correspondence between $\hat{\Phi}$ and $\Pi$ is given by

$$ \hat{\Phi} = W - \frac{\Pi}{\hat{\mu}^2 \hat{r}^2} \Leftrightarrow \Pi = \hat{r}^2 (V - \hat{\mu}^2 \hat{\Phi}) \qquad (B.8) $$

where $W/\hat{\mu}^2 = V = \hat{\rho}_b(\hat{r})$ in vacuum or for a prescribed source and

$$ W = W \left[ \hat{\rho}_0 e^{\frac{\Pi}{\hat{\mu}^2 \hat{r}^2} + \hat{\Phi}_0} \right] / \hat{\mu}^2 \qquad V = \hat{\rho}_0 e^{\hat{\Phi}_0 - \hat{\Phi}} \qquad (B.9) $$

for the isothermal sphere, where $W(x)$ is the Lambert (product logarithm) function.

We have compared the numerical solutions of the two-component system and the reduced Hamiltonian system for the vacuum case in figure 1. Here, figure 10 shows the excellent agreement between the two-component system and the reduced Hamiltonian system for the uniform density sphere and for the hydrostatic isothermal gas also.
<figure><img src="image_19.png" /><figcaption>**Figure 10.** Comparison between the numerical solutions for $|\Phi'|/a_0$ versus dimensionless radius $\hat{r}$ obtained with the two-component system (dashed-dotted orange line) and with the reduced Hamiltonian system (red line). We show the case of a prescribed source (uniform density sphere, left figure) and the case of the hydrostatic isothermal gas (right figure) demonstrating excellent agreement.</figcaption></figure>

### C Low-density RAR

In low-density systems it is possible for the AeST RAR to display a bifurcation, that is, have two possible values of $\Phi'$ for a given $\Phi'_N$ as seen in figure 11. Such a behaviour is not possible in MOND.

The explanation is as follows. A force may be low, and hence depart from the Newtonian expectation, for two reasons: either little mass has been enclosed, or the distance to the centre is large. In MOND, the two cases are on the same point on the RAR, but in AeST, the presence of the mass term $\mu$ distinguishes the two cases. It is helpful to think in terms of apparent $\chi$ condensates: one either encloses much $\chi$ condensate or little $\chi$ condensate, more $\chi$ condensate being enclosed when one is far away from the source.

### D Pure MOND isothermal sphere

#### 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 $\hat{P}'_\Phi \ge 0$, i.e. $\hat{P}_\Phi$ cannot decrease, and moreover integrating (3.35) leads to $GM(r) = P_\Phi(r)$, that is $P_\Phi \ge 0$ with $P_\Phi = 0$ only at $r = 0$. With these considerations, (3.8) implies that $\Phi' \ge 0$ (with strict equality only at $r = 0$), and thus (3.27) implies that $\rho' < 0$ which means that the density will always decrease as $r \to \infty$. Given that also $\rho \ge 0$, then $\rho$ must asymptote to a positive constant as $r \to \infty$, implying also that $\rho' \to 0$ as $r \to \infty$. Given that $\Phi' \ne 0$ except at $r = 0$, this then implies the stronger condition that $\rho \to 0$ as $r \to \infty$. From (3.35), we then get that $\hat{P}'_\Phi \to 0$ as $r \to \infty$, hence, $\hat{P}_\Phi$ is a monotonically increasing function which asymptotes to a constant as $r \to \infty$.
<figure><img src="image_20.png" /><figcaption>**Figure 11.** The RAR in AeST for the hydrostatic isothermal gas when densities are low, $\rho_0/\rho_I = 1$ (in low-acceleration regime throughout). When densities are low there can be two (bifurcation) accelerations $\Phi'$ in AeST (seen in both red and green lines) corresponding to a given Newtonian acceleration $\Phi'_N$, which happens 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.</figcaption></figure>

With the above considerations and since $\hat{r} > 0$ we may set $t = \ln \hat{r}$ so that (3.42) becomes an autonomous second order ODE

$$ \frac{d^2 \hat{P}_\Phi}{dt^2} + \left( \sqrt{\hat{P}_\Phi} - 3 \right) \frac{d \hat{P}_\Phi}{dt} = 0. \qquad (D.1) $$

Now set

$$ v = \frac{d\hat{P}_\Phi}{dt} \qquad (D.2) $$

so that

$$ \frac{dv}{dt} + \left( \sqrt{\hat{P}_{\Phi}} - 3 \right) v = 0 \qquad (D.3) $$

and if $v \neq 0$ we combine (D.2) and (D.3) to get

$$ \frac{dv}{d\hat{P}_{\Phi}} + \sqrt{\hat{P}_{\Phi}} - 3 = 0 \qquad (D.4) $$

with the solution

$$ v = 9v_0 + 3\hat{P}_{\Phi} - \frac{2}{3}\hat{P}_{\Phi}^{3/2} \qquad (D.5) $$

where $v_0$ is a constant. Thus, from (D.2) we find the first order separable ODE

$$ \frac{d\hat{P}_{\Phi}}{dt} = 9v_0 + 3\hat{P}_{\Phi} - \frac{2}{3}\hat{P}_{\Phi}^{3/2} \qquad (D.6) $$
which integrates to

$$t = t_0 + \int \frac{d\hat{P}_{\Phi}}{9v_0 + 3\hat{P}_{\Phi} - \frac{2}{3}\hat{P}_{\Phi}^{3/2}}. \qquad (D.7)$$

We require that $P_{\Phi} \ge 0$ and $P'_{\Phi} \ge 0$. If $v_0 = 0$, then $0 \le \hat{P}_{\Phi} \le \frac{81}{4}$ with $\hat{P}_{\Phi}(0) = 0$ and $\hat{P}_{\Phi}^{(\infty)} = \hat{P}_{\Phi}(\hat{r} \to \infty) = 81/4$. Increasing $v_0$ has no impact on the lower boundary while $\hat{P}_{\Phi}^{(\infty)}$ increases to higher values. However, a negative $v_0$ introduces a minimum $\hat{P}_{\Phi}$ which cannot happen at $\hat{r} = 0$. Thus for the $v_0 < 0$ to be physical, they must connect to a new regime (e.g., Newtonian) which should be valid at small radii down to $\hat{r} = 0$. For $v_0 \ge -1$, $P'_{\Phi} < 0$ except at a point, hence, there are no physical solutions of any kind in this case.

Now let us perform the integral. Since $\hat{P}_{\Phi} > 0$ by definition, we set $\hat{P}_{\Phi} = (q + 3/2)^2$, with $q \ge -3/2$, so that

$$t = t_0 - 3 \int \frac{q + \frac{3}{2}}{q^3 - \frac{27}{4}q + \tilde{v}_0} dq \qquad (D.8)$$

where $\tilde{v}_0 = -\frac{27}{4}(1 + 2v_0)$. We need to factorize the denominator so that we can expand the integrand in partial fractions and for this we need the roots of the cubic $q^3 - \frac{27}{4}q + \tilde{v}_0 = 0$. Clearly $\tilde{v}_0$, and thus $v_0$, play a fundamental role in how many real roots there can be. The discriminant of the cubic is $\Delta_q = -\frac{729}{4}v_0(1 + v_0)$ and if $v_0 > 0$ then there is only one real root, $q_1$, given by [61]

$$q_1 = 3 \cosh \left[ \frac{1}{3} \operatorname{arccosh} (1 + 2v_0) \right] \qquad (D.9)$$

and we have that $q_1 > 3$. For $v_0 = 0$ the roots are $q_I = \{3, -\frac{3}{2}, -\frac{3}{2}\}$ and for $-1 < v_0 < 0$ there are three real roots $q_I$.

## D.2 The case $v_0 = 0$

The case $v_0 = 0$ can be easily integrated as we have

$$t = t_0 - 3 \int \frac{1}{(q-3)(q+\frac{3}{2})} dq = t_0 - \frac{2}{3} \ln \frac{|q-3|}{|q+\frac{3}{2}|} \qquad (D.10)$$

leading after changing back to the original variables to

$$\hat{P}_{\Phi} = \frac{81}{4} \frac{\hat{r}^3 / \hat{r}_0^3}{\left(\sqrt{\hat{r}^3 / \hat{r}_0^3 + 1}\right)^2} \quad \forall \quad \hat{r} > 0 \qquad (D.11)$$

and

$$\hat{P}_{\Phi} = \frac{81}{4} \frac{\hat{r}^3 / \hat{r}_0^3}{\left(\sqrt{\hat{r}^3 / \hat{r}_0^3 - 1}\right)^2} \quad \forall \quad \hat{r} > \hat{r}_0. \qquad (D.12)$$

Solution (D.11) is a cored isothermal sphere valid for all $\hat{r} > 0$ while solution (D.12) is a singular isothermal sphere valid for $\hat{r} > \hat{r}_0$, with the singularity occurring exactly at the lower boundary.
### D.3 The case $v_0 > 0$

In the case $v_0 > 0$, $P_\Phi$ takes values in the range $0 \ge P_\Phi \ge \hat{P}_\Phi^{(\infty)}$ where $\hat{P}_\Phi^{(\infty)} > 81/4$. At $\hat{r} = 0$, the derivative $P'_\Phi > 0$ is strictly positive. We write the integral as

$$t = t_0 - 3 \int \frac{q + \frac{3}{2}}{(q - q_1) \left(q^2 + q_1 q + q_1^2 - \frac{27}{4}\right)} dq \quad (D.13)$$

$$= t_0 + \ln \left( \frac{\sqrt{q^2 + q_1 q + q_1^2 - \frac{27}{4}}}{|q - q_1|} \right)^{\frac{1}{q_1 - \frac{3}{2}}} - \frac{1}{q_1 - \frac{3}{2}} \sqrt{3 \frac{q_1 - 3}{q_1 + 3}} \arctan \left( \frac{2q + q_1}{\sqrt{3\sqrt{q_1^2 - 9}}} \right) \quad (D.14)$$

where we have used that $q^2 + q_1 q + q_1^2 - \frac{27}{4} \ge 0$ for $q \ge -3/2$, and where the second line follows after expanding the integrand in partial fractions, and integrating some of them to get the logarithms while completing a square in one term, leading to arctan. Back to the original variables we find a branch $0 \le P_\Phi \le P_\Phi^{(\infty)}$ (corresponding to $-3/2 \le q \le q_1$) with $P_\Phi^{(\infty)} = (q_1 + 3/2)^2$, which is the analogue of (D.11) and valid for $\hat{r} \ge 0$, and a branch $P_\Phi \ge P_\Phi^{(\infty)}$ which is the analogue of (D.12) and valid for $\hat{r} \ge \hat{r}_0$. Unfortunately, (D.14) cannot be inverted to find $P_\Phi(r)$ in terms of elementary functions, but we can consider the large $r$ expansion of $P_\Phi$. We find

$$\hat{P}_{\Phi} \approx \left(q_1 + \frac{3}{2}\right) \left\{ q_1 + \frac{3}{2} - 2 \left[ \left(\frac{81}{4}\right)^{1/3} \frac{\hat{r}_0}{\hat{r}} \right]^{q_1 - \frac{3}{2}} + \dots \right\}. \quad (D.15)$$

Considering (D.15) further we find the large $r$ limit of the density

$$4\pi G\rho \approx \frac{4}{9} \left(\frac{9}{2}\right)^{\frac{2q_1}{3}} \frac{a_0^2 \xi (q_1^2 - \frac{9}{4})}{\hat{r}_0^3} \left(\frac{\hat{r}_0}{a_0 \xi r}\right)^{q_1 + \frac{3}{2}} \quad (D.16)$$

potential

$$\hat{\Phi} \approx \hat{\Phi}_0 + \frac{1}{\xi} \ln \left[ \frac{9 \hat{\rho}_0 \hat{r}_0^3}{4 q_1^2 - 9} \left(\frac{2}{9}\right)^{\frac{2q_1}{3}} \right] + \frac{q_1 + \frac{3}{2}}{\xi} \ln \left(\frac{\hat{r}}{\hat{r}_0}\right) \quad (D.17)$$

and force

$$\frac{d\hat{\Phi}}{d\hat{r}} \approx \frac{q_1 + \frac{3}{2}}{\xi} \frac{1}{\hat{r}}. \quad (D.18)$$

In the limit $q_1 \to 3$ (i.e. $v_0 \to 0$), all approximations above agree with the expansion of (D.11) for small and large $r$ respectively.

### References

[1] V.C. Rubin and W.K. Ford Jr., *Rotation of the Andromeda nebula from a spectroscopic survey of emission regions*, Astrophys. J. 159 (1970) 379 [INSPIRE].

[2] V.C. Rubin, N. Thonnard and W.K. Ford Jr., *Rotational properties of 21 SC galaxies with a large range of luminosities and radii*, from NGC 4605 $R = 4$ kpc to UGC 2885 $R = 122$ kpc, Astrophys. J. 238 (1980) 471 [INSPIRE].

[3] A. Bosma, 21-cm line studies of spiral galaxies. 2. The distribution and kinematics of neutral hydrogen in spiral galaxies of various morphological types, Astron. J. 86 (1981) 1825 [INSPIRE].
[4] F. Zwicky, *Die Rotverschiebung von extragalaktischen Nebeln* (in German), *Helv. Phys. Acta* **6** (1933) 110 [INSPIRE].

[5] S.D.M. White, J.F. Navarro, A.E. Evrard and C.S. Frenk, *The baryon content of galaxy clusters: a challenge to cosmological orthodoxy*, *Nature* **366** (1993) 429 [INSPIRE].

[6] S.T. Myers et al., *Measurements of the Sunyaev-Zeldovich effect in the nearby clusters A478, A2142 and A2256*, *Astrophys. J.* **485** (1997) 1 [astro-ph/9703123] [INSPIRE].

[7] L. Grego et al., *Galaxy cluster gas mass fractions from Sunyaev-Zel'dovich effect measurements: constraints on $\Omega_M$, Astrophys. J.* **552** (2001) 2 [astro-ph/0012067] [INSPIRE].

[8] S.A. Grossman and R. Narayan, *Gravitationally lensed images in Abell 370*, *Astrophys. J.* **344** (1989) 637.

[9] 2DFGRS collaboration, *The 2dF galaxy redshift survey: the power spectrum and the matter content of the universe*, Mon. Not. Roy. Astron. Soc. **327** (2001) 1297 [astro-ph/0105252] [INSPIRE].

[10] M.M. Ivanov, M. Simonović and M. Zaldarriaga, *Cosmological parameters from the BOSS galaxy power spectrum*, JCAP **05** (2020) 042 [arXiv:1909.05277] [INSPIRE].

[11] BOOMERANG collaboration, *A flat universe from high resolution maps of the cosmic microwave background radiation*, *Nature* **404** (2000) 955 [astro-ph/0004404] [INSPIRE].

[12] WMAP collaboration, *First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: preliminary maps and basic results*, Astrophys. J. Suppl. **148** (2003) 1 [astro-ph/0302207] [INSPIRE].

[13] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. **641** (2020) A6 [Erratum ibid. **652** (2021) C4] [arXiv:1807.06209] [INSPIRE].

[14] PARTICLE DATA GROUP collaboration, *Review of particle physics*, PTEP **2022** (2022) 083C01 [INSPIRE].

[15] M. Milgrom, *A modification of the Newtonian dynamics: implications for galaxies*, Astrophys. J. **270** (1983) 371 [INSPIRE].

[16] M. Milgrom, *A modification of the Newtonian dynamics: implications for galaxy systems*, Astrophys. J. **270** (1983) 384 [INSPIRE].

[17] M. Milgrom, *A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis*, Astrophys. J. **270** (1983) 365 [INSPIRE].

[18] F. Lelli et al., *The baryonic Tully-Fisher relation for different velocity definitions and implications for galaxy angular momentum*, Mon. Not. Roy. Astron. Soc. **484** (2019) 3267 [arXiv:1901.05966] [INSPIRE].

[19] J. Bekenstein and M. Milgrom, *Does the missing mass problem signal the breakdown of Newtonian gravity?*, Astrophys. J. **286** (1984) 7 [INSPIRE].

[20] M. Milgrom, *Models of a modified-inertia formulation of MOND*, Phys. Rev. D **106** (2022) 064060 [arXiv:2208.07073] [INSPIRE].

[21] F. Lelli, S.S. McGaugh, J.M. Schombert and M.S. Pawlowski, *One law to rule them all: the radial acceleration relation of galaxies*, Astrophys. J. **836** (2017) 152 [arXiv:1610.08981] [INSPIRE].

[22] H. Desmond, D.J. Bartlett and P.G. Ferreira, *On the functional form of the radial acceleration relation*, Mon. Not. Roy. Astron. Soc. **521** (2023) 1817 [arXiv:2301.04368] [INSPIRE].
[23] H. Desmond, *The underlying radial acceleration relation*, Mon. Not. Roy. Astron. Soc. 526 (2023) 3342 [arXiv:2303.11314] [INSPIRE].

[24] M.M. Brouwer et al., *The weak lensing radial acceleration relation: constraining modified gravity and cold dark matter theories with KiDS-1000*, Astron. Astrophys. 650 (2021) A113 [arXiv:2106.11677] [INSPIRE].

[25] T. Mistele et al., *Radial acceleration relation of galaxies with joint kinematic and weak-lensing data*, arXiv:2310.15248 [INSPIRE].

[26] R. Stiskalek and H. Desmond, *On the fundamentality of the radial acceleration relation for late-type galaxy dynamics*, Mon. Not. Roy. Astron. Soc. 525 (2023) 6130 [arXiv:2305.19978] [INSPIRE].

[27] R.H. Sanders, *Clusters of galaxies with modified Newtonian dynamics (MOND)*, Mon. Not. Roy. Astron. Soc. 342 (2003) 901 [astro-ph/0212293] [INSPIRE].

[28] Y. Tian et al., *The radial acceleration relation in CLASH galaxy clusters*, Astrophys. J. 896 (2020) 70 [arXiv:2001.08340] [INSPIRE].

[29] P. Li et al., *Measuring galaxy cluster mass profiles into the low-acceleration regime with galaxy kinematics*, Astron. Astrophys. 677 (2023) A24 [arXiv:2303.10175] [INSPIRE].

[30] D. Eckert et al., *The gravitational field of X-COP galaxy clusters*, Astron. Astrophys. 662 (2022) A123 [arXiv:2205.01110] [INSPIRE].

[31] J.D. Bekenstein, *Phase coupling gravitation: symmetries and gauge fields*, Phys. Lett. B 202 (1988) 497 [INSPIRE].

[32] R.H. Sanders, *A stratified framework for scalar-tensor theories of modified dynamics*, Astrophys. J. 480 (1997) 492 [astro-ph/9612099] [INSPIRE].

[33] J.D. Bekenstein, *Relativistic gravitation theory for the MOND paradigm*, Phys. Rev. D 70 (2004) 083509 [Erratum ibid. 71 (2005) 069901] [astro-ph/0403694] [INSPIRE].

[34] R.H. Sanders, *A tensor-vector-scalar framework for modified dynamics and cosmic dark matter*, Mon. Not. Roy. Astron. Soc. 363 (2005) 459 [astro-ph/0502222] [INSPIRE].

[35] C. Skordis, *Generalizing tensor-vector-scalar cosmology*, Phys. Rev. D 77 (2008) 123502 [arXiv:0801.1985] [INSPIRE].

[36] M. Milgrom, *Bimetric MOND gravity*, Phys. Rev. D 80 (2009) 123536 [arXiv:0912.0790] [INSPIRE].

[37] L. Blanchet and S. Marsat, *Relativistic MOND theory based on the Khronon scalar field*, in the proceedings of the 10th International Conference on Gravitation, Astrophysics and Cosmology, (2012) [arXiv:1205.0400] [INSPIRE].

[38] B. Famaey and S. McGaugh, *Modified newtonian dynamics (MOND): observational phenomenology and relativistic extensions*, Living Rev. Rel. 15 (2012) 10 [arXiv:1112.3960] [INSPIRE].

[39] C. Skordis, D.F. Mota, P.G. Ferreira and C. Boehm, *Large scale structure in Bekenstein's theory of relativistic modified newtonian dynamics*, Phys. Rev. Lett. 96 (2006) 011301 [astro-ph/0505519] [INSPIRE].

[40] S. Hou and Y. Gong, *Gravitational waves in Einstein-Aether theory and generalized TeVeS theory after GW170817*, Universe 4 (2018) 84 [arXiv:1806.02564] [INSPIRE].

[41] C. Skordis and T. Złośnik, *Gravitational alternatives to dark matter with tensor mode speed equaling the speed of light*, Phys. Rev. D 100 (2019) 104013 [arXiv:1905.09465] [INSPIRE].
[42] C. Skordis and T. Zlosnik, New relativistic theory for modified newtonian dynamics, Phys. Rev. Lett. 127 (2021) 161302 [arXiv:2007.00082] [INSPIRE].

[43] C. Skordis and T. Zlosnik, *Aether scalar tensor theory: linear stability on Minkowski space*, Phys. Rev. D 106 (2022) 104041 [arXiv:2109.13287] [INSPIRE].

[44] R.C. Bernardo and C.-Y. Chen, *Dressed black holes in the new tensor-vector-scalar theory*, Gen. Rel. Grav. 55 (2023) 23 [arXiv:2202.08460] [INSPIRE].

[45] T. Kashfi and M. Roshan, *Cosmological dynamics of relativistic MOND*, JCAP 10 (2022) 029 [arXiv:2204.05672] [INSPIRE].

[46] T. Mistele, *Cherenkov radiation from stars constrains hybrid MOND dark matter models*, JCAP 11 (2022) 008 [arXiv:2103.16954] [INSPIRE].

[47] T. Mistele, S. McGaugh and S. Hossenfelder, *Æther scalar tensor theory confronted with weak lensing data at small accelerations*, Astron. Astrophys. 676 (2023) A100 [arXiv:2301.03499] [INSPIRE].

[48] S. Tian, S. Hou, S. Cao and Z.-H. Zhu, *Time evolution of the local gravitational parameters and gravitational wave polarizations in a relativistic MOND theory*, Phys. Rev. D 107 (2023) 044062 [arXiv:2302.13304] [INSPIRE].

[49] C. Llinares, *Extension of general relativity with MOND limit predicts novel orbital structure in and around galaxies*, arXiv:2302.12032 [INSPIRE].

[50] P. Verwayen, C. Skordis and C. Bœhm, *Æther Scalar Tensor (AeST) theory: quasistatic spherical solutions and their phenomenology*, arXiv:2304.05134 [INSPIRE].

[51] T. Mistele, *The two quasi-static limits of Æther scalar tensor theory*, arXiv:2305.07742 [INSPIRE].

[52] M. Bataki, C. Skordis and T. Zlosnik, *Æther scalar tensor theory: Hamiltonian formalism*, arXiv:2307.15126 [INSPIRE].

[53] R.J. Scherrer, *Purely kinetic k-essence as unified dark matter*, Phys. Rev. Lett. 93 (2004) 011301 [astro-ph/0402316] [INSPIRE].

[54] N. Arkani-Hamed, H.-C. Cheng, M.A. Luty and S. Mukohyama, *Ghost condensation and a consistent infrared modification of gravity*, JHEP 05 (2004) 074 [hep-th/0312099] [INSPIRE].

[55] N. Arkani-Hamed et al., *Dynamics of gravity in a Higgs phase*, JHEP 01 (2007) 036 [hep-ph/0507120] [INSPIRE].

[56] J. Binney and S. Tremaine, *Galactic dynamics*, second edition, (2008).

[57] J. Freundlich et al., *Probing the radial acceleration relation and the strong equivalence principle with the Coma cluster ultra-diffuse galaxies*, Astron. Astrophys. 658 (2022) A26 [arXiv:2109.04487] [INSPIRE].

[58] J.R. Dormand and P.J. Prince, *A family of embedded Runge-Kutta formulae*, J. Comput. Appl. Math. 6 (1980) 19.

[59] P. Virtanen et al., *SciPy 1.0 — fundamental algorithms for scientific computing in python*, Nature Meth. 17 (2020) 261 [arXiv:1907.10121] [INSPIRE].

[60] E. Hairer, S. Nørsett and G. Wanner, *Solving ordinary differential equations I: nonstiff problems*, Springer, Berlin, Heidelberg, Germany (1993) [DOI:10.1007/978-3-540-78862-1].

[61] G.C. Holmes, *The use of hyperbolic cosines in solving cubic polynomials*, Math. Gazette 86 (2002) 473.