© The Authors 2023

1. Introduction

Recently, various models have been proposed that combine the successes of modified Newtonian dynamics (MOND, Milgrom 1983a,b,c; Bekenstein & Milgrom 1984) on galactic scales with those of the Λ cold dark matter model (ΛCDM) on cosmological scales. Examples are superfluid dark matter (SFDM, Berezhiani & Khoury 2015; Berezhiani et al. 2018), the aether scalar tensor (AeST) model (Skordis & Złosnik 2021, 2022), and the neutrino-based model by Angus (2009; νHDM). The focus of the present paper is on the AeST model. An accompanying paper will look at SFDM.

The AeST model is a relativistic model that can reproduce MOND in the vicinity of galaxies and fits the fluctuations in the cosmic microwave background (CMB) as well as the matter power spectrum. In addition, it has a tensor mode that propagates at the speed of light which avoids difficulties matching the observations associated with GW170817 (Sanders 2018; Boran et al. 2018).

An important ingredient in the AeST model is a so-called ghost condensate (Arkani-Hamed et al. 2004, 2007). This ghost condensate is the major difference between the action of the AeST model in the static limit and the standard MOND-type action for multifield theories (Famaey & McGaugh 2012). The ghost condensate has an energy density that acts as an additional source for the gravitational field equations.

The integrated mass of the ghost condensate is generally negligible close to galaxies, where rotation curves are measured. Beyond a few hundred kiloparsecs, however, the ghost condensate mass is no longer negligible compared to the baryonic mass which leads to deviations from MOND.

This is a desired feature for galaxy clusters where observations require accelerations larger than what MOND predicts (Aguirre et al. 2001; Sanders 2003; Eckert et al. 2022). It may, however, be in conflict with unprecedented recent observations at large radii around galaxies: The analysis of weak-gravitational lensing data from Brouwer et al. (2021) found MOND-like behavior around galaxies up to ∼1 Mpc. Here, we explore whether this finding is compatible with the AeST model.

Since MOND is known to fit these observations well, we adopt an indirect approach to comparing the AeST model to observations. First, we introduce a method to quantify the deviation of the AeST model from MOND and analyze for which solutions these deviations are minimal. We then compare these optimal solutions – as well as slightly suboptimal ones – to the observational data.

2. Equations of motion and chemical potential

For galaxies we can use the quasi-static weak-field limit of the AeST model. In this limit, the model can be described by two fields, $\widehat{\mathrm{\Phi }}$ and φ, whose equations of motion are (at least in the spherically symmetric case that interests us here, see Appendix A),

$$\begin{array}{cc}\hfill \mathrm{\Delta }\widehat{\mathrm{\Phi }}& =f_G·4\pi G_{\mathrm{N}}({\rho }_{\mathrm{b}}+{\rho }_{\mathrm{c}}),\hfill \end{array}$$(1a)

$$\begin{array}{cc}\hfill \mathbf{\nabla }\left(\stackrel{\sim }{\mu }\left(\frac{|\mathbf{\nabla }\phi |}{a_0}\right)\mathbf{\nabla }\phi \right)& =f_G·4\pi G_{\mathrm{N}}({\rho }_{\mathrm{b}}+{\rho }_{\mathrm{c}}),\hfill \end{array}$$(1b)

where a₀ is the MOND acceleration scale and f_G is the conversion factor between Newton’s gravitational constant G_N and the constant $\widehat{G}$ that appears in the Lagrangian (see Appendix A). We use a₀ = 1.2 × 10⁻¹⁰ m s⁻² (Lelli et al. 2017). Both fields are sourced by the baryonic energy density ρ_b and the ghost condensate density ρ_c. The function $\stackrel{\sim }{\mu }$ can be freely chosen in this model and corresponds to an interpolation function of MOND. That is, it determines how the model interpolates between Newtonian gravity at large accelerations and MOND-like gravity at small accelerations. The constant f_G determines how much of the total acceleration (see below) in the Newtonian limit is due to φ and how much is due to $\widehat{\mathrm{\Phi }}$.

In the AeST model, ordinary matter couples to the metric g_μν in the usual way. The metric has the same form as in general relativity (GR), but with the Newtonian potential Φ being a combination of $\widehat{\mathrm{\Phi }}$ and φ, namely $\mathrm{\Phi }\equiv \widehat{\mathrm{\Phi }}+\phi$. Thus, the total acceleration felt by matter is ${\mathit{a}}_{\mathrm{tot}}\equiv {\mathit{a}}_{\widehat{\mathrm{\Phi }}}+{\mathit{a}}_{\phi }\equiv -\mathbf{\nabla }\widehat{\mathrm{\Phi }}-\mathbf{\nabla }\phi$.

The ghost condensate energy density ρ_c is given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{c}}=\frac{m^2}{4\pi G_{\mathrm{N}}{f}_G}(\frac{\dot{\phi }}{Q_0}-\widehat{\mathrm{\Phi }}-\phi ),\end{array}$$(2)

where m and Q₀ are constants. For cosmology, Skordis & Złosnik (2021) considered nonlinear corrections to this. We do not include these here for reasons discussed in Appendix I.

We did not set $\dot{\phi }=0$ in Eq. (2) because any constant $\dot{\phi }=-\dot{\widehat{\mathrm{\Phi }}}$ still gives time-independent equations of motion. Indeed, $\dot{\phi }$ represents the chemical potential of the condensate. To see this, one first observes that the model is shift-symmetric under $\phi \to \phi +\stackrel{\sim }{c}$, $\widehat{\mathrm{\Phi }}\to \widehat{\mathrm{\Phi }}-\stackrel{\sim }{c}$ for any constant $\stackrel{\sim }{c}$. In general, to describe equilibrium states, one introduces a chemical potential μ for each symmetry by shifting the Hamiltonian H by H → H − μ Q where Q is the conserved quantity associated with the symmetry. In the AeST model and on the level of the Lagrangian, this corresponds to shifting $\dot{\phi }\to \dot{\phi }+\mu$ and $\dot{\widehat{\mathrm{\Phi }}}\to \dot{\widehat{\mathrm{\Phi }}}-\mu$ (Mistele 2019; Kapusta 1981; Haber & Weldon 1982; Bilic 2008) or equivalently to considering solutions with $\dot{\phi }=\mu$ and $\dot{\widehat{\mathrm{\Phi }}}=-\mu$. (We note that the parameter m was called μ in Skordis & Złosnik 2021. We use μ instead to denote the chemical potential).

Consequently, the behavior of the AeST model around galaxies depends on the choice of this chemical potential of the ghost condensate. This corresponds to a choice of boundary condition for the combination μ/Q₀ − Φ which is a gauge-invariant variable as shown in Skordis & Złosnik (2022). This is an important difference to MOND where ρ_b alone determines the phenomenology around galaxies.

For real galaxies, these chemical potentials are ultimately set by galaxy formation. Since no simulations of nonlinear structure formation in the AeST model are available, we treated the chemical potential of each galaxy as a free parameter. In this regard, the AeST model is similar to SFDM. Both models require a choice of chemical potential to make predictions in galaxies (Berezhiani & Khoury 2015; Berezhiani et al. 2018; Mistele 2019, 2021; Hossenfelder & Mistele 2020; Mistele et al. 2022).

We now assume spherical symmetry. Then, solutions have the same form as in MOND, just with the baryonic mass replaced by an effective mass that includes the condensate mass,

$$\begin{array}{c}\hfill M_{\mathrm{eff}}(r)\equiv M_{\mathrm{b}}(r)+M_{\mathrm{c}}(r)\equiv 4\pi {\displaystyle {\int }_0^r\mathrm{d}r\prime r{\prime }^2({\rho }_{\mathrm{b}}(r\prime )+{\rho }_{\mathrm{c}}(r\prime )).}\end{array}$$(3)

Thus, with the Newtonian acceleration a_N in the negative radial direction,

$$\begin{array}{c}\hfill a_{\mathrm{N}}(r)\equiv a_{\mathrm{b}}(r)+a_{\mathrm{c}}(r)\equiv \frac{G_{\mathrm{N}}{M}_{\mathrm{b}}(r)}{r^2}+\frac{G_{\mathrm{N}}{M}_{\mathrm{c}}(r)}{r^2},\end{array}$$(4)

we can write the total acceleration in the negative radial direction a_tot as a_tot = a_N(r)⋅ν(|a_N(r)|/a₀) where the interpolation function ν is determined by the free function $\stackrel{\sim }{\mu }$ in the AeST model (see Eq. (1)). Unlike in MOND, a_N can be negative because the condensate mass M_c can be negative. We discuss this in more detail below.

For simplicity, we chose the interpolation function such that

$$\begin{array}{c}\hfill a_{\mathrm{tot}}(r)=s_{\mathrm{eff}}·\left(|a_{\mathrm{N}}(r)|+\sqrt{a_0|a_{\mathrm{N}}(r)|}\right),\end{array}$$(5)

where s_eff is the sign of M_eff. We recover MOND by leaving out the contributions from the condensate, M_c = 0. Thus, when comparing the AeST model to MOND below, we assumed the following acceleration for MOND,

$$\begin{array}{c}\hfill a_{\mathrm{MOND}}(r)\equiv a_{\mathrm{b}}(r)+\sqrt{a_0{a}_{\mathrm{b}}(r)}.\end{array}$$(6)

This assumes the same interpolation function for MOND and the AeST model.

This interpolation function has the correct limits, namely a_MOND → a_b for accelerations much larger than a₀ and $a_{\mathrm{MOND}}\to \sqrt{a_0{a}_{\mathrm{b}}}$ for accelerations much smaller than a₀. Still, in general, this choice is too simplistic. But it suffices for our purposes because we are only interested in the small-acceleration regime where all interpolation functions give $a_{\mathrm{MOND}}\approx \sqrt{a_0{a}_{\mathrm{b}}}$.

The total acceleration a_tot in the AeST model does not explicitly depend on f_G and is the sum of $a_{\widehat{\mathrm{\Phi }}}=f_G{a}_{\mathrm{N}}$ and $a_{\phi }=s_{\text{eff}}((1-f_G)|a_{\mathrm{N}}|+\sqrt{a_0|a_{\mathrm{N}}|})$. Below, we refer to its MOND-like part $\sqrt{a_0|a_{\mathrm{N}}|}$ and its Newton-like part a_N, respectively, as

$$\begin{array}{c}\hfill a_{\stackrel{\sim }{\mathrm{\Phi }}}\equiv s_{\mathrm{eff}}|a_{\mathrm{N}}|,a_{\stackrel{\sim }{\phi }}\equiv s_{\mathrm{eff}}\sqrt{a_0|a_{\mathrm{N}}|}.\end{array}$$(7)

In the following we assumed point particle baryonic masses for simplicity. This suffices for our purposes, because the details of the baryonic mass distribution do not matter much at the large galactocentric radii we consider. See Appendix D for an approximate analytical solution for this case and Appendix E for how we calculate numerical solutions.

3. Deviations from MOND

The AeST model reproduces MOND so long as the condensate’s total mass is small compared to the baryonic mass. This condition is usually fulfilled in the inner parts of galaxies but not far away from the galaxy. Indeed, given a maximum allowed fractional deviation δ from MOND, there is an optimal boundary condition for which the MOND-like behavior extends to a finite maximum radius r_max. For all other boundary conditions, deviations from MOND set in earlier. This is illustrated in Fig. 1.

Specifically, we imposed the maximum allowed deviation δ as

$$\begin{array}{c}\hfill |\frac{a_{\stackrel{\sim }{\phi }}(r)}{\sqrt{a_0{a}_{\mathrm{b}}(r)}}-1|=|\frac{s_{\mathrm{eff}}·\sqrt{a_0|a_{\mathrm{N}}(r)|}}{\sqrt{a_0{a}_{\mathrm{b}}(r)}}-1|\le \delta .\end{array}$$(8)

That is, we compared the accelerations in the AeST model and MOND, focusing on the MOND-like contributions $\sqrt{a_0{a}_{\mathrm{N}}}$ and $\sqrt{a_0{a}_{\mathrm{b}}}$. Alternatively, one could compare the total acceleration in both models, that is, $a_{\mathrm{N}}+\sqrt{a_0{a}_{\mathrm{N}}}$ and $a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}}$. But at the large radii we consider here, the difference is negligible (see Appendix G) and we used Eq. (8) for simplicity.

The maximal radius r_max, up to which we allowed accelerations to deviate by less than a fraction δ from MOND, is given by

$$\begin{array}{c}\hfill \frac{r_{\max }}{r_{\mathrm{MOND}}}\approx 1.53{\left(9\frac{{(1+\delta )}^2-1}{r_{\mathrm{MOND}}^2{m}^2/f_G}\right)}^{1/3}.\end{array}$$(9)

Here, r_MOND is a constant known as the MOND-radius $\sqrt{G_{\mathrm{N}}{M}_{\mathrm{b}}/a_0}$. See Appendix G for a derivation of those properties and an analysis of the conditions under which this estimate holds.

