<figure><img src="image_20.png" /><figcaption>**Figure 19.** The median $\tilde{v}$ as a function of the median $r_{\text{sky}}/r_{\text{M}}$ for the nominal sample of Chae (2023), which requires both proper motion components of both stars in each WB to be more precise than 1 per cent (left-hand panel). This is tightened to 0.3 per cent in the right-hand panel, as done in their appendix B and Chae (2024). The sample size and number of bins are shown at the top. Similarly to Fig. 11, each point corresponds to the same number of WBs apart from the last point, which has a slightly different number due to rounding.</figcaption></figure>

A critical consideration for the WBT is the accuracy of $v_{\text{sky}}$, which is the main source of uncertainty in $\tilde{v}$. We quantify its uncertainty very carefully in order to limit the $\tilde{v}$ uncertainty of the WBs in our sample (Section 2.4.6). As argued there, since the main peak to the $\tilde{v}$ distribution lies at $\tilde{v} \approx 0.5$ (Fig. 12) and $\tilde{v}$ uncertainties would add in quadrature to the intrinsic width of the $\tilde{v}$ distribution, we need the uncertainty in $\tilde{v}$ to be $\lesssim 0.1$ in the regime most relevant to the WBT. An uncertainty of 0.1 would inflate the dispersion from $0.5$ to $\sqrt{0.5^2 + 0.1^2} = 0.51$, implying only a 2 per cent broadening – too little to affect the inferred gravity law. While this approach might seem somewhat conservative, we have seen that it still leaves a reasonable sample size, partly because we impose a less strict limit to the $\tilde{v}$ uncertainty outside the main peak region [equation (10)].

Unfortunately, Chae (2023) do not estimate the uncertainty in $\tilde{v}$ or in $v_{\text{sky}}$. Their approach is to require a maximum 1 per cent uncertainty on the heliocentric proper motions $\mu_{\alpha,*}^A, \mu_{\alpha,*}^B, \mu_{\delta}^A$, and $\mu_{\delta}^B$ of the stars in each WB. To see why this approach is not sufficient, suppose that both stars have a heliocentric velocity of $30 \text{ km s}^{-1}$ in the East-West and North-South directions (the small difference due to WB orbital motion is not relevant here). The sky-projected heliocentric velocity of each star could then have an uncertainty of $300 \text{ m s}^{-1}$ along each direction. Since we need to combine four proper motion components to calculate $v_{\text{sky}}$ [equation (32)], this has an uncertainty of $\approx 300\sqrt{4}$ or $600 \text{ m s}^{-1}$. However, the Newtonian $v_c$ of two Sun-like stars with a $10 \text{ kAU}$ separation is only $420 \text{ m s}^{-1}$, implying a very poorly constrained $\tilde{v}$. This completely undermines the WBT because its main idea is that some typical measure of $\tilde{v} \propto \sqrt{g/g_N}$, making it absolutely essential to have precise constraints on $\tilde{v}$.

To try and address this issue, Chae (2023) conduct an extension analysis in their appendix B where the 1 per cent requirement mentioned above is tightened to 0.3 per cent. This is unlikely to be enough because in the previous argument, we would have a maximum allowed $\sigma(v_{\text{sky}})$ of $180 \text{ m s}^{-1}$ rather than $600 \text{ m s}^{-1}$, where $\sigma(Q)$ is the uncertainty on any quantity $Q$. However, if the Newtonian $v_c$ is $420 \text{ m s}^{-1}$ and the mode of the $\tilde{v}$ distribution is at 0.5, the intrinsic velocity dispersion is $210 \text{ m s}^{-1}$, which is not much larger than the upper limit to $\sigma(v_{\text{sky}})$. Alternatively, the plausible scenario described above would lead to $\sigma(\tilde{v}) = 180/420 = 0.43$, which is still too large for the WBT. Even so, we repeat the median $\tilde{v}$ analysis with all four relevant heliocentric proper motions required to have an accuracy of at least 0.3 per cent. The sample size is almost halved

by this much stricter cut, so we use only half as many bins to reduce the scatter. Despite this, we still achieve a good range in $r_{\text{sky}}/r_{\text{M}}$ and probe well into the asymptotic regime (right-hand panel of Fig. 19). A rising trend is evident here, though the agreement with the MOND expectation is somewhat less good. Importantly, it is not really possible to fit the results with a flat line, as would be expected in Newtonian gravity.

A more rigorous but still very simple way to estimate the uncertainty in $v_{\text{sky}}$ is provided by equation 5 of El-Badry, Rix & Heintz (2021). If we neglect other sources of error like perspective effects and the small error in $\bar{d}$, we can approximate that the only uncertain term is $\mu_{\text{rel}}$. Its uncertainty is

$$ \sigma(\mu_{\text{rel}}) = \frac{\sqrt{\sigma^2(\mu_{\alpha,*}^{\text{rel}}) + \sigma^2(\mu_{\delta}^{\text{rel}})}}{\mu_{\text{rel}}}, \quad (33) $$

$$ \sigma(\mu_{\alpha,*}^{\text{rel}}) = \sqrt{\sigma^2(\mu_{\alpha,*}^A) + \sigma^2(\mu_{\alpha,*}^B)}, \quad (34) $$

$$ \sigma(\mu_{\delta}^{\text{rel}}) = \sqrt{\sigma^2(\mu_{\delta}^A) + \sigma^2(\mu_{\delta}^B)}. \quad (35) $$

This result is based on the fact that at first order in $\sigma(\mu_{\text{rel}})/\mu_{\text{rel}}$, the uncertainty in $\mu_{\text{rel}}$ arises only from the proper motion uncertainty along the direction $(\mu_{\alpha,*}^{\text{rel}}, \mu_{\delta}^{\text{rel}})$ because uncertainty in $\mu_{\text{rel}}$ in the orthogonal direction mainly affects its direction rather than its magnitude. This approximation becomes inaccurate if $\sigma(\mu_{\text{rel}}) \gtrsim \mu_{\text{rel}}$, so a Monte Carlo approach might be preferable even if we could assume that the uncertainty arises entirely from uncorrelated astrometric errors in the proper motions.

With the above approximations, the uncertainty in $v_{\text{sky}}$ is simply

$$ \sigma(v_{\text{sky}}) = \bar{d}\sigma(\mu_{\text{rel}}). \quad (36) $$

Dividing this by the Newtonian $v_c$ then yields an estimated $\tilde{v}$ uncertainty. We suggest that such an approach would capture the major uncertainties affecting $\tilde{v}$ because the mass and systemic distance should both have small uncertainties, while perspective effects should likewise have little effect (Section 2.3.1).

We can now check whether the result of Chae (2023) is robust against a quality cut based on the estimated $\tilde{v}$ uncertainty of each WB. We begin by quantifying the proportion of WBs in each $r_{\text{sky}}/r_{\text{M}}$ bin with $\sigma(\tilde{v}) > 0.1 \text{ max}(1, \tilde{v}/2)$, the cut imposed in our own analysis [equation (10)]. The result is shown in Fig. 20. It is clear that most WBs in the asymptotic regime fail to pass our quality cut on $\sigma(\tilde{v})$.
<figure><img src="image_21.png" /><figcaption>**Figure 20.** The fraction of WBs in each $r_{sky}/r_M$ bin of the nominal Chae (2023) sample with $\sigma(\tilde{v}) > 0.1 \max(1, \tilde{v}/2)$. This applies to most of the WBs in the asymptotic regime critical to testing MOND.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Similar to the left-hand panel of Fig. 19, but with the Chae (2023) nominal sample restricted to WBs with an estimated $\sigma(\tilde{v}) < 0.1 \max(1, \tilde{v}/2)$ according to equations (33) and 36. The median $\tilde{v}$ is now flat with respect to our proxy for the WB acceleration thanks to this cut on the estimated $\tilde{v}$ uncertainty, which is critical to reliably conducting the WBT.</figcaption></figure>

This raises the possibility that the MOND-like trend in Fig. 19 is caused by measurement errors broadening the $\tilde{v}$ distribution at low accelerations.

To analyse the data set of Chae (2023) more robustly, we restrict their nominal sample to the same limiting $\sigma(\tilde{v})$ as we used in our own analysis [equation (10)]. The result is shown in Fig. 21. The sample size only decreases by about 1/5, allowing us to continue using 50 bins in $r_{sky}/r_M$, albeit with a somewhat reduced range given that the sample is reduced by half at the MOND radius and slightly more beyond that. Even so, we still have many WBs well into the asymptotic regime. These show that the median $\tilde{v}$ is now almost completely flat with respect to $r_{sky}/r_M$. There is clearly a major disagreement with the dashed grey line showing the MOND expectation calibrated to the high-acceleration end. This dramatic difference is caused by removing only 1/5 of the WBs, indicating that the issue we identified is indeed extremely serious.

The MOND-like rising median $\tilde{v}$ trend in the nominal sample of Chae (2023) persists despite an apparently tight 0.3 per cent cut on the proper motion errors, which leads to an even smaller sample size than in Fig. 21. This highlights that the important quantity for the

WBT is the relative velocity of the stars in each WB rather than the heliocentric velocity of each star. While one can envisage that a sufficiently restrictive cut on the latter would be sufficient, this is simply not targeting the relevant quantity, so the sample size is likely to become insufficient before the cut has been adequately tightened for the purpose of the WBT. As an exercise in this direction, we found that requiring a 0.1 per cent precision reduces the sample to only 4977 WBs and makes the median $\tilde{v}$ nearly flat with respect to $r_{sky}/r_M$, so perhaps this cut is tight enough. However, the sample size can be much larger while still maintaining sufficient accuracy on $\tilde{v}$ provided that the quality cut is targeted at $\tilde{v}$, or at least at $\mu_{rel}$.

The importance of this issue stems from the fact that for the same $\sigma(v_{sky})$ in m s⁻¹, $\sigma(\tilde{v}) \propto \sqrt{r_{sky}}$ because the Newtonian $v_c \propto 1/\sqrt{r_{sky}}$, making it increasingly difficult to pass a strict quality cut on $\sigma(\tilde{v})$ (Fig. 20). Without such a cut, as we get to large $r_{sky}$ and thus also $r_{sky}/r_M$, the $\tilde{v}$ distribution might be broadened by measurement errors to a greater extent, pushing up the median $\tilde{v}$ and creating a fake Milgromian signal. This raises the more general issue that any systematic trend with $r_{sky}$ can appear to be an acceleration-dependent trend. In principle, one could disentangle the two because our proxy for the WB internal acceleration depends on both $r_{sky}$ and $M$. However, the former is by far the main factor because of its large range and the fact that $r_M$ scales only as $\sqrt{M}$, making the WB mass relatively less important. There is also an issue with the data set of Chae (2023) pertaining to the masses that makes it difficult to use a broad range (Appendix E). One therefore has to be cautious before concluding that any trend discovered in the WB data is best understood as causally related to their internal acceleration.

Given our results in this section, it seems extremely unlikely that local WBs prefer Milgromian over Newtonian gravity once the WB sample is selected carefully to ensure that $\tilde{v}$ is accurately known. It should be straightforward to alter the quality cut related to the astrometric accuracy by focusing on how precisely we know the relative velocity within each WB, since this is what enters into the calculation of the physically relevant $\tilde{v}$ parameter.<fnref n="16" /> We stress that even if $\tilde{v}$ is not calculated explicitly, the WBT is ultimately about the ratio between the relative velocity and the Newtonian circular orbit prediction because the WB orbits are too long to directly measure the acceleration (Banik & Kroupa 2019a). The uncertainty on this velocity ratio is therefore of critical importance to the WBT and must be quantified, with systems rejected if the uncertainty is too large.