This maximum radius r_max scales as $M_{\mathrm{b}}^{1/6}$. This means that more massive galaxies can stay close to MOND up to larger radii than less massive galaxies. In MOND, however, accelerations a_b = G_NM_b/r² are more relevant than radii r. For example, the Radial Acceleration Relation (RAR, Lelli et al. 2017) relates the total acceleration a_tot and the Newtonian baryonic acceleration a_b.

The maximum radius r_max corresponds to a minimum acceleration $a_{\text{b},\text{min}}=G_{\text{N}}{M}_{\text{b}}/r_{\text{max}}^2$ which scales as $M_{\mathrm{b}}^{2/3}$. Thus, in acceleration space, less massive galaxies can stay close to MOND for longer than more massive galaxies. We show a_b, min in Fig. 2 and illustrate the RAR for various baryonic masses in Fig. 3.

The scale of r_max is set by the combination m²/f_G. Skordis & Złosnik (2021) require m²/f_G ≲ 1 Mpc⁻². As we show below, this does not guarantee MOND-like behavior for weak lensing which probes radii up to ∼1 Mpc. One might therefore want to choose an even smaller m²/f_G. But this is not easily possible. Indeed, galaxy clusters require more acceleration than MOND predicts (Aguirre et al. 2001; Sanders 2003; Eckert et al. 2022). To naturally explain this in the AeST model, m²/f_G cannot be much smaller than 1 Mpc⁻² and we assumed

$$\begin{array}{c}\hfill m^2/f_G\sim {(1\mathrm{Mpc})}^{-2}.\end{array}$$(10)

We obtained this estimate by requiring that the minimum acceleration a_b, min in AeST matches where observed clusters deviate from MOND. For example, if we require that the total acceleration a_tot in a cluster with M_b = 10¹⁴ M_⊙ deviates by at least δ = 10% from MOND at a_b, min = 10^−10.5 m s⁻², we find m²/f_G > 2.5 Mpc⁻² (see also Table 1)¹.

Beyond the radius where the AeST model deviates from MOND, the gravitational force eventually becomes oscillatory (Skordis & Złosnik 2021), see also Appendix F, as is typical for condensate models (Arkani-Hamed et al. 2007). The reason is that M_eff(r) and the condensate density ρ_c oscillate. However, condensates with negative energy-density are unstable (or at least we would expect a good reason for why they are not unstable). We therefore expect the AeST model to be unstable in this oscillatory regime. We would then no longer have a macroscopically coherent condensate and the quasi-static action Eq. (A.1) is no longer good to use.

In superfluid dark matter models, for example, we assume that when the energy-density begins to oscillate that we have to continue the condensate density by a standard (nonsuperfluid) phase (Berezhiani & Khoury 2015; Berezhiani et al. 2018). Something similar might happen in the AeST model.

However, so far, a stability analysis for the AeST model in a galactic background has not been done, so maybe the oscillatory regime turns out to be stable after all. Below, we therefore keep in mind both possibilities and discuss where our results depend on whether or not negative condensate densities are stable. For example, Fig. 4 shows the solutions from Fig. 3 but truncated where the condensate density first drops to zero.

4. Weak lensing

The AeST model usually reproduces MOND at the radii probed by rotation curves because the condensate density M_c is negligible there. This means that probing the effects of the condensate requires a different approach. The option we pursued here is to use the recent weak-lensing analysis of Brouwer et al. (2021) who find that accelerations are MOND-like down to at least a_b ∼ 10⁻¹³ m s⁻².

In the AeST model, matter is coupled to the fields φ and $\widehat{\mathrm{\Phi }}$ through the metric g_μν. In the weak-field limit, this metric has the same form as in GR, just with the Newtonian potential replaced by $\mathrm{\Phi }=\widehat{\mathrm{\Phi }}+\phi$. We can therefore use the standard formalism for weak lensing just by taking into account $\mathrm{\Phi }=\widehat{\mathrm{\Phi }}+\phi$.

Most galaxies in the weak-lensing sample from Brouwer et al. (2021) have baryonic masses between 10¹⁰ M_⊙ and 10¹¹ M_⊙. For the AeST model with f_G/m² = 0.99 Mpc², Fig. 2 shows that MOND-like behavior up to a few 10% is possible down to a_b ∼ 10⁻¹³ m s⁻² for those galaxies in the sample with baryonic masses ∼10¹⁰ M_⊙.

Such 𝒪(10%) deviations may be sufficient to match observations. But it is important to keep in mind that in this regime where deviations from MOND start to become important, the details depend on the precise baryonic masses and boundary conditions of the galaxy sample as well as the precise value of the model parameter m²/f_G. For example, Fig. 2 shows solutions for boundary conditions that are optimal for reproducing MOND. But there is no reason why galaxy formation should result in such optimal boundary conditions. One would therefore expect deviations of AeST from MOND to generally be larger than in the optimal case we depict.

We can derive a rough upper bound on the model parameter m²/f_G by requiring that AeST can reproduce MOND in the regime probed by weak lensing. For this, we used the minimum acceleration $a_{\text{b},\text{min}}=G_{\text{N}}{M}_{\text{b}}/r_{\text{max}}^2$ with r_max from Eq. (9). Many of the galaxies in the sample used by Brouwer et al. (2021) have baryonic masses close to 10¹¹ M_⊙. If we require that such galaxies can reproduce MOND down to a_b, min = 10⁻¹⁵ m s⁻² and up to a fraction δ = 10%, we find m²/f_G < 0.001 Mpc⁻². There is a lot of uncertainty in this upper bound. For example, m²/f_G = 0.001 Mpc⁻² is small enough that the ghost condensate density in galaxies as given by Eq. (2) is typically smaller than the cosmological background density. So, at least in principle, there could be corrections from the fact that we should expand around a cosmological background, not around empty Minkowski space (which is how Eq. (1) was derived)². Also, if we disregard the data below a_b = 10⁻¹³ m s⁻² because the isolation criterion used in the weak-lensing analysis is less reliable there (Brouwer et al. 2021), we obtain the weaker bound m²/f_G < 1 Mpc⁻². Still, a tension between the value of m²/f_G required by weak lensing and that required by galaxy clusters, m²/f_G ≳ 1 Mpc⁻², seems likely (see Table 1).

Another aspect to take into account is that, in practice, the weak-lensing RAR is not known for individual galaxies. It is known only in an averaged sense for a large sample of stacked galaxies. It is possible that stacking gives a MOND-like RAR even if most galaxies individually do not. As we show in Appendix H, in our case, stacking simply means calculating a weighted average in acceleration space. So, indeed, accelerations larger than MOND from some galaxies can cancel accelerations smaller than MOND from other galaxies.

To illustrate this, we considered a sample of galaxies with a fixed baryonic mass M_b but various boundary conditions. For simplicity, we weighed all galaxies equally when averaging. One example is shown in Fig. 5. We see how a stacked RAR can be MOND-like even when most individual stacked galaxies are not. Of course, one does not always get a MOND-like RAR from stacking. This works only when accelerations below and above the MOND prediction cancel each other. Whether or not it works for real lensing galaxies depends on which boundary conditions are picked by galaxy formation.

In addition, stacking galaxies with different boundary conditions should lead to increased uncertainties at small a_b where different boundary conditions lead to significantly different accelerations. In principle, such uncertainties might be visible in the error bars of the observed weak-lensing RAR.

However, we expect that in practice there is often no visible effect on the error bars. To see this, we first note that the error bars shown in Figs. 3 and 4 correspond to the uncertainty of the mean value of a_tot, obtained by stacking a large number N of galaxies. They do not represent the scatter from galaxy-by-galaxy variation. The galaxy-by-galaxy variation is larger than the uncertainty of the mean by a factor of about $\sqrt{N}$. In AeST, different boundary conditions induce a form of galaxy-by-galaxy variation. Thus, the scatter from solutions with different boundary conditions should be compared to the larger galaxy-by-galaxy variation and not to the uncertainty of the mean. Put differently, the scatter induced by different boundary conditions should be scaled by a factor $1/\sqrt{N}$ before comparing to the error bars shown in Figs. 3 and 4. Here, N is on the order of 10⁵ (Brouwer et al. 2021). Thus, one would expect a visible effect on the error bars only in extreme cases.

4.1. Negative condensate densities

Whether or not stacked galaxies can follow a MOND RAR down to smaller accelerations than individual galaxies also depends on whether or not negative condensate densities are stable. To see this, we note that accelerations are smaller than in MOND if and only if the effective mass M_eff = M_b + M_c is smaller than in MOND. This is only possible for negative M_c which requires negative condensate densities. Therefore, if negative densities are unstable, cancelling accelerations that are larger against those that are smaller than in MOND does not work. This is simply because all accelerations are larger than in MOND if we do not allow for negative densities. This is illustrated in Fig. 6.

Thus, if negative densities are unstable, one might expect that the AeST model always gives larger accelerations than MOND. Moreover, our estimate for a_b, min (see Fig. 2) suggests that these deviations should set in earlier for larger baryonic masses M_b. And indeed, there are hints of such behavior in the observed weak-lensing data, see Fig. 7 which shows the weak-lensing data separately for early-type galaxies (ETGs) and late-type galaxies (LTGs). We see that the weak-lensing RAR for LTGs follows the MOND prediction even for a_b < 10⁻¹³ m s⁻² while ETGs tend toward larger accelerations than MOND. In general, ETGs have larger baryonic masses than LTGs. So this seems to fit with the AeST model expectations if negative densities are unstable.

However, this M_b-dependence is not a plausible explanation for the difference between the observed weak-lensing RARs for ETGs and LTGs. This is for three reasons.

First, the ETGs and LTGs do not sufficiently differ in baryonic mass. What would be required is a difference in M_b of more than a factor 10^3/2. To see this, we note that LTGs follow the MOND prediction for at least one more order of magnitude in a_b compared to ETGs. This translates into a factor > 10^3/2 in terms of baryonic mass according to our estimate $a_{b,\min }\propto M_{\mathrm{b}}^{2/3}$. In contrast, Brouwer et al. (2021) selected LTGs and ETGs to have the same stellar mass distribution. With the stellar M/L scale corrected to be consistent with that of the observed kinematic RAR, ETG stellar masses would still be larger by a factor 1.4 (McGaugh, priv. comm.). But this is not sufficient here.

Of course, the total baryonic mass should take into account gas, but this is unlikely to account for the required factor of 10^3/2 or more, at least with the simple cold gas mass estimate used in Brouwer et al. (2021). Brouwer et al. (2021) also consider a scenario where ETGs have significantly more hot gas than LTGs. But even in that scenario, the baryonic masses of ETGs differ from those of LTGs by only a factor of two or a bit more. In addition, adopting this scenario means adopting a different observed lensing RAR. Indeed, as discussed in Brouwer et al. (2021), in this scenario there might not even be a discrepancy between ETGs and LTGs left to explain.

Second, the isolation criterion needed to obtain the weak-lensing RAR might fail for ETGs sooner than it does for LTGs. One might naturally expect this to be the case since ETGs are known to be more clustered than LTGs (Dressler 1980). At what point this comes into play here we cannot judge, but mention it as a logical possibility.

Third and finally, even if negative densities are indeed unstable, the AeST model does not necessarily predict larger accelerations than MOND. Indeed, it makes no physical sense to stop looking when the condensate becomes unstable. In a real galaxy, something else must follow after the condensate phase. For example, the macroscopically coherent ghost condensate might be replaced by something closer to a ΛCDM-like collisionless fluid which the AeST model postulates on cosmological scales (Skordis & Złosnik 2021). In principle, whatever replaces the ghost condensate might lead to smaller accelerations than MOND.

That said, the prediction of larger-than-MOND accelerations remains valid if whatever replaces the ghost condensate has as its only effect to replace the ghost condensate density by some other positive density, ρ_c → ρ_replace. This is because then one still has solutions of the same form as before, just with a different effective mass M_eff  →  M_b + M_replace > M_b.

In order to get smaller-than-MOND accelerations, the general structure of the solutions must be modified. That is, the left-hand sides of the equations of motions must be modified, $\mathrm{\Delta }\widehat{\mathrm{\Phi }}=\cdots ,\mathbf{\nabla }(\stackrel{\sim }{\mu }(|\mathbf{\nabla }\phi |/a_0)\mathbf{\nabla }\phi )=\cdots \to ?$. It is not implausible that this indeed happens since the field φ plays a role for both the ghost condensate (for example, it carries the chemical potential $\dot{\phi }=\mu$) as well as the gravitational force. So outside the condensate phase both could be modified.

4.2. External field effect

Another concern at the small accelerations probed by weak lensing is the external field effect (EFE) of MOND (Milgrom 1983b; Famaey & McGaugh 2012). The EFE is a consequence of the specific nonlinear form ∇(|∇φ|∇φ)∝ρ_b of the gravitational field equations in the small-acceleration limit. The crucial nonlinearity is the same in the AeST model, so there is probably a similar effect there, at least within the ghost condensate³.

The EFE generally reduces the observationally inferred a_tot, while a condensate mass M_c > 0 in the AeST model enhances this acceleration. In principle, these two effects could cancel each other to give a MOND-like acceleration even at very large galactocentric distances, but there is no reason to expect such a cancellation to generally happen.

And in any case, if negative densities are unstable and the condensate is replaced by something else at large radii, then any potential EFE depends on the details of what replaces the ghost condensate. There could be a modified nonlinear effect that still allows neighboring galaxies to affect each other in a way that violates the Strong Equivalence Principle (SEP) like the MOND EFE does. Or there could be no such effect, possibly restoring the SEP at large scales.

A restored SEP would fit with the fact that the observed weak-lensing RAR from Brouwer et al. (2021) shows no signs of an EFE. But here one must be careful. The EFE pertains to nonisolated galaxies, while the analysis of Brouwer et al. (2021) requires isolated galaxies. Thus, any EFE effects may be masked by violations of this assumption. In addition, the environment-dependence of the EFE is quite complicated (Llinares et al. 2008; Chae et al. 2021). So it is not even clear for MOND whether or not a significant EFE is expected here. Still, the EFE is something to keep in mind as observations and theoretical predictions are improved.

5. Conclusion

We have explored whether or not the AeST model can explain the observed MOND-like weak-lensing RAR which probes unprecedentedly small accelerations. We find that deviations from MOND start to set in already in the range of the new measurements, creating a tension with data. The model parameter m²/f_G can be adjusted to avoid this tension, but that likely creates a tension with observations of galaxy clusters instead.

It seems that keeping the model in agreement with data would require specific values of boundary conditions for a large variety of galaxies. While this is possible, we do not know of any mechanism that would result in these particular boundary conditions. Thus, while we cannot rule out the model, it does seem that weak-lensing observations pose a challenge for AeST.

---

Acknowledgments

We thank Margot Brouwer, Andrej Dvornik, and Kyle Oman for helpful correspondence. This work was supported by the DFG (German Research Foundation) under grant number HO 2601/8-1 together with the joint NSF grant PHY-1911909. This work was supported by the DFG (German Research Foundation) – 514562826.

References

Appendix A: Action and equations of motion

For galaxies, the quasi-static weak-field limit of the AeST model is relevant. The action in this limit is (Skordis & Złosnik 2021)

$$\begin{array}{cc}\hfill S=-{\displaystyle \int }{d}^{4}x& \{\frac{1}{8\pi \widehat{G}}[{(\mathbf{\nabla }\mathrm{\Phi })}^2-2\mathbf{\nabla }\mathrm{\Phi }\mathbf{\nabla }\phi +{(\mathbf{\nabla }\phi )}^2\hfill \\ \hfill & -m^2{(\frac{\dot{\phi }}{Q_0}-\mathrm{\Phi })}^2+\mathcal{J}\left({(\mathbf{\nabla }\phi )}^2\right)]+\mathrm{\Phi }{\rho }_{\mathrm{b}}\},\hfill \end{array}$$(A.1)

where $\widehat{G}$ is a constant and the function 𝒥 determines the MOND interpolation function. We discuss potential higher-order corrections to the m² term which produces the ghost condensate density in Appendix I.

The AeST model also contains a unit vector field A^μ. The form Eq. (A.1) of the action assumes A = 0 (Skordis & Złosnik 2021). In Appendix B we explain why this assumption is correct in spherical symmetry – which is what we are mainly interested in here – but not in general.

The quasi-static weak-field limit equations of motion derived from the action Eq. (A.1) are

$$\begin{array}{cc}\hfill \mathrm{\Delta }\widehat{\mathrm{\Phi }}& =4\pi \widehat{G}{\rho }_{\mathrm{b}}+m^2(\mu /Q_0-\widehat{\mathrm{\Phi }}-\phi ),\hfill \end{array}$$(A.2a)

$$\begin{array}{cc}\hfill \mathbf{\nabla }\left(\stackrel{\sim }{\mu }\left(\frac{|\mathbf{\nabla }\phi |}{a_0}\right)\mathbf{\nabla }\phi \right)& =4\pi \widehat{G}{\rho }_{\mathrm{b}}+m^2(\mu /Q_0-\widehat{\mathrm{\Phi }}-\phi ),\hfill \end{array}$$(A.2b)

where, following Skordis & Złosnik (2021), we have introduced $\widehat{\mathrm{\Phi }}$ through $\mathrm{\Phi }\equiv \widehat{\mathrm{\Phi }}+\phi$ and $\stackrel{\sim }{\mu }(|\mathbf{\nabla }\phi |/a_0)=\mathcal{J}\prime ({(\mathbf{\nabla }\phi )}^2)$. In order to reproduce MOND-like behavior for small accelerations and Newton-like behavior for large accelerations, the function $\stackrel{\sim }{\mu }(|\mathbf{\nabla }\phi |/a_0)$ must be proportional to |∇φ| for small arguments and it must be a constant for large arguments. One can parametrize these limits by a parameter λ_s (Skordis & Złosnik 2021),

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mu }\left(\frac{|\mathbf{\nabla }\phi |}{a_0}\right)|}_{|\mathbf{\nabla }\phi |\to 0}=\frac{{\lambda }_s}{1+{\lambda }_s}\frac{|\mathbf{\nabla }\phi |}{a_0},{\stackrel{\sim }{\mu }\left(\frac{|\mathbf{\nabla }\phi |}{a_0}\right)|}_{|\mathbf{\nabla }\phi |\to \infty }={\lambda }_s.\end{array}$$(A.3)

This ensures both a standard Newton regime at large accelerations and a standard MOND regime at small accelerations with Newtonian gravitational constant

$$\begin{array}{c}\hfill G_N\equiv \widehat{G}\frac{1+{\lambda }_s}{{\lambda }_s}\equiv \widehat{G}·f_G^{-1}.\end{array}$$(A.4)

How exactly the function $\stackrel{\sim }{\mu }$ interpolates between these two limits is not specified by the AeST model. Various choices are possible. A choice of $\stackrel{\sim }{\mu }$ corresponds to a choice of the so-called interpolation function in MOND (Famaey & McGaugh 2012). Up to the condensate density ρ_c, these equations are typical of multifield MOND models (Famaey & McGaugh 2012).

In the small-acceleration limit |∇φ|≪a₀ and without any baryonic density, ρ_b = 0, these equations approximately become ∇(|∇φ|∇φ) = c₁ + c₂φ with constants c₁ and c₂. This is a special case of the deep-MOND polytropes studied in Milgrom (2021). The results of Milgrom (2021) are not directly useful here, however, because the baryonic density ρ_b plays a very important role for the phenomenology of the AeST model as we see below.

Appendix B: The quasi-static limit more generally

Had we not set A to zero by hand, the action in the quasi-static limit Eq. (A.1) would read,

$$\begin{array}{cc}\hfill S& =-{\displaystyle \int }{d}^{4}x\{\frac{1}{8\pi \widehat{G}}[{(\mathbf{\nabla }\mathrm{\Phi })}^2-2\mathbf{\nabla }\mathrm{\Phi }(\mathbf{\nabla }\phi +Q_0\mathit{A})+{(\mathbf{\nabla }\phi +Q_0\mathit{A})}^2\hfill \\ \hfill & -m^2{(\frac{\dot{\phi }}{Q_0}-\mathrm{\Phi })}^2+\mathcal{J}\left({(\mathbf{\nabla }\phi +Q_0\mathit{A})}^2\right)+\frac{2K_B}{2-K_B}{\mathbf{\nabla }}_{[i}{\mathit{A}}_{j]}{\mathbf{\nabla }}^{[i}{\mathit{A}}^{j]}]+\mathrm{\Phi }{\rho }_{\mathrm{b}}\}.\hfill \end{array}$$(B.1)

We can decompose A into a divergence-less and a curl-less part, A ≡ ∇ × A_c + ∇A_d. For time-independent fields – that is, time-independent up to the chemical potential $\dot{\phi }=\mu =\mathrm{const}.$ – the results of Skordis & Złosnik (2022) show that we can set A_d = 0 by a gauge transformation.

In spherical symmetry, the curl term ∇ × A_c vanishes and the A equation of motion and the φ equation of motion are equivalent. Thus, setting A to zero and using Eq. (A.1) is justified. In general, however, setting A to zero is inconsistent. This can be seen from the A equation of motion,

$$\begin{array}{cc}& \mathbf{\nabla }\mathrm{\Phi }+\frac{1}{2Q_0}\frac{2K_B}{2-K_B}(\mathrm{\Delta }\mathit{A}-\mathbf{\nabla }(\mathbf{\nabla }·\mathit{A}))=\hfill \\ \hfill & (\mathbf{\nabla }\phi +Q_0\mathit{A})(1+{\mathcal{J}}^{\prime }\left({(\mathbf{\nabla }\phi +Q_0\mathit{A})}^2\right)).\hfill \end{array}$$(B.2)

If A were zero, we could infer

$$\begin{array}{c}\hfill \mathbf{\nabla }|\mathbf{\nabla }\mathrm{\Phi }|\times \mathbf{\nabla }\mathrm{\Phi }=0,\end{array}$$(B.3)

by algebraically solving for ∇φ and then taking the curl. But Eq. (B.3) holds only in very special situations, see for example Brada & Milgrom (1995). Thus, except in a few special cases, we must not set A = 0 in the quasi-static limit.

Of course, even when Eq. (B.3) does not hold, setting A to zero might still be a reasonable approximation akin to how it is often reasonable to neglect a curl term in standard models of MOND (Famaey & McGaugh 2012). Investigating this is left for future work.

Appendix C: General structure of the solutions

We now assume spherical symmetry. Then, solutions have the same form as in standard multifield MOND models except that the baryonic mass M_b is replaced by the effective mass M_eff which includes the condensate mass M_c in addition to M_b. The equations of motion are then solved by

$$\begin{array}{cc}\hfill a_{\widehat{\mathrm{\Phi }}}\equiv {\widehat{\mathrm{\Phi }}}^{\prime }(r)& =f_G{a}_N(r),\hfill \end{array}$$(C.1a)

$$\begin{array}{cc}\hfill a_{\phi }\equiv {\phi }^{\prime }(r)& =a_N(r)\stackrel{\sim }{\nu }\left(\frac{|a_N(r)|}{a_0}\right),\hfill \end{array}$$(C.1b)

where a_N = G_NM_eff(r)/r² and the function $\stackrel{\sim }{\nu }$ is determined by $\stackrel{\sim }{\mu }$. See, for example, Famaey & McGaugh (2012) for how these two are related. The total acceleration in the negative radial direction felt by matter is then

$$\begin{array}{cc}\hfill a_{\mathrm{tot}}(r)& ={\widehat{\mathrm{\Phi }}}^{\prime }(r)+{\phi }^{\prime }(r)=a_N(r)·(f_G+\stackrel{\sim }{\nu }\left(\frac{|a_N(r)|}{a_0}\right))\hfill \\ \hfill & \equiv a_N(r)·\nu \left(\frac{|a_N(r)|}{a_0}\right),\hfill \end{array}$$(C.2)

where ν is a MOND interpolation function (Famaey & McGaugh 2012).

Below, we are mainly interested in the deep-MOND regime a_N ≪ a₀. In this regime, the total acceleration does not depend much on the choice of interpolation function. Thus, we chose an interpolation function that is easy to handle,

$$\begin{array}{c}\hfill \nu (y)\equiv 1+1/\sqrt{y}.\end{array}$$(C.3)

This interpolation function is not, in general, suited to fit galaxies (Famaey & McGaugh 2012). But it is sufficient in the small-acceleration limit in which we are interested here. This implies for the accelerations due to the field $\widehat{\mathrm{\Phi }}$ and φ,

$$\begin{array}{cc}\hfill {\widehat{\mathrm{\Phi }}}^{\prime }(r)& =s_{\mathrm{eff}}·f_G|a_N(r)|,\hfill \end{array}$$(C.4a)

$$\begin{array}{cc}\hfill {\phi }^{\prime }(r)& =s_{\mathrm{eff}}·((1-f_G)|a_N(r)|+\sqrt{a_0|a_N(r)|}),\hfill \end{array}$$(C.4b)

and therefore for the total acceleration

$$\begin{array}{cc}\hfill {\widehat{\mathrm{\Phi }}}^{\prime }(r)+{\phi }^{\prime }(r)& =s_{\mathrm{eff}}·\left(|a_N(r)|+\sqrt{a_0|a_N(r)|}\right)\hfill \\ \hfill & =s_{\mathrm{eff}}·(\frac{G_N|M_{\mathrm{eff}}(r)|}{r^2}+\frac{\sqrt{a_0{G}_N|M_{\mathrm{eff}}(r)|}}{r}).\hfill \end{array}$$(C.5)

For later use, we define

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathrm{\Phi }}}^{\prime }(r)\equiv s_{\mathrm{eff}}·|a_N(r)|,{\stackrel{\sim }{\phi }}^{\prime }(r)\equiv s_{\mathrm{eff}}·\sqrt{a_0|a_N(r)|}.\end{array}$$(C.6)

That is, $\stackrel{\sim }{\mathrm{\Phi }}$ carries the Newton-like part of the total acceleration and $\stackrel{\sim }{\phi }$ carries the MOND-like part. The total acceleration can be calculated from either $\stackrel{\sim }{\mathrm{\Phi }}+\stackrel{\sim }{\phi }$ or $\widehat{\mathrm{\Phi }}+\phi$ since their sum is the same,