This concern likely also applies to the analysis of Hernandez (2023): its section 4 states that $v_{sky}$ should exceed 1.5× its uncertainty. This translates to a fractional $\tilde{v}$ accuracy of 2/3, so a system with $\tilde{v} = 1$ might have $\sigma(\tilde{v}) = 0.6$ (this concern also applies to Hernandez et al. 2024). Addressing this issue similarly to our study would reduce the already small sample size. However, our work suggests that the loss might be modest, especially given the high overall quality of the WB sample used in their study.

## 5.3 The MOND interpolating function

Our results for the WBT depend on the assumed MOND interpolating function because local WBs are subject to the Galactic external field of magnitude $1.8 a_0$ (Section 3.1; see also section 7.1 of BZ18). Thus, our falsification of MOND could in principle be avoided with a sufficiently rapid transition to Newtonian behaviour when $g_N > a_0$.

<footnote n="16">¹⁶This would allow to relax the requirement for the heliocentric proper motion components to have a 1 per cent precision, increasing the sample size somewhat.</footnote>
<table><caption>**Table 4.** The enhancement factor $\eta$ to the radial gravity [equation (15)] and our corresponding estimate for $\alpha_{grav}$ [equation (39)] with different MOND formulations and interpolating functions. Results differ little between AQUAL and QUMOND and between the MLS and simple functions.</caption><thead><tr><th rowspan="2">Interpolating function</th><th colspan="2">AQUAL</th><th colspan="2">QUMOND</th></tr><tr><th>η</th><th>α<sub>grav</sub></th><th>η</th><th>α<sub>grav</sub></th></tr></thead><tbody><tr><td>Simple</td><td>1.4056</td><td>0.96</td><td>1.4228</td><td>1</td></tr><tr><td>MLS</td><td>1.3508</td><td>0.84</td><td>1.3692</td><td>0.88</td></tr><tr><td>Standard</td><td>1.0661</td><td>0.17</td><td>1.0726</td><td>0.18</td></tr><tr><td>Sharp</td><td>1</td><td>0</td><td>1</td><td>0</td></tr></tbody></table>

However, this is in strong tension with rotation curve constraints. This is because the enhancement to $g_N$ needed to fit the host galaxy rotation curve is very similar to the predicted enhancement to the gravity binding a WB in the asymptotic regime. This remains true in both AQUAL and QUMOND, which give almost the same boost for the same interpolating function (Fig. 7). We clarify this by briefly giving the main analytic results relevant to EFE-dominated WBs and thereby estimate the predicted $\alpha_{grav}$ with different MOND formulations and interpolating functions (for further details, see sections 2.2 and 7.1 of BZ18).

The WBT is mostly sensitive to $\eta$ [equation (15)] in the asymptotic regime. In QUMOND, equation 16 of BZ18 tells us that

$$ \eta = v_e \left(1 + \frac{K_e}{3}\right), \quad K_e \equiv \frac{\partial \ln v_e}{\partial \ln g_{N,e}}, \qquad (37) $$

where $v_e$ is the QUMOND interpolating function if we set its argument to the Newtonian-equivalent Galactic gravity $g_{N,e}$ [equation (14)]. This result matches equation 57 of Milgrom (2010).

AQUAL does not directly use $v$ but instead uses the interpolating function $\mu$, which is defined so that $\mu g = g_N$ in spherical symmetry. We can talk of AQUAL and QUMOND theories as having the 'same' interpolating function if they give the same relation between $g$ and $g_N$ in spherical symmetry, which requires that $\mu(g)v(g_N) = 1$. According to equation 35 of BZ18, the AQUAL version of equation (37) is

$$ \eta = \frac{\tan^{-1}\sqrt{L_e}}{\mu_e\sqrt{L_e}}, \quad L_e \equiv \frac{\partial \ln \mu_e}{\partial \ln g_e}, \qquad (38) $$

where $L_e$ is the AQUAL interpolating function considering only the Galactic gravity $g_e$. We can find $L_e$ from the QUMOND $K_e$ using the relation $(1 + L_e)(1 + K_e) = 1$ (equation 38 of BZ18). The angle-averaged enhancement to $g_N$ shown above matches that stated in equation 65 of Milgrom (2010) once we use the relation $\sin^2\theta = \tan^2\theta/(\tan^2\theta + 1)$.

The observable in the WBT is the relative velocity and this is related to the circular velocity, which scales as $\sqrt{\eta}$. Moreover, the parameter $\alpha_{grav}$ only captures deviations from the Newtonian result because it is defined to be 0 in the Newtonian case, when $\eta \equiv 1$. Thus, we can approximate that the predicted $\alpha_{grav} \propto \sqrt{\eta} - 1$. Since the normalization must be chosen such that $\alpha_{grav} = 1$ if we use QUMOND with the simple interpolating function,

$$ \alpha_{grav} = \frac{\sqrt{\eta} - 1}{0.193}. \qquad (39) $$

Table 4 shows $\eta$ and the estimated $\alpha_{grav}$ for different MOND formulations and interpolating functions. Differences between AQUAL and QUMOND are very small, as stated in section 5.1 of Milgrom (2010) and explained further in section 7.2 of BZ18, which is devoted to a comparison of the two MOND formulations. The simple and MLS ($v^{-1} = 1 - \exp(-\sqrt{g_N}/a_0)$; McGaugh 2008; Famaey &

McGaugh 2012; McGaugh, Lelli & Schombert 2016) interpolating functions give similar predictions for the WBT because they are numerically very similar. The main difference is that the MLS function converges to Newtonian gravity much faster when $g \gtrsim 10a_0$ due to the exponential term, but this is not important for the WBT because WB observations are sensitive to the gravity binding it rather than merely how much that deviates from $g_N$.<fnref n="17" /> Compared to the simple and MLS cases, the standard interpolating function significantly reduces the predicted signal for the WBT, though even then the predicted $\alpha_{grav} = 0.18$ is still much higher than inferred by our nominal MCMC analysis (bottom right panel of Fig. 14). But since some variations to the modelling assumptions do increase the inferred $\alpha_{grav}$ somewhat (Table 3), it can be argued that the WBT is marginally consistent with the standard interpolating function

$$ v^2 = \frac{1}{2} + \sqrt{\frac{1}{4} + \left(\frac{a_0}{g_N}\right)^2}. \qquad (40) $$

However, this is in tension with rotation curve constraints (section 7.1 of BZ18). We revisit their arguments in more detail below.

Equations (37) and (38) show that $\eta \approx v_e = 1/\mu_e$ up to a modest correction for azimuthal averaging. This is because $v_e$ can be arbitrarily large, but $K_e$ is always between $-1/2$ and 0 while $L_e$ is between 1 and 2, with low $K_e$ and high $L_e$ corresponding to the deep-MOND limit ($K_e = -0.26$ in the Solar neighbourhood). As a result, the enhancement to the orbital velocities of WBs in the asymptotic regime should be only slightly less than the enhancement to the Newtonian baryonic rotation curve of the parent galaxy in the vicinity of the WBs. This means a very stringent null detection of MOND effects in the WBT can only be reconciled with MOND if the Galactic rotation curve at the Solar circle has almost the same amplitude as it would do in Newtonian gravity with baryons alone.

We test this in Fig. 22, where the cyan squares with error bars show binned observational results and the dotted blue vertical line shows the Solar circle radius of 8.2 kpc (the uncertainty on this is very small; Gravity Collaboration 2019). The largest uncertainty in the kinematics comes from the peculiar (non-circular) velocity of the Sun and hence on the actual value of $v_e, \circ$, which normalizes the entire curve. However, the resulting uncertainty has been substantially mitigated by the direct measurement of the Solar System's acceleration relative to distant quasars (Gaia Collaboration 2021b), which is a direct probe of $g_e$. Besides the kinematics, there is also some systematic uncertainty in the Newtonian rotation curve of the Galactic disc because we do not perfectly know its actual mass-to-light ratio and gas content. With these small caveats in mind, the Newtonian rotation curve of the baryons alone (orange line) is clearly below the actual rotation curve at the Solar circle, which is also much flatter. Since scaling the Newtonian rotation curve (e.g. to raise the mass-to-light ratio) would make the predicted decline even steeper, it is clear that MOND can only work with an interpolating function that significantly enhances Newtonian gravity at the Solar circle. Indeed, the solid green line shows that QUMOND with the simple interpolating function fits rather well – even the bumps and wiggles apparent in the observations are reproduced nicely if we consider similar features in the Galactic surface density profile (McGaugh 2018). However, the required significant Milgromian enhancement to the local amplitude of the rotation curve contradicts

<footnote n="17">¹⁷The simple interpolating function provides a good match to the velocity dispersion profiles of elliptical galaxies up to almost 100 $a_0$ (Chae & Gong 2015; Chae et al. 2019, 2020a), beyond which the MOND enhancement to $g_N$ must be very small.</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure 22.** The MW rotation curve as measured by averaging different observational determinations (cyan squares with error bars; Eilers et al. 2019; Mróz et al. 2019; Chrobáková et al. 2020) and as predicted in different gravity theories using the baryonic mass distribution from Wang, Hammer & Yang (2022). Notice that the Newtonian model with only baryons (orange line) falls well short of the observed rotation curve at the Solar circle radius of 8.2 kpc (dotted blue vertical line). This deficiency can be rectified by adding dark matter (blue line), but the halo parameters must be inferred from the data and are not predicted *a priori*. An alternative solution with very little flexibility is QUMOND with the simple interpolating function [green line; equation (13)], but this only works because this function predicts a significant enhancement to Newtonian gravity at the relevant acceleration (compare the simple and sharp interpolating functions in Table 4). This enhancement should also be detectable in local WBS. Reproduced from fig. 1 of Zhu et al. (2023), modified by Haixia Ma.</figcaption></figure>

our stringent null detection of MOND effects in the WBT, which implies an enhancement by at most a few percent. This contradiction is underpinned by the directly measured acceleration of the Solar System relative to distant quasars (Gaia Collaboration 2021b), which precludes substantial deviations from the kinematically deduced $g_e$ (Section 3.1). If we very conservatively assume that $v_{c, \odot}$ is at most 10 km s$^{-1}$ higher than our adopted 232.8 km s$^{-1}$, then the stronger EFE would only reduce the QUMOND $\alpha_{grav}$ for the simple interpolating function from 1 to 0.92, while an even higher $v_{c, \odot} = 250$ km s$^{-1}$ would give $\alpha_{grav} = 0.87$. Such a large deviation of $v_{c, \odot}$ from our adopted value is almost inconceivable given the many decades of research on this issue. Turning instead to the possibility of a much sharper interpolating function, the required enhancement to the baryonic surface density would need to be truly substantial, which seems rather unlikely.

Another interesting aspect of Fig. 22 is the failure of Moffat gravity (MOG; Moffat 2006), whose prediction is shown as the red line. MOG modifies gravity only beyond a certain distance, so it passes the WBT (see section 2.3 of Roshan et al. 2021a). However, MOG underpredicts the Galactic rotation curve (as shown previously by Negrelli et al. 2018). Moreover, a joint fit to the velocity dispersion profile and star formation history of Dragonfly 44 rules out MOG at 5.5$\sigma$ confidence (Haghi et al. 2019). This shows how difficult it is for a modified gravity theory to remain consistent with the WBT and simultaneously explain galaxy dynamics, even if we neglect extragalactic data.

Of course, one should also consider the rather precise constraints on MOND available nowadays from extragalactic rotation curves. In particular, the Spitzer Photometry and Accurate Rotation Curves (SPARC; Lelli, McGaugh & Schombert 2016) catalogue reveals a tight RAR, which we show in Fig. 23 as the ratio between the observed $g$ at fixed $g_N$ and the prediction $g_v$ with different interpo-