$$\begin{array}{c}\hfill \widehat{\mathrm{\Phi }}+\phi =\stackrel{\sim }{\mathrm{\Phi }}+\stackrel{\sim }{\phi }.\end{array}$$(C.7)

Appendix D: Approximate analytical solution

The formal solutions Eq. (C.4) are not directly useful since the effective mass M_eff in a_N depends on the value of the fields themselves through $\widehat{\mathrm{\Phi }}+\phi$. Equation (C.4) can, however, be used to recursively solve for M_eff starting at small radii where the baryonic mass M_b dominates. We first derive the recursion formula and then use it to obtain a first order approximation for M_eff. For simplicity, we assume a point-particle baryonic mass distribution, ρ_b(x) = M_bδ(x).

In the following, we split the fields into a part due only to baryons, that is, ${\widehat{\mathrm{\Phi }}}_{\mathrm{b}}$ and φ_b, and the rest, φ_c and ${\widehat{\mathrm{\Phi }}}_c$ which includes the effects of the ghost condensate,

$$\begin{array}{c}\hfill \widehat{\mathrm{\Phi }}\equiv {\widehat{\mathrm{\Phi }}}_{\mathrm{b}}+{\widehat{\mathrm{\Phi }}}_c,\end{array}$$(D.1a)

$$\begin{array}{c}\hfill \phi \equiv {\phi }_{\mathrm{b}}+{\phi }_c,\end{array}$$(D.1b)

where

$$\begin{array}{cc}\hfill {\widehat{\mathrm{\Phi }}}_{\mathrm{b}}& \equiv -f_G\frac{G_N{M}_{\mathrm{b}}}{r},\hfill \end{array}$$(D.2a)

$$\begin{array}{cc}\hfill {\phi }_{\mathrm{b}}& \equiv -(1-f_G)\frac{G_N{M}_{\mathrm{b}}}{r}+\sqrt{G_N{a}_0{M}_{\mathrm{b}}}ln(r/l),\hfill \end{array}$$(D.2b)

for some l. This split depends on the additive constants chosen for ${\widehat{\mathrm{\Phi }}}_{\mathrm{b}}$ and φ_b. In particular, it depends on the choice of l which parametrizes the additive constant in φ_b.

These additive constants can be shuffled around arbitrarily within the physical combination

$$\begin{array}{c}\hfill \mu /Q_0-{\widehat{\mathrm{\Phi }}}_{\mathrm{b}}-{\widehat{\mathrm{\Phi }}}_c-{\phi }_{\mathrm{b}}-{\phi }_c.\end{array}$$(D.3)

To solve the equations of motions, we need to impose a boundary condition for this combination. In practice, we first fixed a value of l, that is, we fixed the additive constant in φ_b. Then, it is equivalent to impose a boundary condition for $\mu /Q_0-{\widehat{\mathrm{\Phi }}}_c-{\phi }_c$. Here, we chose to impose a value for this combination at r = 0,

$$\begin{array}{c}\hfill \mu /Q_0-{\widehat{\mathrm{\Phi }}}_c(0)-{\phi }_c(0).\end{array}$$(D.4)

The effective mass M_eff is the sum of the baryonic mass M_b and the integrated ghost condensate density ρ_c which depends on the combination $\mu /Q_0-\phi -\widehat{\mathrm{\Phi }}$, see Eq. (2). In turn, Eq. (C.5) allows to calculate the derivative of this combination from M_eff. One can use this to derive a recursion formula for M_eff by equating two different expressions for $\mu /Q_0-\phi -\widehat{\mathrm{\Phi }}$.

In order to get $\mu /Q_0-\phi -\widehat{\mathrm{\Phi }}$ from the derivatives from Eq. (C.5), one must integrate once,

$$\begin{array}{cc}\hfill \mu /Q_0-\phi (r)-\widehat{\mathrm{\Phi }}(r)& =\mu /Q_0-{\phi }_{\mathrm{b}}(r)-{\widehat{\mathrm{\Phi }}}_{\mathrm{b}}(r)-{\phi }_c(0)-{\widehat{\mathrm{\Phi }}}_c(0)\hfill \\ \hfill & -{\displaystyle {\int }_0^{r}d{r}^{\prime }({\phi }^{\prime }(r^{\prime })+{\widehat{\mathrm{\Phi }}}^{\prime }(r)-{\phi }_{\mathrm{b}}^{\prime }(r^{\prime })-{\widehat{\mathrm{\Phi }}}_{\mathrm{b}}^{\prime }(r^{\prime })).}\hfill \end{array}$$(D.5)

On the right-hand side, we can plug in Eq. (C.5), that is, $\phi \prime +\widehat{\mathrm{\Phi }}\prime =s_{\text{eff}}(G_N|M_{\mathrm{eff}}|/r^2+\sqrt{a_0{G}_N|M_{\mathrm{eff}}|}/r)$, and the same expression but with M_b instead of M_eff for ${\phi }_{b\prime }+{\widehat{\mathrm{\Phi }}}_{b\prime }$. The left-hand side is proportional to the condensate density ρ_c. Thus, after multiplying by r² and integrating once more, the left-hand side is M_eff. We find,

$$\begin{array}{cc}& \frac{M_{\mathrm{eff}}(x)}{M_{\mathrm{b}}}=1+\alpha {\displaystyle {\int }_0^{x}d{x}^{\prime }{x}^{\prime 2}\{p_l-\frac{1}{a_0{r}_{\mathrm{MOND}}}({\phi }_{\mathrm{b}}(x^{\prime })+{\widehat{\mathrm{\Phi }}}_{\mathrm{b}}(x^{\prime }))}\hfill \\ \hfill & -{\displaystyle {\int }_0^{x^{\prime }}}d{x}^{″}[\frac{1}{x^{″2}}(\frac{M_{\mathrm{eff}}(x^{″})}{M_{\mathrm{b}}}-1)\hfill \\ \hfill & +\frac{1}{x^{″}}(s_{\mathrm{eff}}(x^{″})\sqrt{\frac{|M_{\mathrm{eff}}(x^{″})|}{M_{\mathrm{b}}}}-1)]\},\hfill \end{array}$$(D.6)

where $r_{\mathrm{MOND}}=\sqrt{G_N{M}_{\mathrm{b}}/a_0}$ is the MOND radius. We further defined

$$\begin{array}{c}\hfill x\equiv \frac{r}{r_{\mathrm{MOND}}},\alpha \equiv \frac{m^2}{f_G}{r}_{\mathrm{MOND}}^2,\end{array}$$(D.7)

and the parameter p_l encodes the boundary condition,

$$\begin{array}{c}\hfill p_l\equiv \frac{1}{a_0{r}_{\mathrm{MOND}}}(\frac{\mu }{Q_0}-{\phi }_c(0)-{\widehat{\mathrm{\Phi }}}_c(0)).\end{array}$$(D.8)

The subscript l indicates that the split between φ_b, ${\widehat{\mathrm{\Phi }}}_{\mathrm{b}}$ and φ_c, ${\widehat{\mathrm{\Phi }}}_c$ depends on l.

Equation (D.6) can be used as a recursion formula to iteratively solve for M_eff. The zeroth order approximation is to forget about the ghost condensate density, corresponding to α = 0, and set

$$\begin{array}{c}\hfill M_{\mathrm{eff},0}\equiv M_{\mathrm{b}}.\end{array}$$(D.9)

The first order approximation is obtained from Eq. (D.6) by using the zeroth order approximation on the right-hand side (i.e., by setting M_eff = M_eff, 0 = M_b there),

$$\begin{array}{cc}\hfill M_{\mathrm{eff},1}(x)& \equiv M_{\mathrm{b}}[1+\alpha {\displaystyle {\int }_0^{x}d{x}^{\prime }{x}^{\prime 2}(p_l-\frac{1}{a_0{r}_{\mathrm{MOND}}}({\phi }_{\mathrm{b}}(x^{\prime })+{\widehat{\mathrm{\Phi }}}_{\mathrm{b}}(x^{\prime })))}]\hfill \\ \hfill & =M_{\mathrm{b}}[1+\alpha x^2(\frac{1}{2}+\frac{1}{9}x(3p+1-3ln(x)))],\hfill \end{array}$$(D.10)

where

$$\begin{array}{c}\hfill p\equiv p_l-ln\left(\frac{r_{\mathrm{MOND}}}{l}\right)=p_{l=r_{\mathrm{MOND}}}.\end{array}$$(D.11)

That is, p is the boundary condition for $\mu /Q_0-{\phi }_c-{\widehat{\mathrm{\Phi }}}_c$ for the choice l = r_MOND.

This first order estimate can also be obtained more directly by setting φ = φ_b and $\widehat{\mathrm{\Phi }}={\widehat{\mathrm{\Phi }}}_{\mathrm{b}}$ in the ghost condensate density Eq. (2). The advantage of the recursion formula Eq. (D.6) is that it is, at least conceptually, straightforward to improve on this approximation. And it allows to analytically estimate when the first order approximation breaks down by going to the next higher order, see Appendix J.5.

In the first order approximation, deviations from MOND are proportional to α. This is typically a small number for galaxies. To see this, first note that we typically have $\sqrt{f_G}/m\gtrsim \mathrm{Mpc}$ (Skordis & Złosnik 2021). This is much larger than the MOND radius of galaxies, which typically satisfies r_MOND ≲ 10 kpc. Thus, α is typically smaller than 10⁻⁴ for galaxies,

$$\begin{array}{c}\hfill {\alpha |}_{\mathrm{galaxies}}=\frac{m^2}{f_G}{r}_{\mathrm{MOND}}^2\lesssim {10}^{-4}\ll 1.\end{array}$$(D.12)

In contrast, for galaxy clusters, the MOND radius can be large enough to give α = 𝒪(1).

Appendix E: Numerical solutions

In addition to the analytical approximation discussed above, we also made use of numerical solutions. To obtain these, we used the Julia package ‘OrdinaryDiffEq.jl‘ with the ‘AutoTsit5(Rosenbrock23())‘ method (Rackauckas & Nie 2017; Tsitouras 2011).

We again used the splits $\widehat{\mathrm{\Phi }}={\widehat{\mathrm{\Phi }}}_{\mathrm{b}}+{\widehat{\mathrm{\Phi }}}_c$ and φ = φ_b + φ_c and numerically solved for ${\widehat{\mathrm{\Phi }}}_c$ and φ_c. The equations of motion are

$$\begin{array}{cc}& {\widehat{\mathrm{\Phi }}}_c^{″}+\frac{2{\widehat{\mathrm{\Phi }}}_c^{\prime }}{r}=S\hfill \end{array}$$(E.1)

$$\begin{array}{cc}\hfill & \frac{{\phi }_c^{\prime }}{r}+{\phi }_c^{″}\frac{1+F({\phi }_{\mathrm{b}}+{\phi }_c)}{2}+{\phi }_{\mathrm{b}}^{″}\frac{F({\phi }_{\mathrm{b}}+{\phi }_c)-F({\phi }_{\mathrm{b}})}{2}\hfill \\ \hfill & =\frac{S}{2\stackrel{\sim }{\mu }(|{\phi }_{\mathrm{b}}^{\prime }+{\phi }_c^{\prime }|/a_0)},\hfill \end{array}$$(E.2)

where

$$\begin{array}{cc}\hfill S& \equiv m^2(\mu /Q_0-{\phi }_{\mathrm{b}}-{\phi }_c-{\widehat{\mathrm{\Phi }}}_{\mathrm{b}}-{\widehat{\mathrm{\Phi }}}_c),\hfill \end{array}$$(E.3)

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mu }(s)& =f_G\frac{1+2s(1-f_G)-\sqrt{1+4s(1-f_G)}}{2s{(1-f_G)}^2},\hfill \end{array}$$(E.4)

$$\begin{array}{cc}\hfill F(\phi )& \equiv \frac{{\stackrel{\sim }{\mu }}^{\prime }\left(\frac{|\mathbf{\nabla }\phi |}{a_0}\right)\frac{|\mathbf{\nabla }\phi |}{a_0}}{\stackrel{\sim }{\mu }\left(\frac{|\mathbf{\nabla }\phi |}{a_0}\right)}=\frac{1}{\sqrt{1+4(1-f_G)\frac{|\mathbf{\nabla }\phi |}{a_0}}}.\hfill \end{array}$$(E.5)

On the right-hand side of the φ_c equation, there is in general an additional term proportional to

$$\begin{array}{c}\hfill 4\pi G_N{f}_G{\rho }_{\mathrm{b}}(\frac{1}{2\stackrel{\sim }{\mu }(|{\phi }_{\mathrm{b}}^{\prime }+{\phi }_c^{\prime }|/a_0)}-\frac{1}{2\stackrel{\sim }{\mu }({\phi }_{\mathrm{b}}^{\prime }/a_0)}).\end{array}$$(E.6)

We left out this term since it vanishes for a baryonic point particle, ρ_b = M_bδ(x), which is what we consider here. To see that it vanishes, we first note that it vanishes outside r = 0 because of the factor ρ_b. For r → 0, we have ρ_b ∝ δ(r)/r², ${\phi }_{\text{b}}^{\prime }\propto 1/r$, and ${\phi }_c^{\prime }\propto r$. The behavior of ${\phi }_c^{\prime }$ can, for example, be read off from our approximate analytical solution which is valid at r → 0. Thus, we can expand $\stackrel{\sim }{\mu }(({\phi }_{b\prime }+{\phi }_{c\prime })/a_0)$ around φ_c′ = 0 and Eq. (E.6) becomes proportional to

$$\begin{array}{c}\hfill {\rho }_{\mathrm{b}}\frac{{\stackrel{\sim }{\mu }}^{\prime }({\phi }_{\mathrm{b}}^{\prime }/a_0)}{\stackrel{\sim }{\mu }{({\phi }_{\mathrm{b}}^{\prime }/a_0)}^2}{\phi }_c^{\prime }\to (1/r^2)\delta (r){(1/r)}^{-3/2}·r=\sqrt{r}\delta (r)=0,\end{array}$$(E.7)

where we used that $\stackrel{\sim }{\mu }\prime (s)$ scales as s^−3/2 at large s while $\stackrel{\sim }{\mu }(s)$ becomes constant.

We used a dimensionless length variable y = r/l with l = 1 kpc. We solved the equations in terms of v and u, which are obtained from φ_c and ${\widehat{\mathrm{\Phi }}}_c$, respectively, by rescaling and absorbing the constant μ/Q₀,

$$\begin{array}{c}\hfill \mu /Q_0-{\phi }_c-{\widehat{\mathrm{\Phi }}}_c\equiv -A(u+v),\end{array}$$(E.8)

with A ≡ 10⁻⁷. We used the boundary conditions

$$\begin{array}{cc}\hfill u^{\prime }(0)=v^{\prime }(0)& =0,\hfill \end{array}$$(E.9)

$$\begin{array}{cc}\hfill u(0)+v(0)& =\mathrm{const}.\hfill \end{array}$$(E.10)

The first follows from spherical symmetry, the second corresponds to a choice of chemical potential for the ghost condensate. When comparing the numerical solution to our analytical approximation, the following relation is useful,

$$\begin{array}{c}\hfill u(0)+v(0)=-A^{-1}\sqrt{G_N{M}_{\mathrm{b}}{a}_0}(p+ln(r_{\mathrm{MOND}}/l)).\end{array}$$(E.11)

The logarithm accounts for the fact that p is defined for l = r_MOND while we used l = 1 kpc for the numerical solution. To avoid numerical complications from the 1/r factor in the equations of motion, we did not impose these boundary conditions at r = 0 but at r = 10⁻⁵ kpc, corresponding to y = 10⁻⁵.

Appendix F: Oscillations, m²-f_G-degeneracy

Above, we saw that the effective mass M_eff – and thus also the acceleration – drops to zero at a finite radius. Beyond this radius, the acceleration begins to oscillate (Skordis & Złosnik 2021; Arkani-Hamed et al. 2007), see for example Fig. F.1. As discussed above, this oscillatory regime is potentially unstable since it involved negative energy densities and even negative masses. In this Appendix, we ignore this, keeping in mind that the oscillations might not be physical.

The oscillations depend strongly on the combination m²/f_G which multiplies the condensate density. This is illustrated in Fig. F.1 which shows multiple solutions with the same boundary condition and the same mass M_b but different values of the parameter m²/f_G.

Indeed, in spherical symmetry, the total acceleration in the AeST model depends not on m and f_G separately but only on the combination m²/f_G. This is also illustrated in Fig. F.1 and can be understood from the integro-differential equation Eq. (C.5) that determines the total acceleration. Crucially, Eq. (C.5) depends only on the sum $\widehat{\mathrm{\Phi }}+\phi$ but not on $\widehat{\mathrm{\Phi }}$ and φ separately.

Appendix G: Maximum radius of MOND-like behavior, optimal boundary conditions, relation to maximum of M_eff

The AeST model deviates from MOND at large radii. Here, we discuss quantitatively at which radii these deviations occur and how large they are.

First, we point out that deviations from MOND are inevitable. By this we mean that, for given model parameters and a given baryonic mass M_b, one cannot push the deviations to arbitrarily large radii by judiciously adjusting the boundary condition. This can be seen by considering the effective mass M_eff. If deviations from MOND are small, the zeroth order approximation M_eff, 0 = M_b should be close to our first order approximation M_eff, 1,

$$\begin{array}{c}\hfill \frac{M_{\mathrm{eff},1}(r)}{M_{\mathrm{b}}}=1+\alpha x^2(\frac{1}{2}+\frac{1}{9}x(3p+1-3ln(x))),\end{array}$$(G.1)

where, again, p parametrizes the boundary condition, x = r/r_MOND, and $\alpha =m^2{r}_{\text{MOND}}^2/f_G$. The general form of this function is that it approaches one for r → 0, has a maximum at some finite radius, and then drops to zero (see e.g., Fig. G.1). This definitely deviates from the zeroth order approximation (i.e., from MOND) when it drops to zero. And we can obviously make the radius where this happens arbitrarily large by making the boundary condition p arbitrarily large. So one might wonder whether we can avoid deviations from MOND by just making p arbitrarily large. But this is not possible because the larger p is, the more M_eff deviates from MOND at its maximum. So there must be a finite optimal boundary condition p. For this optimal boundary condition, the radius up to which deviations from MOND stay small is maximized. We refer to this radius as r_max.

More concretely, we considered the ratio of the acceleration $a_{\stackrel{\sim }{\phi }}$ and its value without deviations from MOND, that is, we considered $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$. If one allows this ratio to deviate from MOND by at most a fraction δ,

$$\begin{array}{c}\hfill |\frac{a_{\stackrel{\sim }{\phi }}(r)}{\sqrt{a_0{a}_{\mathrm{b}}(r)}}-1|=|\mathrm{sign}(M_{\mathrm{eff}}(r))\sqrt{\frac{|M_{\mathrm{eff}}(r)|}{M_{\mathrm{b}}}}-1|\stackrel{!}{<}\delta .\end{array}$$(G.2)

That is, if one allows $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}$ to deviate from 1 by at most a fraction δ. Then, there is an optimal boundary condition p that allows this condition to be fulfilled up to a maximum possible radius r_max.

As we show in Appendix J.3, this optimal boundary condition is that for which $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}=a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$ has the value 1 + δ at its maximum. We refer to this radius where $a_{\stackrel{\sim }{\phi }}$ reaches its maximum as r_maxratio. We illustrate the meaning of the quantities r_max, r_maxratio, and δ in Fig. 1.

Instead of $a_{\stackrel{\sim }{\phi }}$, we can also consider the total acceleration $a_{\mathrm{tot}}=a_{\stackrel{\sim }{\phi }}+a_{\stackrel{\sim }{\mathrm{\Phi }}}$. This gives an analogous result: If we allow the total acceleration a_tot to deviate from MOND by at most a fraction δ_tot, then the optimal boundary condition is that for which $a_{\mathrm{tot}}/(a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}})$ has the value 1 + δ_tot at its maximum.

Strictly speaking, the optimal boundary conditions for $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$ and $a_{\mathrm{tot}}/(a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}})$ differ. At least for galaxies, however, the optimal boundary conditions are almost identical for both cases. Moreover, the maximizing radius and the value at this maximum are almost identical. This is shown in Appendix J.2. Thus, for our purposes, it does not matter much whether we consider the optimal boundary condition for $a_{\stackrel{\sim }{\phi }}$ or for a_tot.

In Fig. G.2 we show the relation between the maximum allowed deviation from MOND δ and the maximum radius r_max up to which this condition can be fulfilled. For the numerical solutions, we found r_max by maximizing the radius where a solution first deviates by more than a fraction δ from MOND using the Julia package ‘Optim.jl‘ (Mogensen & Riseth 2018).

As discussed above, for the optimal boundary conditions, the ratio $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$ has the value 1 + δ at its first maximum. The radius r_maxratio where this maximum occurs is closely related to the radius r_max. This is illustrated in Fig G.3 which shows r_maxratio as a function of δ. This is very similar to Fig. G.2 just with the values on the y-axis a bit larger. In particular, r_maxratio corresponds to $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}=1+\delta$, while r_max corresponds to $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}=1-\delta$.

In Appendix J.2 and Appendix J.4, we give an analytical estimate for the relation between δ, r_max, and r_maxratio. For galaxies, a good approximation is

$$\begin{array}{c}\hfill \frac{r_{\mathrm{maxratio}}}{r_{\mathrm{MOND}}}\approx {\left(9\frac{{(1+\delta )}^2-1}{r_{\mathrm{MOND}}^2{m}^2/f_G}\right)}^{1/3},r_{\max }\approx 1.53r_{\mathrm{maxratio}}.\end{array}$$(G.3)

Thus, where this approximation is valid, the only difference between Fig. G.3 and Fig. G.2 is a factor 1.53 in the y-axis values.

Our estimate for r_max can be compared to a related estimate from Skordis & Złosnik (2021). There, the authors estimate that the AeST model acceleration is MOND-like up to a critical radius r_c,

$$\begin{array}{c}\hfill \frac{r_c}{r_{\mathrm{MOND}}}\sim {\left(\frac{1}{m^2{r}_{\mathrm{MOND}}^2}\right)}^{1/3}.\end{array}$$(G.4)

This has the same scaling in r_MOND and m as our estimate. However, our estimate improves on this in a few ways. First, the m² factor in r_c should be m²/f_G. Otherwise, the estimate does not correctly take into account the difference between G_N and $\widehat{G}$. Second, our version comes with worked out prefactors, including a parameter that controls how big deviations actually are.

We emphasize again that the maximum radius r_max – or the critical radius r_c – corresponds to choosing boundary conditions that are optimal for reproducing MOND (given an allowed deviation δ). Galaxies are not guaranteed to reproduce MOND up to that radius. In general, galaxies will end up with boundary conditions different from the optimal one and deviate from MOND already at smaller radii.

In Fig. G.2 and Fig. G.3, we show our analytical approximation only where we analytically estimate that it deviates by at most q = 10% from the full solution. This validity estimate works by going to the next-higher order in our method for approximating M_eff and then comparing to our first-order approximation, see Appendix J.5. Essentially, our approximation is valid for small deviations from MOND δ, but not for larger ones, as can already be guessed from Fig. G.1.

From Fig. G.2 and Fig. G.3 we see that more massive galaxies can remain close to MOND for longer than less massive galaxies. However, in the context of MOND, accelerations are often more important than radii. Thus, in Fig. 2, we show the minimum acceleration $a_{b,\text{min}}=G_N{M}_{\text{b}}/r_{\text{max}}^2$ corresponding to the maximum radius r_max. We see that the trend is now inverted due to the additional factor of M_b in a_b. Less massive galaxies can have MOND-like behavior down to smaller accelerations than more massive galaxies.

Appendix H: Stacking for weak lensing

In the weak-lensing analysis of Brouwer et al. (2021), the central quantity is the stacked excess surface density (ESD) profile ΔΣ_stacked. Stacking here means taking a weighted average over the galaxy sample. More specifically,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{stacked}}=\frac{{\mathrm{\Sigma }}_i{W}_i\left(\frac{1}{1+\mu }{ϵ}_{t,i}{\mathrm{\Sigma }}_{\mathrm{crit},i}\right)}{{\mathrm{\Sigma }}_i{W}_i},\end{array}$$(H.1)

where the W_i are weights, Σ_crit is the critical surface density, ϵ_t is the ellipticity, and the factor 1 + μ calibrates the shear estimates. The ellipticity ϵ_t of a galaxy is the sum of its intrinsic ellipticity ${ϵ}_t^{\text{int}}$ and the tangential shear γ_t caused by weak lensing. The intrinsic ellipticities average to zero in a large sample, so that ΔΣ_stacked measures the tangential shear γ_t. For an individual lens and up to the calibration factor 1 + μ, the combination γ_tΣ_crit is given by the ESD ΔΣ,

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)=\frac{2\pi {\int }_0^{R}d{R}^{\prime }\mathrm{\Sigma }(R^{\prime })}{\pi R^2}-\mathrm{\Sigma }(R),\end{array}$$(H.2)

where Σ(R) is the surface density corresponding to the lensing mass M_lens. Thus, for a galaxy sample with known surface densities, we can calculate ΔΣ_stacked by simply averaging the individual ESD profiles of each galaxy,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{stacked}}=\frac{{\mathrm{\Sigma }}_i{W}_i\mathrm{\Delta }{\mathrm{\Sigma }}_i}{{\mathrm{\Sigma }}_i{W}_i}.\end{array}$$(H.3)

The stacked ESD profile is linear in the surface densities Σ_i. These surface densities are calculated linearly from the density ρ_lens that produces the lensing mass M_lens. This lensing mass is defined by Φ′(r) = G_NM_lens/r². In our case,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{lens}}(r)=\frac{1}{4\pi G_N{r}^2}{\partial }_r\left(r^2{\partial }_r(\widehat{\mathrm{\Phi }}+\phi )\right).\end{array}$$(H.4)