<figure><img src="image_24.png" /><figcaption>**Figure 23.** The ratio between the observed $g$ at fixed baryonic $g_N$ (McGaugh, Lelli & Schombert 2016) and the prediction $g_v$ according to different MOND interpolating functions for the canonical $a_0 = 1.2 \times 10^{-10}$ m s$^{-2}$ (Begeman, Broeils & Sanders 1991; Gentile, Famaey & de Blok 2011; Lelli et al. 2017b). Observational uncertainties are found by dividing the dispersion in $\log_{10} g$ at each $g_N$ by $\sqrt{N}$, where $N$ is the number of data points in the bin. Due to the numerical similarity of the simple and MLS functions (black and blue lines, respectively); error bars are omitted for the latter, but they are shown on the standard interpolating function [red; equation (40)] and for the case of an infinitely sharp transition between the Newtonian and MOND regimes [green; equation (41)]. The dashed vertical line shows the Solar circle $g_N$ assuming the simple interpolating function adopted for this study [equation (13)].</figcaption></figure>

lating functions.<fnref n="18" /> The error bars show the logarithmic dispersion in $g_N$ about a smooth relation scaled down by $\sqrt{N}$, where $N$ is the number of data points in the $g_N$ bin. It is clear that the MLS and simple functions are numerically rather similar and both match the observations fairly well. However, the standard interpolating function [equation (40)] deviates very substantially from observations (this is also apparent in fig. 1 of Tian & Ko 2019). As with other analyses, the crucial assumption is the stellar mass-to-light ratio from stellar population synthesis models (Schombert, McGaugh & Lelli 2022), which would need substantial modifications to allow the standard interpolating function to match the data.

This discrepancy is worsened by using an infinitely sharp transition between the Newtonian and deep-MOND regimes, i.e.

$$v = \max \left( 1, \sqrt{\frac{a_0}{g_N}} \right). \qquad (41)$$

This yields no enhancement to gravity at the Solar circle ($\alpha_{grav} = 0$), so the results of the WBT work best with this function given that our posterior inference on $\alpha_{grav} = -0.021^{+0.065}_{-0.045}$. Even the modest enhancement predicted by the standard interpolating function causes 3$\sigma$ tension with our WB results. At the same time, this enhancement is too little to explain extragalactic rotation curves across almost the full range of $g_N$ that they probe (Fig. 23). Similar conclusions can be drawn from the Galactic rotation curve, which barring major systematics shows that $g > g_N$ by about 50 per cent at the Solar circle even though both slightly exceed $a_0$ (Fig. 22). Therefore, it seems impossible for MOND to simultaneously match galaxy dynamics

<footnote n="18">¹⁸A unique relation is not expected in modified gravity theories because disc galaxies are not spherically symmetric, but deviations from the spherically symmetric relation should be very small (Jones-Smith et al. 2018; Chae 2022).</footnote>
and Gaia data on local WBs. This falsifies classical modified gravity versions of MOND at high significance.

## 5.4 Broader implications

So far, we have assumed that MOND should be understood as a modification to Newtonian gravity in the weak-field regime [equation (1)]. Given the failure of this approach in local WBs, it is helpful to consider MOND as a modification to inertia at low accelerations while preserving Newtonian gravity (Milgrom 1994, 2011, 2022). This may avoid significant MOND effects in local WBs because their internal orbital motion has a much higher frequency than their Galactocentric orbit (e.g. see fig. 14 of BZ18). We can consider this in terms of an effective Newtonian external field $g_{N,e}^{eff}$, defined such that the boost to the gravity binding a WB is still given by equation (37), but with argument $g_{N,e}^{eff}$ instead of the Newtonian gravity $g_{N,e}$ sourced by the rest of the Galaxy. Since the MOND model considered in the main part of our study predicts a 42 per cent enhancement to the average radial gravity (Table 4) and our results show that $\alpha_{grav} < 0.2$ at quite high confidence (Table 3), we can interpret the WBT as constraining the average gravity binding a WB to within 8 per cent of the Newtonian prediction in the asymptotic regime where the Galactic gravity dominates over the internal gravity. The corresponding lower limit on $g_{N,e}^{eff}$ is 7.8 $a_0$ (4.6 $a_0$) for the simple (MLS) interpolating function. In both cases, this is several times larger than the expected value of $\approx 1.2 a_0$ (Section 3.1.1). The high $g_{N,e}^{eff}$ could indicate that the Galactic EFE suppresses the gravity of a local WB by much more than we have assumed, making it very nearly Newtonian.

Modified inertia theories of MOND are still at an early stage of development, so it is not yet possible to conduct detailed simulations of situations like interacting galaxies with a low degree of symmetry. Even so, focusing on systems where the predictions are clear led to the first attempt to distinguish whether MOND is better understood as a modification to gravity or inertia (Petersen & Lelli 2020). The results mildly prefer modified gravity at $\approx 1.5\sigma$ confidence. Modified inertia theories make the very strong prediction that test particles on circular orbits will have a unique relation between their kinematic acceleration and the gravity $g_N$ that they experience. As a result, the inner and outer parts of galaxy rotation curves should fall on the same RAR. A difference is expected in modified gravity, essentially because the vertical gravity just outside the disc plane is higher nearer the galactic centre, reducing the MOND enhancement (Banik, Milgrom & Zhao 2018). Recently, a stacked analysis of SPARC rotation curves identified a $6.9\sigma$ difference between the RAR traced by data points from the inner and outer parts of galaxy rotation curves if the observational uncertainties are taken at face value (Chae 2022). The difference in the RAR evident in their fig. 6 is similar to that expected in modified gravity formulations of MOND. This could be a coincidence caused by unknown systematic errors, but this seems unlikely because almost the same signal is detected when excluding galaxies with a significant bulge or with a high luminosity, which could be signs of a greater degree of pressure support and thus possibly mean a larger uncertainty on the data points at low radii. This suggests that we should take seriously the observed difference in the RAR traced by the inner and outer parts of rotation curves. The high formal significance of the difference evident in their fig. 7 is a priori not expected in the modified inertia interpretation of MOND. One caveat is that gas motions are not perfectly circular, so until we have a fully fledged formulation of modified inertia theories capable of handling somewhat eccentric orbits, it is impossible to fully conclude on whether the results of Chae (2022) indeed rule out

such theories. But taking their results at face value and given also that the WBT rules out MOND as modified gravity (though see Milgrom 2023a, b), we most likely need to fundamentally change MOND if it is to survive the latest constraints, which questions its validity on all scales.

WBs are not the only bound systems that challenge MOND. Its predicted gravity in galaxy clusters generally falls short of the observed value, though one can assume a non-baryonic dark mass component that makes up about half of the gravitating mass (e.g. Sanders 2003; Angus, Famaey & Diaferio 2010). This discrepancy was recently illustrated in terms of the cluster RAR, as shown in fig. 8 of Eckert et al. (2022a) and fig. 5 of Li et al. (2023). The RAR traced by galaxy clusters lies about 0.3 dex (2$\times$) above the RAR traced by galaxies when $g_N \gtrsim 0.1 a_0$ (McGaugh 2020, and references therein). This gap narrows further and almost vanishes at the low acceleration end. The discrepancy with the Newtonian expectation without DM is of course larger at about 0.8 dex or 6$\times$, which matches the ratio between total and baryonic mass in $\Delta$CDM fits to the power spectrum of anisotropies in the cosmic microwave background radiation (Planck Collaboration VI 2020).

The WBT tests MOND on a much smaller scale than the traditional galaxy-scale tests with rotation curves, velocity dispersions, and weak lensing (section 3 of Banik & Zhao 2022, and references therein). This raises the possibility of suppressing MOND effects at short range, for instance with a length-dependent cutoff to any MOND enhancement to gravity below a scale of 0.1 pc. Babichev, Deffayet & Esposito-Farèse (2011) proposed an extended version of the Vainshtein screening mechanism in massive gravity (Vainshtein 1972; Babichev, Deffayet & Ziour 2010) in which MOND effects are suppressed below a length $r_B \propto M^{1/4}$ around a point mass $M$. Since $r_B$ rises slower than the MOND radius ($r_M \propto M^{1/2}$), we get that in a sufficiently massive system, the suppression of MOND effects arises only in the Newtonian regime and thus has no effect. But in a very low mass system, a new regime appears where the distance $r$ lies in the range $r_M < r < r_B$. In this regime, MOND as classically formulated predicts a departure from Newtonian gravity, but the behaviour becomes Newtonian with the extended screening mechanism (Babichev, Deffayet & Esposito-Farèse 2011).

Another way to think of this is in terms of the phantom dark matter (PDM) density $\rho_{pdm}$ generated by MOND, which is defined such that

$$ \nabla \cdot \mathbf{g} \equiv -4\pi G (\rho + \rho_{pdm}), \qquad (42) $$

where $g$ must be found by solving equation (12) and $\rho$ is the physical mass density. Since $g \propto \sqrt{M/r}$ in the deep-MOND limit [equation (1)], $\rho_{pdm} \propto \sqrt{M/r^2}$. If we focus on the PDM density at $r = r_B \propto M^{1/4}$, we get that $\rho_{pdm}$ is a constant. Thus, we can think of the model as providing a maximum limit to $|\rho_{pdm}|$ (the modulus is needed because the PDM density can be negative in more complicated geometries; see Milgrom 1986b; Oria et al. 2021). Such an upper limit arises in some attempts to unify the acceleration discrepancies in galaxies with the late-time accelerated expansion of the Universe (Zhao 2007). We can use the WBT to estimate an upper bound on $\rho_0$, the maximum possible $|\rho_{pdm}|$. Since our analysis indicates that $\alpha_{grav} < 0.2$ at rather high confidence (Table 3), we can assume that there is at most a 4 per cent enhancement to the orbital velocities of WBs in the asymptotic regime. This corresponds to an enclosed phantom mass equal to 8 per cent of the actual mass in the stars. Assuming our results are sensitive down to $M = 1 M_\odot$ and up to separations of $r = 20$ kAU or 0.1 pc, we can estimate that

$$ \rho_0 < 0.08 M_\odot \div \frac{4\pi r^3}{3} = 20 M_\odot/\text{pc}^3. \qquad (43) $$
This is three orders of magnitude above the local DM density inferred from a Newtonian dynamical analysis of Galactic data (Read 2014; Hagen & Helmi 2018; Salomon et al. 2020; de Salas & Widmark 2021). While we expect $ρ_{pdm}$ to be somewhat higher in a dwarf galaxy due to its lower mass, it is clear that the application of MOND to galaxies would be unaffected by a maximum limit to $|ρ_{pdm}|$ of order $20 M_{⊙}/\text{pc}^3$. The impact on larger scales would be even less significant, leaving open the issue of how MOND might be reconciled with observations of galaxy clusters and larger structures.