Thus, the stacked ESD profile is linear in $\widehat{\mathrm{\Phi }}+\phi$. As a consequence, the total acceleration a_tot inferred from the stacked ESD profile is just the weighted average of the total accelerations a_tot, i of each stacked galaxy,

$$\begin{array}{c}\hfill a_{\mathrm{tot},\mathrm{stacked}}(R)=\frac{{\mathrm{\Sigma }}_i{W}_i{a}_{\mathrm{tot},i}(R)}{{\mathrm{\Sigma }}_i{W}_i}.\end{array}$$(H.5)

In order to calculate a stacked RAR, we should not stack at a fixed position R but at a fixed acceleration a_b = G_NM_b/R². So we should instead use

$$\begin{array}{c}\hfill a_{\mathrm{tot},\mathrm{stacked}}(a_{\mathrm{b}})=\frac{{\mathrm{\Sigma }}_i{W}_i{a}_{\mathrm{tot},i}\left(\sqrt{\frac{G_N{M}_{b,i}}{a_{\mathrm{b}}}}\right)}{{\mathrm{\Sigma }}_i{W}_i}.\end{array}$$(H.6)

For a fixed baryonic mass M_b, stacking in position space and acceleration space is equivalent.

We now illustrate how a stacked weak-lensing RAR in the AeST model might look like. For simplicity, we assume that all weights W_i are the same and we further assume that all galaxies have the same baryonic mass M_b. Essentially, we consider stacking a sample of galaxies that differ only in their boundary conditions. Then, we have

$$\begin{array}{c}\hfill a_{\mathrm{tot},\mathrm{stacked}}(R)=N^{-1}{\mathrm{\Sigma }}_i{a}_{\mathrm{tot},i}(R),\end{array}$$(H.7)

where N is the number of galaxies in the sample.

One can also consider a large number of galaxies with boundary conditions p distributed uniformly in an interval [p₁, p₂]. Then, we can write

$$\begin{array}{c}\hfill a_{\mathrm{tot},\mathrm{stacked}}(R)=\frac{1}{p_2-p_1}{\displaystyle {\int }_{p_1}^{p_2}dp{a}_{\mathrm{tot},p}(R).}\end{array}$$(H.8)

As long as all accelerations a_tot, p(R) point in the usual direction, that is, as long as M_eff is positive, this integral can be done analytically for our first order analytical approximation,

$$\begin{array}{cc}& {a_{\mathrm{tot},\mathrm{stacked}}(R)|}_{M_{\mathrm{eff}}>0}=a_{\mathrm{b}}\frac{M_{\mathrm{eff}}(x,\overline{p})}{M_{\mathrm{b}}}\hfill \\ \hfill & +\frac{\sqrt{a_0{a}_{\mathrm{b}}}}{p_2-p_1}\frac{2}{\alpha x^3}({\left(\frac{M_{\mathrm{eff}}(x,p_2)}{M_{\mathrm{b}}}\right)}^{3/2}-{\left(\frac{M_{\mathrm{eff}}(x,p_1)}{M_{\mathrm{b}}}\right)}^{3/2}),\hfill \end{array}$$(H.9)

where x = r/r_MOND, $\overline{p}=(p_1+p_2)/2$, and M_eff(x, p) is our first order analytical approximation Eq. (D.10).

Appendix I: Higher-order terms in condensate density

The ghost condensate density from Eq. (2) is linear in the fields φ and $\widehat{\mathrm{\Phi }}$. In contrast, for their cosmological calculations, Skordis & Złosnik (2021) used a ghost condensate density including nonlinear corrections. Here, we discuss whether using the linearized form around galaxies is valid and explain our decision to do so.

In the full action, not assuming the quasi-static weak-field limit, the condensate density corresponds to a term

$$\begin{array}{c}\hfill \mathcal{K}(Q)={\mathcal{K}}_2{(Q-Q_0)}^2,\end{array}$$(I.1)

where Q = A^μ∇_μϕ is a scalar combination of the normalized vector field A^μ and the scalar field ϕ. In the quasi-static weak-field limit, ϕ = Q₀ ⋅ t + φ. The constant 𝒦₂ is related to the mass parameter m by $m=\sqrt{2{\mathcal{K}}_2/(2-K_B)}{Q}_0$.

Our expression for the ghost condensate Eq. (2) is linear because the function 𝒦(Q) is quadratic. Indeed, in the quasi-static weak-field limit,

$$\begin{array}{c}\hfill Q-Q_0\approx Q_0·(\frac{\dot{\phi }}{Q_0}-\phi -\widehat{\mathrm{\Phi }}).\end{array}$$(I.2)

However, a quadratic function 𝒦(Q) cannot simultaneously satisfy cosmological constraints and support a MOND-like regime in galaxies (Skordis & Złosnik 2021). In particular, cosmological observations require that the ghost condensate’s equation of state w satisfies w ≲ 0.02 at a scale factor a = 10⁻⁴ (Ilić et al. 2021). But if we assume m²/f_G ≲ 1/Mpc in order to have a MOND-like phenomenology around galaxies, w is too large. Indeed, using 0 < K_B < 2 and f_G < 1, (Skordis & Złosnik 2021)

$$\begin{array}{c}\hfill w\approx \frac{3H_0^2{\mathrm{\Omega }}_0}{2(2-K_B)f_G(m^2/f_G)a^3}\gtrsim {10}^{-8}{a}^{-3},\end{array}$$(I.3)

where H₀ is the Hubble constant and Ω₀ is the matter density parameter today. This is illustrated in Fig. I.1 which shows the cosmological ghost condensate density as a function of the scale factor a for m = 1 Mpc⁻¹ (see the Q₀/𝒵₀ = 0 line, the parameter 𝒵₀ is discussed below). We see that a dust-like evolution is possible at late times but not around a ∼ 10⁻⁴.

To avoid this problem, Skordis & Złosnik (2021) introduced two alternative forms of 𝒦(Q),

$$\begin{array}{cc}\hfill {\mathcal{K}}_{\exp }(Q)& ={\mathcal{K}}_2{\mathcal{Z}}_0^2(e^{{\mathcal{Z}}^2}-1),\hfill \end{array}$$(I.4)

$$\begin{array}{cc}\hfill {\mathcal{K}}_{\cosh }(Q)& =2{\mathcal{K}}_2{\mathcal{Z}}_0^2(cosh(\mathcal{Z})-1),\hfill \end{array}$$(I.5)

where 𝒵₀ is a constant and

$$\begin{array}{c}\hfill \mathcal{Z}\equiv \frac{Q-Q_0}{{\mathcal{Z}}_0}.\end{array}$$(I.6)

These suppress the equation of state at early times where 𝒵 is large and reduce to the quadratic 𝒦(Q) at small 𝒵.

In our galaxy-scale calculations above, we assumed the quadratic form of 𝒦(Q). This is justified only if 𝒵 is sufficiently small. To see whether or not this is the case, we first note that a typical value of $\dot{\phi }/Q_0-\phi -\widehat{\mathrm{\Phi }}$ around galaxies is 10⁻⁷. Thus, a typical value of 𝒵 around galaxies is $\mathcal{Z}\sim \frac{Q_0}{{\mathcal{Z}}_0}·{10}^{-7}$. That is, we can use the quadratic 𝒦(Q) as long as

$$\begin{array}{c}\hfill \frac{Q_0}{{\mathcal{Z}}_0}\lesssim {10}^7.\end{array}$$(I.7)

Interestingly, this condition is not satisfied, or only barely satisfied, for the explicit cosmological perturbation calculations in Skordis & Złosnik (2021). In particular, Skordis & Złosnik (2021) used Q₀/𝒵₀ = 10⁸ for their example calculation using 𝒦_cosh and Q₀/𝒵₀ = 10¹³ for 𝒦_exp.

We nevertheless assume the quadratic form around galaxies for the following reasons. First, this seems to be what the authors of the model had in mind. After all, they use the quadratic form of 𝒦(Q) when deriving the action for the quasi-static weak-field limit (Skordis & Złosnik 2021).

Second, large values of Q₀/𝒵₀ do not seem to be required in order to satisfy cosmological constraints. Indeed, the constraint w ≲ 0.02 at a = 10⁻⁴ can easily be satisfied with a much smaller value of Q₀/𝒵₀. This is illustrated for 𝒦_exp in Fig. I.2 which shows that Q₀/𝒵₀ = 𝒪(1) is sufficient.

Even if w satisfies w ≲ 0.02 at a = 10⁻⁴, it can be larger at later times and potentially violate observational constraints that apply at these later times. Indeed, w(a) is not a monotonous function, as is illustrated in Fig. I.3. Its maximum is set by Q₀/𝒵₀ and m controls at which scale factor a this maximum occurs. Still, we see from Fig. I.3 that Q₀/𝒵₀ = 𝒪(10) is sufficient for w to satisfy w ≲ 0.02 even at later times, thus easily satisfying the constraints at these times from Ilić et al. (2021). That is, these constraints can easily be satisfied with Q₀/𝒵₀ ≪ 10⁷ which allows using the quadratic 𝒦(Q) around galaxies. We find a similar result for 𝒦_cosh.⁴

A third reason is that, if the nonlinearities of 𝒦_exp or 𝒦_cosh were to be important around galaxies, then either many galaxies would deviate very strongly from MOND or there would be almost no deviations from MOND at all. Both cases do not require a thorough investigation here. Strong deviations from MOND are ruled out by observations while pure-MOND predictions are already discussed elsewhere.

It remains to explain why 𝒦_exp and 𝒦_cosh have the effects described in the previous paragraph. We first consider 𝒦_exp and 𝒦_cosh with the parameter m having roughly the value we assumed above for the quadratic function 𝒦(Q), that is, m ∼ 1 Mpc⁻¹. Then, as soon as nonlinearities become important, the condensate density is enhanced exponentially compared to the case of a quadratic 𝒦(Q). Thus, deviations from MOND set in exponentially earlier, that is, the radius r_max is exponentially smaller. For example, the oscillations discussed in Appendix F start at exponentially smaller radii. This is in conflict with, for example, the observed RAR which is MOND-like as discussed above.

One way out would be to make the prefactor m² of 𝒦_exp exponentially smaller to keep the ghost condensate from becoming large. Indeed, by choosing m² sufficiently small, we can get rid of any significant deviations from MOND around galaxies. The predictions of the AeST model in this case are just those of MOND so do not require a special investigation.

A middle ground would be to make m² exponentially smaller, but by just the right amount so that deviations from MOND set in on roughly galactic scales. However, the maximum radius r_max would still be exponentially sensitive to the galactic potential $\mu /Q_0-\phi -\widehat{\mathrm{\Phi }}$. Thus, for galaxies with slightly larger $\mu /Q_0-\phi -\widehat{\mathrm{\Phi }}$, deviations from MOND still set in at an exponentially smaller radius compared to the quadratic 𝒦(Q). Galaxies with slightly smaller $\mu /Q_0-\phi -\widehat{\mathrm{\Phi }}$ would perfectly follow the MOND predictions up to exponentially large radii. Thus, a significant fraction of galaxies should still deviate strongly from MOND which is not what is observed.

Appendix J: A few useful properties of M_eff, 1(r)

In this appendix, we consider a few properties of the analytical first order approximation M_eff, 1 from Eq. (D.10),

$$\begin{array}{c}\hfill M_{\mathrm{eff},1}(r)=M_{\mathrm{b}}[1+\alpha x^2(\frac{1}{2}+\frac{1}{9}x(3p+1-3ln(x)))].\end{array}$$(J.1)

We sometimes write M_eff instead of M_eff, 1 for brevity.

J.1. Increasing functions of boundary condition

We first show that both $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$ and $a_{\mathrm{tot}}/(a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}})$ as well as their spatial derivatives are increasing functions of the boundary condition p. At least as long as M_eff is positive. This technical result will be useful in Appendix J.3.

We first consider $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}=\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}$ as a function of p. Both M_eff(r) and $M_{\text{eff}}^{\prime }(r)$ are increasing functions of the boundary condition, that is, of p, for all r. This follows from α and r being positive. We then consider $\sqrt{M_{\mathrm{eff}}(r)}$ and ${\partial }_r\sqrt{M_{\mathrm{eff}}(r)}$ for radii where M_eff is positive. In this case, $\sqrt{M_{\mathrm{eff}}(r)}$ is an increasing function of p everywhere because M_eff is,

$$\begin{array}{c}\hfill {\partial }_p\sqrt{M_{\mathrm{eff}}(r)}=\frac{{\partial }_p{M}_{\mathrm{eff}}}{2\sqrt{M_{\mathrm{eff}}(r)}}>0.\end{array}$$(J.2)

The same holds for the spatial derivative ${\partial }_r\sqrt{M_{\mathrm{eff}}(r)}$ but it is a bit harder to see. We find

$$\begin{array}{cc}& {\partial }_p{\partial }_x\sqrt{M_{\mathrm{eff}}(r)/M_{\mathrm{b}}}\hfill \\ \hfill & =\frac{3\alpha x^2}{\sqrt{2}}\frac{9+\alpha x^2(\frac{9}{2}+x(3p+1-3ln(x)))+9+\alpha x^2((6-\frac{9}{2})+x)}{{(18M_{\mathrm{eff}}(r)/M_{\mathrm{b}})}^{3/2}}\hfill \\ \hfill & >\frac{3\alpha x^2}{\sqrt{2}}\frac{9M_{\mathrm{eff}}(r)/M_{\mathrm{b}}}{{(18M_{\mathrm{eff}}(r)/M_{\mathrm{b}})}^{3/2}}>0,\hfill \end{array}$$(J.3)

which follows because r and M_eff are positive.

We consider next, again assuming M_eff > 0, the function

$$\begin{array}{c}\hfill \frac{a_{\stackrel{\sim }{\phi }}+a_{\stackrel{\sim }{\mathrm{\Phi }}}}{a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}}}=\frac{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}\frac{1}{x}+\sqrt{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}}}{\frac{1}{x}+1}\equiv \frac{H/x+\sqrt{H}}{1/x+1},\end{array}$$(J.4)

where we defined H = M_eff/M_b. Since M_eff is an increasing function of the boundary condition p, the same holds for this function. The slope of this function with respect to x is also an increasing function of p, but this is again a bit harder to see. We have

$$\begin{array}{cc}& {\partial }_p{\partial }_x\frac{H/x+\sqrt{H}}{1/x+1}\hfill \\ \hfill & =x^2\alpha \frac{2xH(4+3x)+4(3+2x)H^{3/2}-x^2(1+x){\partial }_{x}H}{12{(1+x)}^2{H}^{3/2}},\hfill \end{array}$$(J.5)

where we already inserted the explicit expressions for ∂_pH and ∂_p∂_xH due to their simple form. It now suffices to show positivity of the following expression,

$$\begin{array}{cc}& 2H(4+3x)-x(1+x){\partial }_{x}H\hfill \\ \hfill & =2H(4+3x)-x^2(1+x)\alpha (1+px-xln(x)).\hfill \end{array}$$(J.6)

In general, this can be negative. However, here it cannot be due to our assumption that M_eff is positive. To see this, we first use the definition of H to rewrite x(p − ln(x)) in terms of H,

$$\begin{array}{c}\hfill x(p-ln(x))=-\frac{3}{2}-\frac{x}{3}+\frac{3}{\alpha x^2}(H-1).\end{array}$$(J.7)

With this, we find

$$\begin{array}{cc}& 2H(4+3x)-x(1+x){\partial }_{x}H\hfill \\ \hfill & =H(5+3x)+\frac{1}{6}(1+x)(18+3x^2\alpha +2x^3\alpha )>0,\hfill \end{array}$$(J.8)

which completes the proof.

J.2. Maximum

We first consider the maximum of $a_{\stackrel{\sim }{\phi }}(r)/\sqrt{a_0{a}_{\mathrm{b}}(r)}=\sqrt{M_{\mathrm{eff}}(r)/M_{\mathrm{b}}}$. The maximum of both M_eff(r) and $\sqrt{M_{\mathrm{eff}}(r)}$ is determined by $M_{\text{eff}}^{\prime }(r)=0$. This condition can be written as

$$\begin{array}{c}\hfill p=ln(x)-\frac{1}{x},\end{array}$$(J.9)

which can be solved in terms of the Lambert W function,

$$\begin{array}{c}\hfill x=\frac{1}{W(e^{-p})}.\end{array}$$(J.10)

Using this condition gives further for the maximum of $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}$ with value 1 + δ,

$$\begin{array}{c}\hfill \frac{x^2}{6}+\frac{x^3}{9}=\frac{{(1+\delta )}^2-1}{\alpha }.\end{array}$$(J.11)

This can be solved exactly but the resulting expression is not particularly illuminating. The size of the right-hand side determines whether the x² or the x³ term dominates. For large values, the x³ term dominates, for small values the x² term dominates. For large right-hand sides, we have for the maximum

$$\begin{array}{c}\hfill x={\left(9\frac{{(1+\delta )}^2-1}{\alpha }\right)}^{1/3}-\frac{1}{2}+\mathcal{O}\left({\left(\frac{\alpha }{{(1+\delta )}^2-1}\right)}^{1/3}\right).\end{array}$$(J.12)

The boundary condition needed to obtain such a maximum is

$$\begin{array}{c}\hfill p=\frac{1}{3}ln\left(9\frac{{(1+\delta )}^2-1}{\alpha }\right)+\mathcal{O}\left({\left(\frac{\alpha }{{(1+\delta )}^2-1}\right)}^{1/3}\right).\end{array}$$(J.13)

We consider next the maximum of $a_{\mathrm{tot}}/(a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}})$, that is, the maximum of the function

$$\begin{array}{c}\hfill \frac{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}\frac{1}{x}+\sqrt{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}}}{\frac{1}{x}+1}\equiv \frac{H/x+\sqrt{H}}{1/x+1},\end{array}$$(J.14)

assuming M_eff > 0. The maximum of this function is both at a similar location and has a similar value as that of $\sqrt{M_{\mathrm{eff}}}$ – at least for galaxies. We first show that the maximum is at a similar location. To this end, we set the derivative with respect to x to zero, plug in x + δ_x where x is the maximizer of $\sqrt{M_{\mathrm{eff}}}$, expand to lowest order in δ_x, use H′(x) = 0, and solve for δ_x. This gives

$$\begin{array}{c}\hfill \frac{{\delta }_x}{x}={(\frac{2x}{1+x}-\frac{9{(1+x)}^2(2+\delta )(2+x+2\delta )}{x(3+2x){(1+\delta )}^2})}^{-1},\end{array}$$(J.15)

where we have additionally used Eq. (J.11), H(x) = (1 + δ)², H″(x) = − (1 + x)α, and we have eliminated α by using Eq. (J.11), giving α = ((1 + δ)² − 1)/(x²/6 + x³/9).

The simplest way to understand this expression is to plot it as a function of x for a number of values of δ. For δ ≪ 1, we have that |δ_x/x| is always smaller than 2.5%. For δ = 𝒪(1), it can become a bit larger. For example, up to δ = 1.5, we have that |δ_x/x| is always smaller than 5%. But even for extremely large values of δ, we have that |δ_x/x| is always smaller than 14.5%. Thus, the maximum of $a_{\mathrm{tot}}/(a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}})$ is indeed at a similar location as that of $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$.

Now we can evaluate the maximum value of the function Eq. (J.14), assuming that the maximizer is close to the one of $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}$,

$$\begin{array}{c}\hfill \frac{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}\frac{1}{x}+\sqrt{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}}}{\frac{1}{x}+1}\approx \frac{{(1+\delta )}^2+x(1+\delta )}{1+x}.\end{array}$$(J.16)

We can write this as

$$\begin{array}{c}\hfill 1+{\delta }_{\mathrm{tot}}\approx \frac{{(1+{\delta }_{\stackrel{\sim }{\phi }})}^2+x(1+{\delta }_{\stackrel{\sim }{\phi }})}{1+x}=(1+{\delta }_{\stackrel{\sim }{\phi }})(1+\frac{{\delta }_{\stackrel{\sim }{\phi }}}{1+x}),\end{array}$$(J.17)

where 1 + δ_tot is the maximum value of Eq. (J.14) and $1+{\delta }_{\stackrel{\sim }{\phi }}$ is the maximum of $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$. We consider first small ${\delta }_{\stackrel{\sim }{\phi }}$. Then, the factors relative to MOND are approximately the same in both cases, that is, $1+{\delta }_{\mathrm{tot}}\approx 1+{\delta }_{\stackrel{\sim }{\phi }}$. (But we note: δ_tot and δ_φ themselves might differ by up to a factor 2 for small x). When ${\delta }_{\stackrel{\sim }{\phi }}$ is not small, x is very large, see Eq. (J.11). At least for galaxies, which have α ≪ 1. Thus, the factors relative to MOND are again approximately the same in both cases – unless ${\delta }_{\stackrel{\sim }{\phi }}$ is very large. But for such large δ, our approximation is anyway no longer valid, see Appendix J.5.

J.3. Optimal boundary condition

We consider the ratio

$$\begin{array}{c}\hfill \frac{a_{\stackrel{\sim }{\phi }}}{\sqrt{a_0{a}_{\mathrm{b}}}}=\sqrt{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}},\end{array}$$(J.18)

and require that it stays between 1 + δ and 1 − δ for some positive δ. One may wonder what the boundary condition is that allows this requirement to be fulfilled up to the maximum possible radius r_max. We refer to this as the optimal boundary condition. Here, we show that it is the boundary condition for which $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}$ has the value 1 + δ at its maximum, which we discussed in Appendix J.2.

To this end, we consider this solution for which $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}$ has the maximum value 1 + δ at some radius r₁. This solution will have the value 1 − δ at some radius r₂ > r₁ since M_eff drops to zero after the maximum. We refer to this solution’s boundary condition p as p₁.

Here, we consider M_eff > 0. When M_eff < 0, the gravitational force in the AeST model points into the opposite direction from that in MOND. We are not interested in this regime here. Thus, we assume δ < 1 which implies M_eff > 0.

Any solution with a boundary condition p larger than p₁ has a larger value $\sqrt{M_{\mathrm{eff}}}$ everywhere (see Appendix J.1). This implies that such solutions must have $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}>1+\delta$ at r = r₁. Thus, boundary conditions p larger than p₁ cannot be optimal.

We then consider solutions with boundary condition p smaller than p₁. Such solutions have a smaller value $\sqrt{M_{\mathrm{eff}}}$ and slope ${\partial }_r\sqrt{M_{\mathrm{eff}}}$ everywhere (see Appendix J.1). Thus, these solutions reach their maximum earlier than r₁ and this maximum has a value smaller than 1 + δ. Starting at this lower maximum, such solutions then go to zero faster (since their slope is more negative). Thus, these solutions reach $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}=1-\delta$ earlier than r₂. Thus, smaller boundary conditions also cannot be optimal.

It follows that the optimal boundary condition is that which gives $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}=1+\delta$ at the maximum of M_eff. That is, we have p = p₁.

The preceding argument implicitly considered negative condensate densities by considering radii beyond the radius where M_eff is maximal. One may wonder whether anything changes when restricting solutions to have positive densities. But this is not the case, at least as long as we just cut off the solutions when the density reaches zero and do not continue them with something else (such as a Navarro-Frenk-White halo). The above argument can be straightforwardly adapted to this case.

Instead of $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$ we could also consider the total acceleration relative to MOND. That is,

$$\begin{array}{c}\hfill \frac{a_{\stackrel{\sim }{\phi }}+a_{\stackrel{\sim }{\mathrm{\Phi }}}}{a_{\mathrm{b}}+\sqrt{a_0{a}_{\mathrm{b}}}}=\frac{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}\frac{1}{x}+\sqrt{\frac{M_{\mathrm{eff}}}{M_{\mathrm{b}}}}}{\frac{1}{x}+1}.\end{array}$$(J.19)

The argument above can be adapted to this case as well. That is, the optimal boundary boundary condition is that for which this total acceleration ratio has the maximum value 1 + δ.

The optimal boundary conditions for $a_{\stackrel{\sim }{\phi }}$ are, in general, not the same as those for the total acceleration $a_{\stackrel{\sim }{\phi }}+a_{\stackrel{\sim }{\mathrm{\Phi }}}$. However, in practice, they are very similar and give very similar deviations from MOND at very similar maximum radius. This follows from the properties of the maximum discussed in Appendix J.2.

J.4. Relation between r_max and r_maxratio

Here, we consider a solution that is optimal for some δ. As discussed above and as illustrated in Fig. 1, the radius r_maxratio is where the ratio $a_{\stackrel{\sim }{\phi }}/\sqrt{a_0{a}_{\mathrm{b}}}$ has its maximum value, namely 1 + δ, and r_max is the radius where this ratio drops to 1 − δ. Here, we consider the relation between these two radii. We argue that, usually, r_max is larger than r_maxratio by a fixed factor 1.53 for small δ.

Based on numerical examples, we first guess a factor around 1.5. Then, we calculate first-order corrections to this factor 1.5 and find that corrections are, in a certain regime, small. This confirms the validity of our initial guess. In addition, we find that the first-order corrections are universal and give a factor 1.53.

To find r_max, we must solve $\sqrt{M_{\mathrm{eff}}/M_{\mathrm{b}}}=1-\delta$. That is,

$$\begin{array}{c}\hfill 1+\alpha x^2(\frac{1}{2}+\frac{1}{9}x(3p+1-3ln(x)))={(1-\delta )}^2.\end{array}$$(J.20)

where x = r_max/r_MOND. We write

$$\begin{array}{c}\hfill x=f{x}_0(1+ϵ),\end{array}$$(J.21)