A maximum limit to $|ρ_{pdm}|$ capable of reconciling MOND with the WBT would substantially reduce MOND effects in the Solar System, invalidating attempts to explain some peculiarities of the Kuiper Belt using MOND (Brown & Mathur 2023; Migaszewski 2023) but improving the agreement with the rather tight bound on non-standard effects provided by Cassini radio tracking measurements of the Earth-Saturn range (Hees et al. 2014, 2016). This bound is already so tight that it formally rules out even the exponentially truncated MLS interpolating function at $8σ$ confidence (see equation C15 of Brown & Mathur 2023), so it could be difficult to reconcile Solar System ephemerides with galaxy dynamics in classical modified gravity theories of MOND (Desmond et al., in preparation). Since the PDM density at the MOND radius of an isolated point mass $M$ scales as $ρ_{pdm} ∝ 1/\sqrt{M}$ and we are roughly at the MOND radius of the MW with $M = 7 × 10^{10} M_{⊙}$ (Banik & Zhao 2018b), adding $ρ_0$ as a new constant of nature would suppress MOND effects only in systems up to six orders of magnitude less massive than the MW, i.e. up to about $M = 10^5 M_{⊙}$. This is in the regime of massive star clusters and globular clusters, which might help to explain why the internal kinematics of the outer halo globular cluster NGC 2419 are consistent with Newtonian expectations despite feeling little gravity from the Galaxy by virtue of its distance (Ibata et al. 2011a, b). While it has been argued that the observations are consistent with MOND as classically understood (Sanders 2012a, b), this requires one to invoke observational systematics like a radially varying polytropic equation of state. Since the mass of NGC 2419 is $9 × 10^5 M_{⊙}$, a clear detection or exclusion of MOND effects here would help to constrain how MOND effects must be suppressed on small scales to ensure consistency with the WBT. Interestingly, MOND correctly predicts the velocity dispersions of isolated LG dwarf galaxies down to $10^6 M_{⊙}$ (fig. 3 of McGaugh et al. 2021), while tidal stability considerations of Fornax Cluster dwarfs reach down to about the same mass (fig. 7 of Asencio et al. 2022). It is therefore possible that MOND as classically formulated breaks down at lower masses and that hints of this are already apparent in NGC 2419. A problem with this scenario is the asymmetric tidal tails of star clusters in the Solar neighbourhood, which have been argued to favour MOND (Kroupa et al. 2022; Pflamm-Altenburg et al. 2023). Further investigation is needed to see if the results can be explained in Newtonian gravity with a more complicated Galactic model that includes bars and spiral arms, though recent work indicates that the bar is not sufficient by itself for any pattern speed (Thomas et al. 2023). If MOND effects persist down to order $10^3 M_{⊙}$ but are not apparent in the WBT, then this would give a tight constraint on where classical MOND breaks down.

A limit to the PDM density might also limit the predicted Milgromian enhancement to the vertical gravity of the Galactic disc, which for now is neither confirmed nor rejected (Zhu et al. 2023). If it becomes clear that the enhancement is smaller than predicted, this would be interesting because the Galactic rotation curve works very well in MOND (McGaugh 2018; Zhu et al. 2023) and probes scales which are not much larger than the disc thickness, thereby

pinning down the boundary between where MOND works and where additional physics must be invoked to suppress MOND effects. While such tests would be useful, for MOND to be valid on kpc scales and ultimately provide the correct explanation for the RAR, we have to invoke at least two new fundamental constants ($a_0$ and $ρ_0$ or a fundamental length). This seems rather contrived and not in line with Occam's Razor, especially as $ρ_0$ is not obviously related to any cosmological density scale (unlike $a_0$) and cannot be substantially reduced below our estimate in equation (43) if we are to preserve the successes of MOND in dwarf galaxies (McGaugh et al. 2021; Asencio et al. 2022). We therefore discuss other theories that may be able to explain galaxy dynamics while passing the WBT.

Since we do not have a fully fledged quantum gravity theory, the idea of emergent gravity (EG; Verlinde 2017) is that gravity is an entropic force arising from underlying microscopic degrees of freedom that are currently not understood – and do not always need to be, much like how a very detailed understanding of molecules is not needed to understand the ideal gas law. EG faces major issues on galaxy cluster scales (Tamosiunas et al. 2019) and fails in galaxies because it predicts hook-shaped deviations from a single universal RAR that cannot be hidden within the uncertainties (Lelli, McGaugh & Schombert 2017a). While its prediction for the WBT is unclear, EG as presently formulated is strongly inconsistent with Solar System ephemerides (Hees, Famaey & Bertone 2017; Chan & Lee 2023). These precise constraints also rule out gravitational dipoles (Hajdukovic 2020; Banik & Kroupa 2020a) and scale-invariant dynamics (Maeder & Gueorguiev 2020; Banik & Kroupa 2020b). The WBT greatly exacerbates the already severe difficulties faced by modified gravity theories in explaining the RAR without DM while leaving no discernible trace in highly precise Solar System ephemerides, which our work effectively extends to 20 kAU at lower precision. We thus explore whether some form of DM on galaxy scales might provide a tight RAR (Famaey & McGaugh 2012) while lessening the very severe difficulties encountered by the standard CDM approach in other respects (Kroupa et al. 2010; Kroupa 2012, 2015; Di Valentino 2021; Banik & Zhao 2022, and references therein). We stress that the failures of $\Lambda$CDM identified in those works are unrelated to the validity of any alternatives proposed as solutions.

A hybrid MOND-DM model is provided by superfluid dark matter (SFDM; Berezhiani & Khoury 2015). Its basic idea is that galaxies are embedded in DM haloes, but their total mass within the region traced by rotation curves is not significant compared to the baryonic mass. Instead, the flat rotation curve problem is solved by postulating additional non-gravitational interactions between the baryons. These are mediated by phonons propagating in the DM halo, which is possible in its central superfluid core. On larger scales, SFDM reduces to $\Lambda$CDM as the superfluid phase only arises at low temperature and high density, so it arises in galaxies but not in galaxy clusters. SFDM should alleviate the problem faced by $\Lambda$CDM with the fast observed rotation speeds of galaxy bars (Roshan et al. 2021b) thanks to reduced dynamical friction on subsonic flows like a rotating galaxy bar (Berezhiani, Elder & Khoury 2019). However, the LG satellite planes are still very difficult to understand in SFDM because their anisotropy strongly suggests a tidal origin but their internal velocity dispersions imply some enhancement to the forces binding the satellites, which is hard to understand as the more distant satellites in these structures would be outside the superfluid portion of the halo (see section 5.6 of Roshan et al. 2021a). Another major issue is that the phonon-mediated forces at the heart of how SFDM reproduces the MOND phenomenology only enhance the forces on baryons, so it is difficult to understand why the RAR inferred
from rotation curves is also evident in strong lensing (Mistele, McGaugh & Hossenfelder 2022) and weak lensing data down to $g_N \approx 10^{-5} a_0$ (Brouwer et al. 2021) — the latter are particularly problematic for SFDM (Mistele, McGaugh & Hossenfelder 2023). These observational difficulties need to be considered alongside theoretical difficulties in making a stable covariant theory of DM superfluidity that reduces to MOND in galaxies (Hertzberg, Litterer & Shah 2021) and avoids significant orbital decay of stars due to Cherenkov-like radiation caused by their orbital velocity exceeding the local sound speed in the superfluid (Mistele 2022). If these difficulties are ultimately overcome, a positive aspect of SFDM is that the superfluid phase reduces to a normal phase near a star due to the steep potential gradient, avoiding any anomalous effects in the Solar System. However, this screening mechanism only works within a few hundred AU of a Sun-like star, beyond which MOND-like behaviour would be recovered (equation 86 of Berezhiani & Khoury 2015). Therefore, the WBT falsifies SFDM as presently understood. To pass the WBT, the force binding each WB should differ from the Newtonian expectation by $\lesssim 8$ per cent out to the separation limit of our WB sample. In the context of SFDM, this would require the normal phase bubble around every local Sun-like star to extend out to $\gtrsim 20$ kAU. This might be difficult to achieve by tuning the model parameters given constraints from galaxies and galaxy clusters (Hodson et al. 2017).

Thus, the extra forces needed to explain flat galaxy rotation curves might come from the mass of the DM rather than through it mediating a MOND-like interaction between the baryons. While the severe difficulties encountered by ΛCDM make this approach very unlikely (e.g. Banik & Zhao 2022, and references therein), the hypothetical DM particles might interact with each other — the standard assumption that they only interact gravitationally for all practical purposes is only made for simplicity. Self-interacting dark matter (SIDM; Spergel & Steinhardt 2000) offers a promising explanation for galaxies following a tight RAR (Ren et al. 2019). The required self-interaction cross-section is much larger than the upper limit imposed by the Newtonian dynamical mass profiles of galaxy clusters (Eckert et al. 2022b), but this could indicate that the cross-section depends strongly on the velocity. SIDM would limit the central DM density in galaxies, reducing dynamical friction on galaxy bars and thus possibly alleviating the $13\sigma$ fast bar tension faced by ΛCDM (Roshan et al. 2021b). SIDM would reduce to ΛCDM on large scales and thus still suffer the same issues with regards to the KBC void and Hubble tension (Haslbauer, Banik & Kroupa 2020) and the early formation of galaxies (Haslbauer et al. 2022) and galaxy clusters (Asencio, Banik & Kroupa 2021, 2023). In addition, SIDM does not provide an obvious explanation for the LG satellite planes because DM self-interactions do not have any obvious effect on the positions and velocities of satellite galaxies around their host — though the frequency of satellite planes in SIDM should be checked. Moreover, it is unclear how SIDM can explain the observed signs of tidal disturbance in Fornax Cluster dwarf galaxies and the lack of low surface brightness dwarfs towards the cluster centre (Asencio et al. 2022). Reducing the central DM density would make the situation less problematic, but the DM fraction within the baryonic extent of the Fornax dwarfs would need to be less than for isolated LG dwarfs with accurately measured internal velocity dispersions — these require a considerable amount of DM in a Newtonian context. Thus, SIDM models face some challenges but also appear to hold some promise, especially as they may explain some of the MOND phenomenology while being consistent with the WBT and Solar System ephemerides due to the lack of any change in the gravity sector. The same is true if we postulate a long-range non-

gravitational interaction between baryons and dark matter beyond their self-interactions (Famaey et al. 2020).

## 5.5 Future prospects

The significant challenge to MOND presented by the WBT is an important result which should be confirmed with additional data. There are good prospects for future improvement: the next Gaia data release (DR4) is anticipated at the end of 2025. DR4 will double the time baseline analysed in DR3, giving an $\approx 2.8 \times$ improvement in proper motion precision and reducing the impact of CBs at AU separations (Appendix B). We estimate that using DR4 with a similar selection to our main analysis but extended to $m_G < 17.6$ and $d < 400$ pc would approximately double the search volume and almost triple the final sample size, bearing in mind that more precise proper motions will lead to more systems passing our cut on the $\tilde{v}$ precision [equation 10].

Systems rejected here due to a lack of either star's RV can also be recovered from ground-based RV observations at a moderate precision of $\approx 1$ km s$^{-1}$. This would increase the sample size by about 10 per cent, but the more useful gain could be that the contamination rate can be reduced if the RV is known for both stars.

There are also interesting prospects for directly constraining the triple population to reduce the degrees of freedom inherent in our modelling (Manchanda, Sutherland & Pittordis 2023). Their simulations have shown that we should be able to detect nearly all triples with a single main sequence star orbited by a CB whose inner period $\gtrsim 3$ years, which is long enough to give a substantial time-averaged perturbation to the photocentre velocity (Appendix B). Such hierarchical systems can be found either by Gaia astrometric accelerations in the 10-yr final data set (for separations $\lesssim 25$ AU); or by direct, speckle, or coronagraphic imaging (for separations $\gtrsim 20$ AU). The required imaging observations are much less time-consuming than high-precision RV measurements, so they may be feasible for a substantial fraction of the WB sample, especially for the high $\tilde{v}$ tail. This would provide an important check on our estimated CB contamination fraction.

In principle, it is possible to refine the WBT via the addition of binary RV differences. This would allow the use of 3D velocities and 2D projected separations, which would reduce the scatter in the $\tilde{v}$ distribution inherent from projected velocities. However, this would require high precision (order $0.05$ km s$^{-1}$) RVs for all stars, which would be very costly in observing time and require corrections for gravitational redshift (Loeb 2022) and convective blueshift (Liebing et al. 2021). A fully 3D version of the WBT would also require much more precise astrometry to obtain kAU-level constraints on the relative heliocentric distances to the stars in each WB. This is not envisaged in the foreseeable future for a statistically large sample of WBs. More precise information is available for some very nearby systems (Kervella et al. 2016; Kervella, Thévenin & Lovis 2017), but an individual WB cannot test MOND using velocities alone even if $\tilde{v} > \sqrt{2}$ because a single snapshot would not prove that the WB is bound (though a chance flyby is unlikely; see section 8.1 of BZ18). A test is possible if the acceleration can be directly measured, which may become feasible with future observatories (Banik & Kroupa 2019a).