where x₀ is the radius where this equation is satisfied with (1 + δ) instead of (1 − δ) on the right-hand side, that is to say it is the radius where M_eff is maximal, that is, x₀ = r_maxratio/r_MOND. We later set f = 3/2 and let ϵ parametrize deviations from a factor 1.5 between r_max and r_maxratio, that is, later we find f(1 + ϵ)≈1.53.

We plug this expression for x back into the equation determining x, expand in ϵ, and solve for ϵ. Then, we eliminate first p and then α using the two conditions we know from x₀ being a maximum, namely Eq. (J.9) and Eq. (J.11). This gives

$$\begin{array}{c}\hfill ϵ=-\frac{1}{3}+\frac{(2-\delta )(3+2x_0)+3f^2(2+\delta )+2f^3{x}_0(2+\delta )}{18f^2(2+\delta )(f-1+f{x}_{0}ln(f))}.\end{array}$$(J.22)

We now consider f = 3/2 and check whether or not ϵ is actually small. For large x₀, we find

$$\begin{array}{c}\hfill {ϵ}_{f=3/2}=\frac{70+19\delta -81(2+\delta )ln(3/2)}{243(2+\delta )ln(3/2)}+\mathcal{O}(1/x_0).\end{array}$$(J.23)

The maximum radius x₀ is large as long as δ/α ≫ 1, see Appendix J.2. For galaxies, we typically have α ≲ 10⁻⁴ (see Appendix G). Thus, this condition is fulfilled unless δ is extremely large. For such large δ, the first order analytical approximation for M_eff is anyway no longer valid, as we show in Appendix J.5. Thus, for galaxies, the large-x₀ result is typically valid.

If, in addition, we have δ ≪ 1, we find

$$\begin{array}{c}\hfill {ϵ}_{f=3/2}=-\frac{1}{3}+\frac{35}{243ln(3/2)}+\mathcal{O}(\delta )+\mathcal{O}(1/x_0).\end{array}$$(J.24)

Thus,

$$\begin{array}{c}\hfill \frac{r_{\max }}{r_{\mathrm{maxratio}}}=\frac{3}{2}(1+{ϵ}_{f=3/2})\approx 1+\frac{35}{162ln(3/2)}\approx 1.53.\end{array}$$(J.25)

J.5. Validity of approximation

We now consider the validity of our first order approximation M_eff, 1 for M_eff from Eq. (D.10). Figure G.1 and other plots suggest that it is often a reasonable approximation but also that it becomes worse for larger deviations from MOND δ. Here, we estimate analytically when this approximation is valid.

To this end, we calculate the 2nd order approximation to M_eff and check when the first and second orders deviate from each other. To facilitate a quantitative analytical estimate, we do not calculate the 2nd order approximation exactly. We make additional approximations and assumptions. The result should still be a reasonable estimate of when the first order approximation is valid.

The 2nd order approximation is obtained by plugging the first order approximation into the right-hand side of Eq. (D.6),

$$\begin{array}{cc}\hfill \frac{M_{\mathrm{eff},2}(x)}{M_{\mathrm{b}}}& =\frac{M_{\mathrm{eff},1}(x)}{M_{\mathrm{b}}}-\alpha {\displaystyle {\int }_0^{x}d{x}^{\prime }{x}^{\prime 2}{\int }_0^{x^{\prime }}d{x}^{″}[\frac{1}{x^{″2}}(\frac{M_{\mathrm{eff},1}(x^{″})}{M_{\mathrm{b}}}-1)}\hfill \\ \hfill & +\frac{1}{x^{″}}(s_{\mathrm{eff},1}(x^{″})\sqrt{\frac{|M_{\mathrm{eff},1}(x^{″})|}{M_{\mathrm{b}}}}-1)],\hfill \end{array}$$(J.26)

where the subscripts 1 and 2 refer to the first and second order approximations, respectively.

The double integrals in this expression can be done analytically for the part linear in M_eff, 1. We find

$$\begin{array}{cc}\hfill \frac{M_{\mathrm{eff},2}^{\mathrm{lin}}(x)}{M_{\mathrm{b}}}& =\frac{M_{\mathrm{eff},1}(x)}{M_{\mathrm{b}}}\hfill \\ \hfill & -\frac{x^4{\alpha }^2}{1800}(225+(62+60p)x-60xln(x)).\hfill \end{array}$$(J.27)

Below we express radii relative to where M_eff, 1 has its maximum,

$$\begin{array}{c}\hfill \widehat{x}\equiv \frac{r}{r_{\mathrm{maxratio}}}\equiv \frac{x}{x_{\mathit{mr}}}.\end{array}$$(J.28)

We can also eliminate p by using the relation p = ln(x_mr)−1/x_mr. This gives

$$\begin{array}{cc}\hfill \frac{M_{\mathrm{eff}}^{\mathrm{lin}}(\widehat{x})}{M_{\mathrm{b}}}& \approx \frac{M_{\mathrm{eff},1}(\widehat{x})}{M_{\mathrm{b}}}\hfill \\ \hfill & -\frac{{\widehat{x}}^4{x}_{\mathrm{mr}}^4{\alpha }^2}{1800}(225-60\widehat{x}+62x_{\mathrm{mr}}\widehat{x}(1-\frac{60}{62}ln(\widehat{x}))).\hfill \end{array}$$(J.29)

We consider $\widehat{x}$ between 0 and 1. The function $\widehat{x}(1-60/62ln(\widehat{x}))$ grows from 0 to 1 in this $\widehat{x}$ interval. Thus, when x_mr is large, that is, when r_maxratio is large compared to r_MOND, this term will dominate, except at very small radii $\widehat{x}$. Otherwise, the $225-60\widehat{x}$ term dominates.

We now consider the remaining terms in M_eff, 2/M_b, that is, those with $\sqrt{M_{\mathrm{eff},1}/M_{\mathrm{b}}}$. Here, the double integrals are not straightforward to do analytically. Therefore, we make an additional approximation and we restrict ourselves to radii r before M_eff, 1 reaches its maximum, that is, we only consider r ≤ r_maxratio. Concretely, we assume that $\sqrt{M_{\mathrm{eff},1}/M_{\mathrm{b}}}=a_{\stackrel{\sim }{\phi },1}/\sqrt{a_0{a}_{\mathrm{b}}}$ grows linearly from 1 at r = 0 to its maximum 1 + δ at r = r_maxratio,

$$\begin{array}{c}\hfill \frac{a_{\stackrel{\sim }{\phi },1}(r)}{\sqrt{a_0{a}_{\mathrm{b}}(r)}}=\sqrt{\frac{M_{\mathrm{eff},1}(r)}{M_{\mathrm{b}}}}=1+\delta \frac{r}{r_{\mathrm{maxratio}}},r\le r_{\mathrm{maxratio}}.\end{array}$$(J.30)

From Fig. 1 we see that this underestimates M_eff, 1 at small radii and then overestimates it at larger radii toward r_maxratio. This is a rough approximation but we have verified numerically that it captures the validity of the first order approximation reasonably well. With this, we have

$$\begin{array}{c}\hfill \frac{M_{\mathrm{eff},2}^{\mathrm{sqrt}}(x)}{M_{\mathrm{b}}}=-\frac{\alpha \delta }{4}{x}_{\mathrm{mr}}^3{\widehat{x}}^4.\end{array}$$(J.31)

The remaining task is to sum $M_{\mathrm{eff},2}^{\mathrm{lin}}$ and $M_{\mathrm{eff},2}^{\mathrm{sqrt}}$ and determine when this sum deviates from M_eff, 1. More specifically, we are interested mostly in $a_{\stackrel{\sim }{\phi }}$ so we are interested in the square root. We parametrize the deviation by a parameter q,

$$\begin{array}{c}\hfill \sqrt{\frac{M_{\mathrm{eff},2}(\widehat{x})}{M_{\mathrm{eff},1}(\widehat{x})}}\stackrel{!}{=}1-q.\end{array}$$(J.32)

This relation determines the radius $\widehat{x}$ where the 2nd order becomes important if we want to trust the first order approximation up to a fraction q.

Unfortunately, this equation is not straightforward to solve analytically. Thus, we make further approximations. First, we use the linear approximation $\sqrt{M_{\mathrm{eff},1}/M_{\mathrm{b}}}=1+\delta \widehat{x}$ also in the denominator of Eq. (J.32) and the first term in $M_{\mathrm{eff},2}^{\mathrm{lin}}$. Second, we neglect the $\widehat{x}(1-(60/62)ln(\widehat{x}))$ term in $M_{\mathrm{eff},2}^{\mathrm{lin}}$.

The reasoning behind the second of these approximations is the following. As argued above, the $\widehat{x}(1-(60/62)ln(\widehat{x}))$ term is important only for large x_mr. But in this case, the double-integral contributions to $M_{\mathrm{eff}}^{\mathrm{lin}}$ are anyway small compared to $M_{\mathrm{eff}}^{\mathrm{sqrt}}$ so it does not matter. Indeed,

$$\begin{array}{c}\hfill {\frac{M_{\mathrm{eff},2}^{\mathrm{lin}}-M_{\mathrm{eff},1}}{M_{\mathrm{eff},2}^{\mathrm{sqrt}}}|}_{x_{\mathrm{mr}}\gg 1}\approx \frac{31}{25}\frac{1}{x_{\mathrm{mr}}}\frac{{(1+\delta )}^2-1}{\delta }\widehat{x}(1-\frac{60}{62}ln(\widehat{x})),\end{array}$$(J.33)

where we used $\alpha =9({(1+\delta )}^2-1)/x_{\text{mr}}^3$ which is valid for x_mr ≫ 1, see Eq. (J.11). Since $\widehat{x}\le 1$ and x_mr ≫ 1, this ratio is small unless δ ≳ x_mr. Thus, for our purposes this ratio is small. We are not interested in very large deviations from MOND δ which would be required for δ ≳ x_mr ≫ 1.

Equation (J.32) then becomes

$$\begin{array}{c}\hfill \frac{-\frac{\alpha \delta }{4}{x}_{\mathrm{mr}}^3{\widehat{x}}^4-\frac{{\alpha }^2}{8}{\widehat{x}}^4{x}_{\mathrm{mr}}^4(1-\frac{4}{15}\widehat{x})}{{(1+\delta \widehat{x})}^2}\stackrel{!}{=}{(1-q)}^2-1.\end{array}$$(J.34)

This is still not solvable analytically. Thus, our final approximation will be to neglect the $-4/15\widehat{x}$ term. This is again a rough approximation but, as mentioned above, we have verified numerically that the result is reasonable. Thus, we need to solve

$$\begin{array}{c}\hfill \sqrt{\frac{1}{4}\alpha x_{\mathrm{mr}}^3\delta +\frac{1}{8}{\alpha }^2{x}_{\mathrm{mr}}^4}\frac{{\widehat{x}}^2}{1+\delta \widehat{x}}\stackrel{!}{=}\sqrt{1-{(1-q)}^2}.\end{array}$$(J.35)

We finally find for the radius r_approx up to which the first order approximation is valid to within a fraction q,

$$\begin{array}{c}\hfill \frac{r_{\mathrm{approx}}}{r_{\mathrm{maxratio}}}=\frac{1}{2}(\beta \delta +\sqrt{{(\beta \delta )}^2+4\beta }),\end{array}$$(J.36)

where

$$\begin{array}{c}\hfill \beta =\frac{\sqrt{2}}{3}\sqrt{\frac{1-{(1-q)}^2}{{\delta }^2(2+\delta )}}\sqrt{\frac{{(3+2x_{\mathrm{mr}})}^2}{(3+2x_{\mathrm{mr}})x_{\mathrm{mr}}+9(2+\delta )}},\end{array}$$(J.37)

where we used Eq. (J.11) to eliminate α.

We note that we assumed r ≤ r_maxratio when deriving this estimate. Thus, if our expression for r_approx/r_maxratio is larger than 1, the first order approximation is valid up to at least r_maxratio but the precise value of r_approx/r_maxratio is not meaningful.

The most important implication of our estimate is that the first order approximation is good for small δ (except for extremely restrictive values of q). For larger δ, the approximation is worse. Explicitly, for small δ the quantity β scales as

$$\begin{array}{c}\hfill {\beta }_{\delta \ll 1}\sim \sqrt{q}/\delta \end{array}$$(J.38)

so that r_approx/r_maxratio scales as

$$\begin{array}{c}\hfill {\frac{r_{\mathrm{approx}}}{r_{\mathrm{maxratio}}}|}_{\delta \ll 1}\sim \sqrt{\beta }\sim {\left(\frac{q}{{\delta }^2}\right)}^{1/4}.\end{array}$$(J.39)

When δ is small, this ratio is large, so the first order approximation is valid. If we demand a better approximation (i.e., a smaller q), the radius r_approx becomes smaller.

In the opposite limit, δ ≫ 1,

$$\begin{array}{c}\hfill {\beta }_{\delta \gg 1}\sim \sqrt{q/{\delta }^3}\mathrm{or}\sqrt{q/{\delta }^4},\end{array}$$(J.40)

where the power 3 applies for sufficiently large x_mr and the power 4 for sufficiently small x_mr. In both cases does r_approx/r_maxratio tend to zero for large δ. So the first order approximation is not valid much.

All Tables

All Figures