Our work has assumed that the potential of a star at kAU distances must be traced by another star. It is possible for a spacecraft to serve as the tracer instead (Banik & Kroupa 2019b). The predicted aspherical shell of PDM concentrated around the Solar MOND radius would also affect Solar System ephemerides (Hees et al. 2014, 2016; Brown & Mathur 2023).
# 6 SUMMARY AND CONCLUSIONS

MOND has enjoyed unparalleled predictive success with regards to galaxy dynamics (Famaey & McGaugh 2012; Banik & Zhao 2022, and references therein). Its central postulate is that the dynamics of a system deviates from Newtonian expectations when the gravity $g \lesssim a_0$. This could in principle occur at any distance from a system provided it has a sufficiently low mass [equation (2)]. Since the MOND radius of a star is $7\sqrt{M/M_{\odot}}$ kAU, MOND ought to have detectable effects on local WBs with kAU separations. Indeed, simple analytic estimates backed up by detailed calculations show a 20 per cent enhancement to the orbital velocity over the Newtonian expectation (BZ18), as also evident in N-body simulations of local star clusters for much the same reason (Kroupa et al. 2022).

We test this prediction using an observational sample of 8611 WBs from Gaia DR3 (Gaia Collaboration 2023) where the uncertainty on $\tilde{v}$ [equation (3)] is very small and other quality cuts have been applied (Sections 2.1 and 2.4). Dividing our sample into ten equally sized subsamples shows no trend in the median $\tilde{v}$ with respect to $r_{sky}/r_M$ once we restrict our attention to $\tilde{v} \lesssim 2$, which should significantly reduce contaminating effects while preserving a genuine MOND signal (Fig. 11). Observational uncertainties would be expected to broaden the $\tilde{v}$ distribution preferentially at large separations because the same uncertainty on the relative velocity implies a larger uncertainty in $\tilde{v}$. This would if anything lead to a broader $\tilde{v}$ distribution at low accelerations. As a result, the flat behaviour of the median $\tilde{v}$ with respect to our proxy for the internal WB acceleration strongly suggests that local WBs are Newtonian.

We then present the most detailed statistical hypothesis test of MOND to date using local WBs. Our model includes a rigorous calculation of the MOND gravitational field under an EFE of strength $g_e = 1.8 a_0$ (Section 3.1), an allowance for an undetected CB companion to one or both of the stars in each WB (Section 3.2), and LOS contamination (Section 3.3). The procedure was fixed in advance as much as possible to minimise biases that could arise when conducting the WBT (Banik, Pittordis & Sutherland 2021). The model parameters relate to those of:

(i) The WBs, for which we need the $\alpha_{grav}$ parameter interpolating between different gravity laws [see equation (39) for its relation to the gravity law] and the distribution of semi-major axes [$a_{break}$ and $\beta$; see equation (16)] and orbital eccentricities [$\gamma$; see equation (17)];

(ii) The CBs, for which we assume the same $\gamma$ parameter but also need the fraction $f_{CB}$ of stars in our WB sample with an undetected CB companion (Section 3.2.4) and the upper limit to the CB semi-major axis, which is defined by $k_{CB}$ [see equation (22)]; and

(iii) The LOS contamination fraction $f_{LOS}$ due to chance alignments (see Section 3.3).

These seven parameters are allowed to vary freely so as to best match the observed distribution of $(r_{sky}, \tilde{v})$ in 540 pixels (Table 1).

Our best-fitting Newtonian model significantly outperforms our best-fitting MOND model: the likelihood ratio of exp(175) implies a preference for Newtonian gravity at $19\sigma$ confidence according to our analyses with fixed $\alpha_{grav}$. This is in line with the result of PS23, who conducted a less detailed version of the WBT using a somewhat different $r_{sky}$ range of 5–20 kAU and considered WBs with $\tilde{v} < 7$. Their analysis also differs from ours in several other respects, including their use of a linear rather than cubic mass-luminosity relation (Fig. 1), their CB model, and their much more limited exploration of the parameter space. Those authors found that the $\chi^2$ of their best Newtonian model is smaller than that of their best MOND model by 525, suggesting a $23\sigma$ preference for

Newtonian gravity. While our best-fitting Newtonian model is not a perfect representation of the WB data set, it is clear that the data do not show the predicted broadening to the $\tilde{v}$ distribution at larger separations (Fig. 12).

In our main analysis, we allow the gravity law to interpolate between Newtonian and Milgromian. For this, we use the gravity law parameter $\alpha_{grav}$, which is 0 in Newtonian gravity and 1 in MOND. Using $10^5$ MCMC trials to explore the parameter space, we find that the inferred $\alpha_{grav} = -0.021_{-0.045}^{+0.065}$. This is very consistent with Newtonian dynamics but rules out MOND at $16\sigma$ confidence, which is in line with the results of PS23 and our results for fixed $\alpha_{grav}$. Our result is robust to various changes in the modelling assumptions and the sample selection (Table 3), including when we use a narrower mass range to minimize possible trends in the mass distribution across the parameter range used for the WBT (Section 5.1.1). We also explain in some detail why, without a drastic change to our understanding of baryonic surface densities, the MOND interpolating function cannot be chosen to simultaneously pass the WBT and constraints from the rotation curves of galaxies, including our own (Section 5.3).

Our conclusion disagrees with two studies that were published while this work was under review (Chae 2023; Hernandez 2023). In Section 5.2, we focus on the study of Chae (2023) due to its much larger sample size and its claim to have detected the predicted MOND enhancement to Newtonian gravity at $10\sigma$ confidence. We find that the nominal sample of 26 615 WBs used in that study does indeed show a clear signal that closely resembles the MOND expectation and appears to rule out Newtonian dynamics (Fig. 19). However, we then identify a major deficiency with the handling of astrometric uncertainties. Despite the importance of the relative velocity between the stars in each WB, the uncertainty on this quantity is not estimated in Chae (2023). Instead, the focus is on ensuring that the heliocentric velocity of each star is very precise. This is insufficient to ensure a precisely known relative velocity and thus a reliable $\tilde{v}$, which however is essential to conducting the WBT. We therefore use equation 5 of El-Badry, Rix & Heintz (2021) to estimate the uncertainty in $v_{sky}$ and thus in $\tilde{v}$. We then impose an additional condition on the $\tilde{v}$ uncertainty matching that used in our own analysis [equation (10)]. This completely removes the apparent MOND signal despite only reducing the sample size by about 1/5 (see Fig. 21). We argue that a similar issue may well have affected the analysis of Hernandez (2023) based on some statements in that paper, though we do not study its WB sample in detail. We also find that the WB sample used by Chae (2023) artificially imposes that $\tilde{v} < \sqrt{5} M_{\odot}/M$, making it not very well suited to the WBT (Appendix E).

We conclude that the gravity law inferred from our analysis of local WBs is consistent with Newtonian expectations but rules out MOND as modified gravity at $\gg 5\sigma$ confidence in both our nominal analysis and a considerable range of variations to it. This conclusion is of course reliant on our modelling approach, which is not a perfect match to the data (Fig. D1).<fnref n="19" /> Even so, neither our main analysis nor any of the variations considered can significantly improve the fit by using a gravity law different to Newtonian. This is also evident from the model-independent Fig. 11, which considers the median $\tilde{v}$ in subsamples with different $r_{sky}/r_M$ as a proxy for the WB internal

<footnote n="19">¹⁹This is similar to the rotation curve predictions in MOND, which often do not match the data within formal uncertainties but visually provide a good fit, which is typically considered sufficient given the inevitable modelling deficiencies (Kroupa et al. 2018; Cameron, Angus & Burgess 2020).</footnote>
acceleration. It has been argued that modified inertia interpretations of MOND could reduce the predicted enhancement to the Newtonian acceleration and thus remain consistent with the WBT (Milgrom 1994, 2011, 2022). This would lead to the inner and outer parts of galaxy rotation curves following the same RAR, but a difference is detected at $6.9\sigma$ confidence (Chae 2022). Although the predictions of modified inertia for slightly non-circular orbits are not known, the magnitude of the above-mentioned difference is consistent with expectations of MOND as modified gravity – but this is ruled out by the WBT. Further modifications to MOND could preserve its successes on galaxy scales, especially with regards to issues like the LG satellite planes, which are difficult to understand any other way (Pawlowski 2021a, b; Banik et al. 2022). This entails at least one new fundamental scale beyond $a_0$, marking the end of MOND as a purely acceleration-dependent modification to standard physics. Limiting the MOND phantom density to $|\rho_{pdm}| \lesssim 20 M_\odot/\text{pc}^3$ would yield consistency with the WBT and also alleviate tensions related to the null detection of MOND effects in Cassini radio tracking data from Saturn (Hees et al. 2014, 2016; Brown & Mathur 2023). It might also improve the agreement with observations of globular clusters and the vertical gravity from the Galactic disc. Limiting $|\rho_{pdm}|$ would have little effect on larger scales, similarly to the Vainshtein mechanism used to screen modified gravity effects in the Solar System (Babichev, Deffayet & Esposito-Farèse 2011). This would preserve MOND's successes in galaxies but leave open the issue of galaxy clusters and large-scale structure, for which the hybrid vHDM paradigm may be a promising approach (section 9.2 of Banik & Zhao 2022, and references therein). While our results falsify MOND as currently understood, given the many problems for $\Lambda$CDM discussed in that work, our results cannot be used to argue that it is the correct model either – both models are clearly incomplete. Hybrid models like SFDM struggle to explain galactic-scale observations like lensing (Mistele, McGaugh & Hossenfelder 2022, 2023) and the WBT, while non-MOND modifications to gravity like EG and MOG usually fail in galaxies (Section 5.4). SIDM may be a promising approach but it is very similar to $\Lambda$CDM on large scales, thus encountering the same difficulties with the KBC void (Keenan, Barger & Cowie 2013; Haslbauer, Banik & Kroupa 2020; Wong et al. 2022) and Hubble tension (Di Valentino 2021). This significant anomaly for standard cosmology appears to persist in the JWST era (Yuan et al. 2022) and must be solved consistently with the ages of the oldest stars (Cimatti & Moresco 2023, and references therein). We hope that our results from local WBs motivate the development of a more complete theory, which is likely to borrow some elements from both $\Lambda$CDM and MOND given their successes in different domains.

## ACKNOWLEDGEMENTS

IB is supported by Science and Technology Facilities Council grant ST/V000861/1, which also partially supports HZ. IB acknowledges support from a “Pathways to Research” fellowship from the University of Bonn, during which the primary statistical analysis was largely coded. BF and RI acknowledge funding from the Agence Nationale de la Recherche (ANR projects ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European Research Council (ERC) under the European Union’s Horizon 2020 Framework programme (grant agreement number 834148).

The authors thank the referee for highlighting some important deficiencies with an earlier version of this work regarding the modelling of undetected companions. They also thank Pavel Kroupa for comments on the article and for lending computational resources from the Stellar Populations and Dynamics Research Group at the

University of Bonn, where all the statistical analyses were conducted. IB is very grateful to Elena Asencio for helping to check every line of the statistical analysis code and for lending algorithms used for making several figures in this paper, especially the triangle plots made using PYGTC (Bocquet & Carter 2016). He thanks Haixia Ma for providing Fig. 22. The authors are grateful for a masters project by Nathan Findlay which gave a much better understanding of how to analyse the data set and what quality cuts would be appropriate. IB acknowledges significant computational resources provided by Kyu-Hyun Chae to simulate the impact of undetected companions. Most of the figures were prepared using MATLAB.

## DATA AVAILABILITY

The WB data set used in this contribution will be made available upon reasonable request, showing Gaia catalogue information for the sample of 19786 WBs described in the text and an extra flag to indicate if each WB is one of the 8611 systems used in our main analysis. The WB data set used in the nominal analysis of Chae (2023) was provided by Kyu-Hyun Chae. Gaia DR3 is publicly available.

## REFERENCES

* Acedo L., 2020, Universe, 6, 209

* Ambartsumian V. A., 1937, Astron. Zh., 14, 207

* Andrews J. J., Chanamé J., Agüeros M. A., 2017, MNRAS, 472, 675

* Angus G. W., 2009, MNRAS, 394, 527

* Angus G. W., Famaey B., Diaferio A., 2010, MNRAS, 402, 395

* Arenou F. et al., 2018, A&A, 616, A17

* Asencio E., Banik I., Kroupa P., 2021, MNRAS, 500, 5249

* Asencio E., Banik I., Mieske S., Venhola A., Kroupa P., Zhao H., 2022, MNRAS, 515, 2981

* Asencio E., Banik I., Kroupa P., 2023, ApJ, 954, 162

* Babichev E., Deffayet C., Ziour R., 2010, Phys. Rev. D, 82, 104008

* Babichev E., Deffayet C., Esposito-Farèse G., 2011, Phys. Rev. D, 84, 061502

* Bahcall J. N., Hut P., Tremaine S., 1985, ApJ, 290, 15

* Banik I., 2019, MNRAS, 487, 5291

* Banik I., Kroupa P., 2019a, MNRAS, 487, 1653

* Banik I., Kroupa P., 2019b, MNRAS, 487, 2665

* Banik I., Kroupa P., 2020a, MNRAS, 495, 3974

* Banik I., Kroupa P., 2020b, MNRAS, 497, L62

* Banik I., Zhao H., 2018a, SciFed J. Astrophys., 1, 1000008

* Banik I., Zhao H., 2018b, MNRAS, 473, 419

* Banik I., Zhao H., 2018c, MNRAS, 480, 2660

* Banik I., Zhao H., 2022, Symmetry, 14, 1331

* Banik I., Milgrom M., Zhao H., 2018, preprint (arXiv:1808.10545)

* Banik I., Thies I., Candlish G., Famaey B., Ibata R., Kroupa P., 2020, ApJ, 905, 135

* Banik I., Pittordis C., Sutherland W., 2021, preprint (arXiv:2109.03827)

* Banik I., Thies I., Truelove R., Candlish G., Famaey B., Pawlowski M. S., Ibata R., Kroupa P., 2022, MNRAS, 513, 129

* Beech M., 2009, MNRAS, 399, L21

* Beech M., 2011, Ap&SS, 333, 419

* Begeman K. G., Broeils A. H., Sanders R. H., 1991, MNRAS, 249, 523

* Bekenstein J., Milgrom M., 1984, ApJ, 286, 7

* Belokurov V. et al., 2020, MNRAS, 496, 1922

* Berezhiani L., Khoury J., 2015, Phys. Rev. D, 92, 103510

* Berezhiani L., Elder B., Khoury J., 2019, J. Cosmol. Astropart. Phys., 2019, 074

* Bflek M., Thies I., Kroupa P., Famaey B., 2018, A&A, 614, A59

* Bocquet S., Carter F. W., 2016, J. Open Source Softw., 1, 46

* Bosma A., 1978, PhD thesis, Groningen University

* Bovy J., Bird J. C., García Pérez A. E., Majewski S. R., Nidever D. L., Zasowski G., 2015, ApJ, 800, 83

* Brouwer M. M. et al., 2021, A&A, 650, A113
Brown K., Mathur H., 2023, AJ, 166, 168

Cameron E., Angus G. W., Burgess J. M., 2020, Nat. Astron., 4, 132

Candlish G. N., Smith R., Jaffé Y., Cortesi A., 2018, MNRAS, 480, 5362

Chabrier G., Baraffe I., 2000, ARA&A, 38, 337

Chae K.-H., 2022, ApJ, 941, 55

Chae K.-H., 2023, ApJ, 952, 128

Chae K.-H., 2024, ApJ, in press (arXiv:2309.10404)

Chae K.-H., Gong I.-T., 2015, MNRAS, 451, 1719

Chae K.-H., Milgrom M., 2022, ApJ, 928, 24

Chae K.-H., Bernardi M., Sheth R. K., Gong I.-T., 2019, ApJ, 877, 18

Chae K.-H., Bernardi M., Domínguez Sánchez H., Sheth R. K., 2020a, ApJ, 903, L31

Chae K.-H., Lelli F., Desmond H., McGaugh S. S., Li P., Schombert J. M., 2020b, ApJ, 904, 51

Chae K.-H., Desmond H., Lelli F., McGaugh S. S., Schombert J. M., 2021, ApJ, 921, 104

Chan M. H., Lee C. M., 2023, MNRAS, 518, 6238

Chrobáková Ž., López-Corredoira M., Sylos Labini F., Wang H. F., Nagy R., 2020, A&A, 642, A95

Cimatti A., Moresco M., 2023, ApJ, 953, 149

Ciotti L., 2022, ApJ, 936, 180

Clarke C. J., 2020, MNRAS, 491, L72

Desmond H., 2023, MNRAS, 526, 3342

de Almeida Á. O. F., Piattella O. F., Rodrigues D. C., 2016, MNRAS, 462, 2706

de Salas P. F., Widmark A., 2021, Rep. Prog. Phys., 84, 104901

Di Valentino E., 2021, MNRAS, 502, 2065

Dinnbier F., Kroupa P., Anderson R. I., 2022, A&A, 660, A61

Eckert D., Ettori S., Pointecouteau E., van der Burg R. F. J., Loubser S. I., 2022a, A&A, 662, A123

Eckert D., Ettori S., Robertson A., Massey R., Pointecouteau E., Harvey D., McCarthy I. G., 2022b, A&A, 666, A41

Efstathiou G., Sutherland W. J., Maddox S. J., 1990, Nature, 348, 705

Eilers A.-C., Hogg D. W., Rix H.-W., Ness M. K., 2019, ApJ, 871, 120

Einstein A., 1915, Sitzungsber. Kgl Preuss. Akad. Wiss. (Berlin), 6, p. 844

Eker Z. et al., 2015, AJ, 149, 131

El-Badry K., 2019, MNRAS, 482, 5018

El-Badry K., Rix H.-W., 2018, MNRAS, 480, 4884

El-Badry K., Rix H.-W., Tian H., Duchêne G., Moe M., 2019, MNRAS, 489, 5822

El-Badry K., Rix H.-W., Heintz T. M., 2021, MNRAS, 506, 2269

Faber S. M., Gallagher J. S., 1979, ARA&A, 17, 135

Famaey B., Binney J., 2005, MNRAS, 363, 603

Famaey B., McGaugh S. S., 2012, Living Rev. Relativ., 15, 10

Famaey B., Bruneton J.-P., Zhao H., 2007, MNRAS, 377, L79

Famaey B., Khoury J., Penco R., Sharma A., 2020, J. Cosmol. Astropart. Phys., 2020, 025

Feng F., Jones H. R. A., 2018, MNRAS, 473, 3185

Fletcher R., Powell M. J. D., 1963, Comput. J., 6, 163

Francis C., Anderson E., 2014, Celest. Mech. Dyn. Astron., 118, 399

Gaia Collaboration, 2016, A&A, 595, A1

Gaia Collaboration, 2018, A&A, 616, A1

Gaia Collaboration, 2021a, A&A, 649, A1

Gaia Collaboration, 2021b, A&A, 649, A9

Gaia Collaboration, 2023, A&A, 674, A1

Gelman A., Gilks W. R., Roberts G. O., 1997, Ann. Appl. Probab., 7, 110

Gentile G., Famaey B., de Blok W. J. G., 2011, A&A, 527, A76

Gordon D., de Witt A., Jacobs C. S., 2023, AJ, 165, 49

Gravity Collaboration, 2019, A&A, 625, L10

Hagen J. H. J., Helmi A., 2018, A&A, 615, A99

Haghi H., Amiri V., Hasani Zonoozi A., Banik I., Kroupa P., Haslbauer M., 2019, ApJ, 884, L25

Hajdukovic D. S., 2020, MNRAS, 491, 4816

Hartman Z. D., Lépine S., Medan I., 2022, ApJ, 934, 72

Haslbauer M., Banik I., Kroupa P., 2020, MNRAS, 499, 2845

Haslbauer M., Kroupa P., Zonoozi A. H., Haghi H., 2022, ApJ, 939, L31

Hayashi C., 1961, PASJ, 13, 450

Hayashi C., Hoshi R., 1961, PASJ, 13, 442

Hees A., Folkner W. M., Jacobson R. A., Park R. S., 2014, Phys. Rev. D, 89, 102002

Hees A., Famaey B., Angus G. W., Gentile G., 2016, MNRAS, 455, 449

Hees A., Famaey B., Bertone G., 2017, Phys. Rev. D, 95, 064019

Hejlesen P. M., 1980, A&AS, 39, 347

Hernandez X., 2023, MNRAS, 525, 1401

Hernandez X., Jiménez M. A., Allen C., 2012, Eur. Phys. J. C, 72, 1884

Hernandez X., Cortés R. A. M., Allen C., Scarpa R., 2019, Int. J. Mod. Phys. D, 28, 1950101

Hernandez X., Cookson S., Cortés R. A. M., 2022, MNRAS, 509, 2304

Hernandez X., Verteletskyi V., Nasser L., Aguayo-Ortiz A., 2024, MNRAS, in press (arXiv:2309.10995)

Hertzberg M. P., Litterer J. A., Shah N., 2021, J. Cosmol. Astropart. Phys., 2021, 015

Hodson A. O., Zhao H., Khoury J., Famaey B., 2017, A&A, 607, A108

Hoof S., Geringer-Sameth A., Trotta R., 2020, J. Cosmol. Astropart. Phys., 2020, 012

Hwang H.-C., Ting Y.-S., Zakamska N. L., 2022, MNRAS, 512, 3383

Ibata R., Sollima A., Nipoti C., Bellazzini M., Chapman S. C., Dalessandro E., 2011a, ApJ, 738, 186

Ibata R., Sollima A., Nipoti C., Bellazzini M., Chapman S. C., Dalessandro E., 2011b, ApJ, 743, 43

Jiao W.-C., Henry T. J., Gies D. R., Hambly N. C., 2018, ApJ, 861, L11

Jeans J. H., 1919, MNRAS, 79, 408

Jiang Y.-F., Tremaine S., 2010, MNRAS, 401, 977

Jones-Smith K., Abraham R., Kell L., Mathur H., 2018, New J. Phys., 20, 063042

Kashfi T., Roshan M., Famaey B., 2023, MNRAS, 519, L21

Keenan R. C., Barger A. J., Cowie L. L., 2013, ApJ, 775, 62

Kervella P., Mignard F., Mérand A., Thévenin F., 2016, A&A, 594, A107

Kervella P., Thévenin F., Lovis C., 2017, A&A, 598, L7

Knox R. A., Hawkins M. R. S., Hambly N. C., 1999, MNRAS, 306, 736

Kortnevčič C., Kaczmarek T., Pfalzner S., 2012, A&A, 543, A126

Kroupa P., 2008, Initial Conditions for Star Clusters. Lecture Notes in Physics Series, Springer, Dordrecht, The Netherlands, p. 181

Kroupa P., 2012, Publ. Astron. Soc. Aust., 29, 395

Kroupa P., 2015, Can. J. Phys., 93, 169

Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545

Kroupa P. et al., 2010, A&A, 523, A32

Kroupa P. et al., 2018, Nat. Astron., 2, 925

Kroupa P. et al., 2022, MNRAS, 517, 3613

LUX Collaboration, 2017, Phys. Rev. Lett., 118, 021303

Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157

Lelli F., McGaugh S. S., Schombert J. M., 2017a, MNRAS, 468, L68

Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017b, ApJ, 836, 152

Lépine S., Bongiorno B., 2007, AJ, 133, 889

Li P., Lelli F., McGaugh S., Schombert J., 2018, A&A, 615, A3

Li P. et al., 2023, A&A, 677, A24

Liebing F., Jeffers S. V., Reiners A., Zechmeister M., 2021, A&A, 654, A168

Lindegren L. et al., 2021, A&A, 649, A4

Loeb A., 2022, Res. Notes AAS, 6, 55

Lux-Zeplin Collaboration, 2023, Phys. Rev. Lett., 131, 041002

Maeder A., Gueorguiev V. G., 2020, MNRAS, 492, 2698

Manchanda D., Sutherland W., Pittordis C., 2023, Open J. Astrophys., 6, E2

Mann A. W., Feiden G. A., Gaidos E., Boyajian T., von Braun K., 2015, ApJ, 804, 64

Mansfield S., Kroupa P., 2021, A&A, 650, A184

McConnachie A. W., 2012, AJ, 144, 4

McCulloch M. E., Lucio J. H., 2019, Ap&SS, 364, 121

McGaugh S. S., 2008, ApJ, 683, 137

McGaugh S. S., 2012, AJ, 143, 40

McGaugh S. S., 2018, Res. Notes AAS, 2, 156

McGaugh S., 2020, Galaxies, 8, 35

McGaugh S. S., de Blok W. J. G., 1998a, ApJ, 499, 66

McGaugh S. S., de Blok W. J. G., 1998b, ApJ, 499, 41

McGaugh S. S., Schombert J. M., Bothun G. D., 1995, AJ, 109, 2019

McGaugh S., Lelli F., Schombert J., 2016, Phys. Rev. Lett., 117, 201101
McGaugh S. S., Lelli F., Schombert J. M., Li P., Visgaitis T., Parker K. S.,
Pawlowski M. S., 2021, AJ, 162, 202

McMillan P. J., 2017, MNRAS, 465, 76

Migaszewski C., 2023, MNRAS, 525, 805

Milgrom M., 1983, ApJ, 270, 365

Milgrom M., 1986a, ApJ, 302, 617

Milgrom M., 1986b, ApJ, 306, 9

Milgrom M., 1994, Ann. Phys., 229, 384

Milgrom M., 2010, MNRAS, 403, 886

Milgrom M., 2011, Acta Phys. Pol. B, 42, 2175

Milgrom M., 2022, Phys. Rev. D, 106, 064060

Milgrom M., 2023a, Phys. Rev. D, 108, 063009

Milgrom M., 2023b, Phys. Rev. D, 108, 084005

Mistele T., 2022, J. Cosmol. Astropart. Phys., 2022, 008

Mistele T., McGaugh S., Hossenfelder S., 2022, A&A, 664, A40

Mistele T., McGaugh S., Hossenfelder S., 2023, J. Cosmol. Astropart. Phys.,
2023, 004

Moffat J. W., 2006, J. Cosmol. Astropart. Phys., 2006, 004

Mróz P. et al., 2019, ApJ, 870, L10

Nagesh S. T. et al., 2023, MNRAS, 519, 5128

Negrelli C., Benito M., Landau S., Iocco F., Kraiselburd L., 2018, Phys. Rev. D, 98, 104061

Nelson T., Ting Y.-S., Hawkins K., Ji A., Kamdar H., El-Badry K., 2021, ApJ, 921, 118

Offner S. S. R., Moe M., Kratter K. M., Sadavoy S. I., Jensen E. L. N., Tobin J. J., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M. eds, ASP Conf. Ser. 534, Protostars and Planets VII, Astron. Soc. Pac., San Francisco, p. 275

Oh S., Price-Whelan A. M., Hogg D. W., Morton T. D., Spergel D. N., 2017, AJ, 153, 257

Öpik E., 1924, Publ. Tartu Astros. Obs., 25, 1

Oria P. A. et al., 2021, ApJ, 923, 68

Ostriker J. P., Steinhardt P. J., 1995, Nature, 377, 600

Pace A. B., Erkal D., Li T. S., 2022, ApJ, 940, 136

Paučo R., Klačka J., 2016, A&A, 589, A63

Pawlowski M. S., 2021a, Nat. Astron., 5, 1185

Pawlowski M. S., 2021b, Galaxies, 9, 66

Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9

Peebles P. J. E., 2017, Nat. Astron., 1, 0057

Penner A. R., 2020, Ap&SS, 365, 65

Penoyre Z., Belokurov V., Evans N. W., 2022, MNRAS, 513, 2437

Petersen J., Lelli F., 2020, A&A, 636, A56

Pflamm-Altenburg J., Kroupa P., Thies I., Jerabkova T., Beccari G., Prusti T.,
Boffin H. M. J., 2023, A&A, 671, A88

Pichon B., 2007, in Bouvier J., Chalabaev A., Charbonnel C. eds, SF2A-2007: Proceedings of the Annual Meeting of the French Society of Astronomy and Astrophysics. Grenoble, France, Société Francaise d'Astronomie et d'Astrophysique, p. 549

Pittordis C., Sutherland W., 2018, MNRAS, 480, 1778

Pittordis C., Sutherland W., 2019, MNRAS, 488, 4740

Pittordis C., Sutherland W., 2023, Open J. Astrophys., 6, 4

Planck Collaboration VI, 2020, A&A, 641, A6

Read J. I., 2014, J. Phys. G Nucl. Phys., 41, 063101

Reid M. J., Brunthaler A., 2004, ApJ, 616, 872

Reid I. N., Gizis J. E., Hawley S. L., 2002, AJ, 124, 2721

Ren T., Kwa A., Kaplinghat M., Yu H.-B., 2019, Phys. Rev. X, 9, 031020

Renaud F., Famaey B., Kroupa P., 2016, MNRAS, 463, 3637

Riddle R. L. et al., 2015, ApJ, 799, 4

Riello M. et al., 2021, A&A, 649, A3

Roshan M., Banik I., Ghafourian N., Thies I., Famaey B., Asencio E., Kroupa P., 2021a, MNRAS, 503, 2833

Roshan M., Ghafourian N., Kashfi T., Banik I., Haslbauer M., Cuomo V., Famaey B., Kroupa P., 2021b, MNRAS, 508, 926

Rowland D. R., 2015, Int. J. Mod. Phys. D, 24, 1550065

Rubin V. C., Ford W. K. J., Thonnard N., 1978, ApJ, 225, L107

Salomon J.-B., Bienaymé O., Reylé C., Robin A. C., Famaey B., 2020, A&A, 643, A75

Sanders R. H., 2003, MNRAS, 342, 901

Sanders R. H., 2012a, MNRAS, 419, L6

Sanders R. H., 2012b, MNRAS, 422, L21

Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525

Schombert J., McGaugh S., Lelli F., 2022, AJ, 163, 154

Schönrich R., 2012, MNRAS, 427, 274

Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829

Sellwood J. A., Shen J., Li Z., 2019, MNRAS, 486, 4710

Shaya E. J., Olling R. P., 2011, ApJS, 192, 2

Skordis C., Złośnik T., 2019, Phys. Rev. D, 100, 104013

Skordis C., Złośnik T., 2021, Phys. Rev. Lett., 127, 161302

Skordis C., Złośnik T., 2022, Phys. Rev. D, 106, 104041

Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760

Tamosiunas A., Bacon D., Koyama K., Nichol R. C., 2019, J. Cosmol. Astropart. Phys., 2019, 053

Thomas G. F., Famaey B., Ibata R., Renaud F., Martin N. F., Kroupa P., 2018, A&A, 609, A44

Thomas G. F., Famaey B., Monari G., Laporte C. F. P., Ibata R., de Laverny P., Hill V., Boily C., 2023, A&A, 678, A180

Tian Y., Ko C.-M., 2019, MNRAS, 488, L41

Tokovinin A., 2017, MNRAS, 468, 3461

Tokovinin A., 2023, AJ, 165, 165

Tokovinin A., Moe M., 2020, MNRAS, 491, 5158

Tokovinin A. A., Smekhov M. G., 2002, A&A, 382, 118

Tokovinin A., Hartung M., Hayward T. L., 2010, AJ, 140, 510

Vainshtein A. I., 1972, Phys. Lett. B, 39, 393

Verlinde E., 2017, SciPost Phys., 2, 016

Wang J., Hammer F., Yang Y., 2022, MNRAS, 510, 2242

Wong J. H. W., Shanks T., Metcalfe N., Whitbourn J. R., 2022, MNRAS, 511, 5742

Yuan W., Riess A. G., Casertano S., Macri L. M., 2022, ApJ, 940, L17

Zhao H., 2007, ApJ, 671, L1

Zhao H., Li B., Bienaymé O., 2010, Phys. Rev. D, 82, 103001

Zhao H., Famaey B., Lüghausen F., Kroupa P., 2013, A&A, 557, L3

Zhou Y., Li X., Huang Y., Zhang H., 2023, ApJ, 946, 73

Zhu Y. et al., 2023, MNRAS, 519, 4479

Zonoozi A. H., Lieberz P., Banik I., Haghi H., Kroupa P., 2021, MNRAS, 506, 5468

APPENDIX A: EXACTLY EQUAL MASS BINARIES

Although we might expect the mass ratio distribution of binaries to be smooth, observations indicate a population of exactly equal mass WBs (El-Badry et al. 2019). We therefore include an allowance for a δ-function in the $\bar{q}$ distribution at 1/2, where $\bar{q}$ is the fraction of the CB mass in the less massive star. To find the likelihood $P_{eqm}$ of a binary having exactly equal mass stars, we plot the cumulative distribution of $0.5 - \bar{q}$ for the WBs in our sample, which we restrict to $\bar{v} < 1$ to try and mitigate CB contamination. The idea is to extrapolate the distribution down to zero and look for a positive intercept. Since a δ-function would be smeared somewhat by observational uncertainties, we need to start our fit at some slightly positive value of $0.5 - \bar{q}$. We also need to avoid extending our fit to very high $0.5 - \bar{q}$ because our polynomial fitting function might become inaccurate.

The results in the top panel of Fig. A1 show real data and a cubic fit over the range $0.05 - 0.2$. The data are shown in black, while the cubic fit is shown in red — a cubic is the lowest degree polynomial which provides an accurate match to the data. The bottom panel of Fig. A1 shows the intercept as a function of the upper limit to $0.5 - \bar{q}$, with each coloured line used to show results for a different lower limit. All analysis variants show a positive intercept, whose value converges at close to 0.04 for a wide range of lower and upper limits to $0.5 - \bar{q}$. We therefore adopt $P_{eqm} = 0.04$ for the CB model in our main analysis. This assumption is not relevant to modelling uncontaminated WBs, whose actual mass ratio distribution is used
<figure><img src="image_25.png" /><figcaption>**Figure A1.** Top: The cumulative distribution of $0.5 - \tilde{q}$ for WBs with $\tilde{v} < 1$ (black line). We only show results over the range $0.05 - 0.2$ because a $\delta$-function in the probability distribution at zero would be smeared somewhat by observational uncertainties. The red line shows a cubic fit over this range, but also extrapolates it down to zero. Notice that the intercept is slightly positive, indicating the presence of a population of WBs with an exactly equal mass (consistent with El-Badry et al. 2019). Bottom: The above-mentioned intercept is plotted as a function of the upper limit to the fitting range. Different lines show results with different lower limits, as indicated in the legend. Notice that for a wide range of fitting ranges, the results converge around 0.04, which we assume is the likelihood $P_{\text{eqm}}$ of a binary having exactly equal mass components.</figcaption></figure>

when needed (Section 3.2.4). An equal mass fraction of a few percent is in line with the double main sequence apparent in the colour-magnitude diagram of our WB sample (Fig. 2).

## APPENDIX B: WHEN ORBITAL MOTION BECOMES NON-LINEAR

We assume uniform rectilinear motion for the WB and any undetected CB companions, implicitly assuming that the orbital periods are very long. This should be a very good approximation for the WB given the kAU separations of even the tightest WBs we consider. However, CB orbital periods can be much shorter. This can substantially reduce the mean motion of the CB over the Gaia DR3 observing baseline of $t_G = 34$ months. In this section, we estimate the CB orbital period below which we can no longer safely assume a very long orbital period P. This is used to estimate a minimum CB semi-major axis $a_{\text{int}}$ when setting up its distribution in Section 3.2.3.

For simplicity, we assume that the CB is on a circular orbit and that the line connecting its components rotates by some angle $\theta$ over

the Gaia observing baseline. We consider the impact of increasing $\theta$ by reducing the CB separation while not altering the other CB orbital parameters. The induced recoil velocity on the star detected as part of the WB depends on the CB orbital velocity $v \propto 1/\sqrt{a_{\text{int}}}$. Bearing in mind Kepler's Third Law that $P \propto a_{\text{int}}^{3/2}$, we get that $v \propto P^{-1/3}$. Since the CB rotation angle $\theta$ over a fixed duration of time varies as $\theta \propto 1/P$, we get that

$$v \propto \theta^{1/3}. \qquad (B1)$$

The positive exponent captures the fact that larger $\theta$ is associated with a tighter CB on a faster orbit.

However, the mean velocity $\bar{v}$ will in general involve some additional shape factor $s(\theta)$ that accounts for the extent to which the orbital arc over the Gaia observing baseline deviates from a straight line. Basic trigonometry tells us that the linear distance between two points on a circle is smaller than the distance between them along the circumference by a factor of $\sin(\theta/2)$. Putting this extra factor into equation (B1) tells us that

$$\bar{v} \propto \theta^{-2/3} \left| \sin \left( \frac{\theta}{2} \right) \right|. \qquad (B2)$$

Treating this as a function of $u \equiv \theta/2$, we find that the maximum occurs when $\tan u = 3u/2$. The first non-trivial solution is $u = 0.97$ or $\theta = 111^{\circ}$. Rotation by this angle over $t_G$ implies $P = 3.25 t_G$, so we expect our linear motion approximation to break down when $P < 110$ months or 9.2 yr.

The corresponding orbital separation depends on the mass of the CB. For a low-mass undetected companion, the CB total mass would only slightly exceed the mass of the contaminated star, which is just one of the two stars detected as a WB. The blue bars in Fig. 5 show that the stars in our WB sample easily reach down to $0.4 M_{\odot}$. At this mass, an orbital period of 9.2 years corresponds to $a_{\text{int}} = 3.2$ AU. While CBs with even smaller separations would still affect the inferred WB orbital velocity, the impact rapidly becomes smaller than implied by equation (B1).

We note that 3.2 AU is much larger than the minimum $a_{\text{int}}$ of 0.1 AU considered in the orbit integrations of PS23 (see their section 3.2). CBs with a separation of 0.1 AU would complete many orbits over $t_G$, so the effect of such tight CBs may have been overestimated in their analysis.

## APPENDIX C: COMPARING THE BEST NEWTONIAN AND MOND MODELS

In Fig. 13, we presented the difference in log-likelihood between our best-fitting Newtonian and Milgromian models. The panels in this figure showed the results for each $r_{\text{sky}}$ but with all $\tilde{v}$ pixels summed over, and vice versa. To give a better understanding of which pixels work better in which theory, we show the full 2D distribution of $\Delta\ln P$ between these models in Fig. C1. Since $\Delta\ln P$ is close to zero in several pixels, we clarify which model does better by adding an open white circle at the centre of a pixel if MOND fits it better. Out of the 540 pixels used in our analysis (Table 1), only 221 pixels (41 per cent) prefer MOND while the remaining 59 per cent prefer Newtonian gravity. However, we also need to consider that the pixels which prefer MOND only do so to a rather small extent, while sometimes pixels which prefer Newtonian gravity do so to a very large extent.
<figure><img src="image_26.png" /><figcaption>**Figure C1.** The relative performance of our best-fitting Newtonian and MOND models, shown here for every pixel. The binomial log-likelihood of each model is found by comparison with the observed number of WBs in each pixel (Fig. 6). Open white circles indicate pixels for which MOND outperforms Newtonian gravity. The combined difference in log-likelihood at each $r_{sky}$ summed across all $\tilde{v}$ (and vice versa) is shown in Fig. 13.</figcaption></figure>

## APPENDIX D: THE BEST FIT MODEL WITH REDUCED CB CONTAMINATION

The only alteration to our nominal analysis which gives an appreciable preference for $\alpha_{grav} > 0$ is the one in which we fix $f_{CB} = 0.3$ (Table 3). We explore this model in more detail to understand if it plausibly fits the WB data with a gravity law closer to Milgromian.

In Fig. D1, we compare the observed $\tilde{v}$ distribution to the prediction of the best model with our nominal assumptions and with $f_{CB} = 0.3$, considering both the final result of the gradient ascent and the whole MCMC chain in each case. It is clear that fixing $f_{CB}$ to such a low value causes a catastrophic disagreement with nearly all aspects of the observations. The lack of sufficient CB contamination causes the peak of the $\tilde{v}$ distribution to be more pronounced than in the Gaia number counts at very low $r_{sky}$, though we see that both models struggle somewhat in the peak region at intermediate $r_{sky}$. The most important differences concern the tail of the distribution, which cannot be adequately fit with so few CBs. The analysis has to try and fit the extended tail with substantially more LOS contamination, quadrupling the inferred $f_{LOS}$. This may help somewhat at low $r_{sky}$, but it leads to a rapidly rising $\tilde{v}$ distribution at high $r_{sky}$, in strong disagreement with the observations (bottom right panel).

These very serious issues cause the overall fit to be poorer than in our best-fitting nominal analysis by $\Delta \ln P = 413$, which implies that the latter is preferred at a confidence equivalent to $\sqrt{2\Delta \ln P} = 29\sigma$ for a 1D Gaussian. Clearly, assuming a substantially lower likelihood of CB contamination is not a viable proposition. In any case, this only shifts the inferred gravity law towards MOND by a small amount, with Newtonian gravity still strongly preferred despite facing just over $3\sigma$ tension (Fig. 18). It is therefore extremely difficult to reconcile our WB results with MOND as usually understood given other constraints.

## APPENDIX E: THE ARTIFICIAL VELOCITY LIMIT IN THE CHAE (2023) ANALYSIS

In Section 5.2, we discussed several problems with the WB analysis of Chae (2023). One additional problem which may be noteworthy is that the sample is based on the WB catalogue of El-Badry, Rix & Heintz (2021). Each WB is assigned a likelihood of being genuine based on several factors, one of the most important being the magnitude of $v_{sky}$ (see their equation 7).<fnref n="20" /> Unlike our limit that $\tilde{v} < 5$, their limit to $v_{sky}$ does not easily translate into a limit on $\tilde{v}$ because there is no allowance for more massive WBs having a faster Newtonian $v_c$. Those authors require that

$$v_{sky} < 2.1 \sqrt{\frac{\text{kAU}}{r_{sky}}} \text{ km s}^{-1}, \qquad (\text{E1})$$

which in terms of $\tilde{v}$ corresponds to

$$\tilde{v} < \sqrt{\frac{5 M_{\odot}}{M}}. \qquad (\text{E2})$$

We illustrate this cutoff in Fig. E1, where we plot $\tilde{v}$ as a function of $M$ for the 26615 WBs in the nominal sample of Chae (2023). Nearly all of these WBs lie below the blue line showing the limit given by equation (E2). This is not a completely strict limit because each WB is required to be consistent with this bound within the $2\sigma$ observational uncertainty on $v_{sky}$ (see section 2 of El-Badry, Rix & Heintz 2021). Even so, it should be clear that a WB sample selected in this way is not ideal for the WBT.<fnref n="21" />

To check if our results in Section 5.2 remain reliable despite this complication, we need to impose an upper limit on the mass so that results remain reliable up to some maximal $\tilde{v}$. This is equivalent to drawing a rectangular cut towards the bottom left of Fig. E1. To retain most of the WB sample while having reliable results beyond the main peak of the $\tilde{v}$ distribution, we impose a limit of $M < 1.25 M_{\odot}$. With this mass limit, regardless of the mass of a WB, there is no truncation to the $\tilde{v}$ distribution for $\tilde{v} < 2$. We then use Fig. E2 to show the median $\tilde{v}$ as a function of $r_{sky}/r_M$ for the case where the sample is further limited to only those WBs where the estimated $\sigma(\tilde{v}) < 0.1 \max(1, \tilde{v}/2)$, as done in Fig. 21. Both figures show very similar results, indicating that the $\tilde{v}$ limit imposed by equation (E2) has little impact on the flat trend of the median $\tilde{v}$ with respect to our proxy for the internal WB acceleration. This is presumably because the main peak of the $\tilde{v}$ distribution lies at $\tilde{v} < 1$, which is not much affected by this cut (Fig. E1). Even so, having a mass-dependent limit to $\tilde{v}$ could in principle bias the WBT because a WB with a lower mass also has a lower acceleration at the same $r_{sky}$. Since the main idea of the WBT is that some typical measure of $\tilde{v} \propto \sqrt{g/g_N}$, one should be cautious about a sample selection that allows lower mass systems to reach a higher $\tilde{v}$ than higher mass systems.

<footnote n="20">²⁰This is estimated using the difference in proper motions without allowing for the systemic RV, which should be a reasonable approximation given the small angular sizes of WBs.</footnote>

<footnote n="21">²¹K. El-Badry, private communication.</footnote>
<figure><img src="image_27.png" /></figure>

**Figure D1.** Comparison between the observed $\tilde{v}$ distribution in four different $r_{\text{sky}}$ ranges (solid red bars in different panels) with the prediction of our best model from our nominal analysis (black lines with error bars) and from our revised analysis with fixed $f_{\text{CB}} = 0.3$ (blue lines with error bars). This revised analysis is the only one to appreciably shift the inferred gravity law towards MOND (Fig. 18). The fit is much poorer in this case. The results shown here are normalized based on the observed WB distribution, helping to highlight that the $f_{\text{CB}} = 0.3$ model predicts too many WBs at high $r_{\text{sky}}$.

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

**Figure E1.** The red dots show the distribution of WB mass and $\tilde{v}$ for the nominal sample of Chae (2023). Notice how the vast majority of WBs lie below the solid blue line [equation (E2)] due to selection effects (equation 7 of El-Badry, Rix & Heintz 2021).

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

**Figure E2.** Similar to Fig. 21, but restricted further to only WBs with $M < 1.25 M_{\odot}$. This ensures that the artificial limit to the $\tilde{v}$ distribution in the Chae (2023) sample (blue line in Fig. E1) does not influence the results for $\tilde{v} < 2$. The mass limit has little impact on the results, which continue to show a flat trend once we impose that $\sigma(\tilde{v}) < 0.1 \max(1, \tilde{v}/2)$, as done here.

This paper has been typeset from a TeX/LATeX file prepared by the author.