pdftexdestination with the same hyperrefOption hyperrefToken pdftex(dest)

# Strong constraints on the gravitational law from _Gaia_ DR3 wide binaries

Indranil Banik1, Charalambos Pittordis2, Will Sutherland2, Benoit Famaey3,
Rodrigo Ibata3, Steffen Mieske4 and Hongsheng Zhao1
1Scottish Universities Physics Alliance, University of Saint Andrews, North
Haugh, Saint Andrews, Fife, KY16 9SS, UK
2The School of Physical and Chemical Sciences, Queen Mary University of
London, Mile End Road, London, E1 4NS, UK
3Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de
Strasbourg, 11 rue de l’Université, 67000 Strasbourg, France
4European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago,
Chile Email<EMAIL_ADDRESS>(Indranil Banik);
<EMAIL_ADDRESS>(Will Sutherland)

###### Abstract

We test Milgromian dynamics (MOND) using wide binary stars (WBs) with
separations of $2-30$ kAU. Locally, the WB orbital velocity in MOND should
exceed the Newtonian prediction by $\approx 20\%$ at asymptotically large
separations given the Galactic external field effect (EFE). We investigate
this with a detailed statistical analysis of _Gaia_ DR3 data on 8611 WBs
within 250 pc of the Sun. Orbits are integrated in a rigorously calculated
gravitational field that directly includes the EFE. We also allow line of
sight contamination and undetected close binary companions to the stars in
each WB. We interpolate between the Newtonian and Milgromian predictions using
the parameter $\alpha_{\rm{grav}}$, with 0 indicating Newtonian gravity and 1
indicating MOND. Directly comparing the best Newtonian and Milgromian models
reveals that Newtonian dynamics is preferred at $19\sigma$ confidence. Using a
complementary Markov Chain Monte Carlo analysis, we find that
$\alpha_{\rm{grav}}=-0.021^{+0.065}_{-0.045}$, which is fully consistent with
Newtonian gravity but excludes MOND at $16\sigma$ confidence. This is in line
with the similar result of Pittordis and Sutherland using a somewhat different
sample selection and less thoroughly explored population model. We show that
although our best-fitting model does not fully reproduce the observations, an
overwhelmingly strong preference for Newtonian gravity remains in a
considerable range of variations to our analysis. Adapting the MOND
interpolating function to explain this result would cause tension with
rotation curve constraints. We discuss the broader implications of our results
in light of other works, concluding that MOND must be substantially modified
on small scales to account for local WBs.

###### keywords:

gravitation – binaries: general – celestial mechanics – stars: kinematics and
dynamics – galaxies: kinematics and dynamics – methods: statistical

††pubyear: 2023††pagerange: Strong constraints on the gravitational law from
_Gaia_ DR3 wide binaries–E

##  1 Introduction

Our current best understanding of gravity is encapsulated by the theory of
General Relativity (GR; Einstein, 1915). This reduces to Newtonian dynamics in
the weak-field non-relativistic limit (Rowland, 2015; de Almeida, Piattella &
Rodrigues, 2016; Ciotti, 2022). According to this theory, the gravity from a
point mass follows an inverse square law, which historically was inferred from
the fact that the rotation velocity $v_{c}$ of a planet around the Sun
declines with its heliocentric distance $r$ as $v_{c}\propto 1/\sqrt{r}$. In
principle, this Keplerian decline should also apply to galaxies as a whole if
one considers tracers of their potential far beyond the bulk of their
detectable baryonic mass (stars and gas). However, it is well known that
galaxy rotation curves are flat (Bosma, 1978; Rubin, Ford & Thonnard, 1978;
Faber & Gallagher, 1979, and references therein) $-$ this is nicely
illustrated in figure 15 of Famaey & McGaugh (2012).

Several explanations have been put forward for this missing gravity problem,
which is also apparent in a number of other ways like dwarf galaxy velocity
dispersions (figure 11 of McConnachie, 2012) and weak gravitational lensing
(Brouwer et al., 2021). The most popular idea is that the extra gravity needed
to bind galaxies and enhance their effect on passing photons comes from an
otherwise undetected halo of particles beyond the well-tested standard model
of particle physics (for a review, see Peebles, 2017). This forms the basis
for the currently prevailing Lambda cold dark matter ($\Lambda$CDM) standard
model of cosmology (Efstathiou, Sutherland & Maddox, 1990; Ostriker &
Steinhardt, 1995). However, DM particles have not been found despite decades
of highly sensitive searches for collisions with terrestrial detectors (LUX
Collaboration, 2017; Lux-Zeplin Collaboration, 2023). Similarly, there is no
sign of $\gamma$-rays from DM annihilation or decay processes in supposedly
DM-dominated Galactic satellites, providing stringent limits on the allowed
properties of the hypothetical particles (Hoof, Geringer-Sameth & Trotta,
2020).

This has led some workers to consider that the missing gravity problem might
not be due to seemingly undetectable particles. The most developed alternative
is Milgromian dynamics (MOND; Milgrom, 1983). It proposes that the gravity $g$
in an isolated spherically symmetric system is asymptotically related to the
Newtonian gravity $g_{{}_{N}}$ of the baryons alone according to

$\displaystyle
g~{}\to~{}\begin{cases}g_{{}_{N}}\,,&\textrm{if}~{}g_{{}_{N}}\gg
a_{{}_{0}}\,,\\\ \sqrt{a_{{}_{0}}g_{{}_{N}}}\,,&\textrm{if}~{}g_{{}_{N}}\ll
a_{{}_{0}}\,.\end{cases}$ (1)

MOND introduces a fundamental acceleration scale $a_{{}_{0}}$ below which the
gravitational field deviates substantially from Newtonian expectations.111This
is true in percentage terms, but the difference in absolute terms is $\la
a_{{}_{0}}$. Just like with the Newtonian constant $G$, we must determine
$a_{{}_{0}}$ empirically. Using the rotation curves of a handful of disc
galaxies with properties similar to our own, it has been found that
$a_{{}_{0}}=1.2\times 10^{-10}$ m/s2 (Begeman, Broeils & Sanders, 1991). This
value has remained very stable over many decades (Gentile, Famaey & de Blok,
2011; McGaugh, Lelli & Schombert, 2016; Chae, 2022).

MOND works remarkably well at predicting the dynamics of galaxies across a
huge range in baryonic mass, surface brightness, and gas fraction (Famaey &
McGaugh, 2012; McGaugh et al., 2016; Lelli et al., 2017b; Kroupa et al., 2018;
Li et al., 2018). Observations show that galaxies follow a remarkably tight
‘radial acceleration relation’ (RAR) between the radial components of $\bm{g}$
and $\bm{g}_{{}_{N}}$ as deduced from rotation curves and the baryonic mass
distribution, respectively (McGaugh, 2020). This conclusion is not much
dependent on assumptions about the stellar mass-to-light ratio because one can
restrict attention to gas-rich galaxies (McGaugh, 2012). The intrinsic scatter
in the RAR is only 0.034 dex, which given the uncertainties is consistent with
zero (Desmond, 2023). The RAR also extends to elliptical galaxies and even to
weak gravitational lensing down to $g_{{}_{N}}\approx 10^{-5}\,a_{{}_{0}}$
(Brouwer et al., 2021). The gravitational forces in virialized galaxies are
thus very well described by MOND (for a review, see section 3 of Banik & Zhao,
2022).

These successes are complemented by recent numerical MOND simulations of
galaxies, which allow an exploration of non-axisymmetric dynamical features
like bars and spiral arms (Banik et al., 2020; Roshan et al., 2021a). MOND has
also been shown to get about the right star formation rates in disc galaxies
with an appropriate size and gas fraction for their stellar mass (Nagesh et
al., 2023). While these works typically target isolated galaxies, MOND has
also been applied to interacting galaxies, in particular to understand the
Antennae (Renaud, Famaey & Kroupa, 2016), the tidal stability of dwarf
spheroidal galaxies in the gravitational field of the Fornax Cluster (Asencio
et al., 2022), the asymmetric tidal tails of the globular cluster Palomar 5
(Thomas et al., 2018) and open star clusters in the Solar neighbourhood
(Kroupa et al., 2022), and the formation of the Local Group (LG) satellite
planes around the Milky Way (MW) and M31 from tidal debris expelled by a past
flyby encounter between their discs around 8 Gyr ago (Zhao et al., 2013; Bílek
et al., 2018; Banik et al., 2022).

Many of these results were motivated by tensions and inconsistencies within
the $\Lambda$CDM paradigm. For instance, it is known to face severe
difficulties explaining the LG satellite galaxy planes (Pawlowski, 2021b, and
references therein). Hydrodynamical MOND simulations of M33 (Banik et al.,
2020) were motivated by difficulties in understanding why it is weakly barred,
bulgeless, and has a two-armed spiral $-$ behaviour which is in many ways the
opposite of what arises in Newtonian simulations with a live halo (Sellwood,
Shen & Li, 2019). Problems understanding galaxy bars in $\Lambda$CDM are
certainly not confined to M33: the model faces a $13\sigma$ tension when
confronted with the observed ratios between bar lengths and corotation radii
(Roshan et al., 2021b). This fast bar tension has been rather directly linked
to dynamical friction between the rotating bar and the hypothetical CDM halo
(figure 12 of Roshan et al., 2021a). Those authors demonstrated that the
problem is alleviated in MOND (see their figure 20), with a similar result
found in subsequent more detailed simulations that include hydrodynamics and
star formation (see figure 16 of Nagesh et al., 2023). The presence of bars in
galaxies with a low surface brightness is also very problematic for
$\Lambda$CDM because these galaxies are supposedly embedded in a dominant dark
halo, which reduces the role of disc self-gravity (Kashfi, Roshan & Famaey,
2023). Milgromian disc galaxies are always completely self-gravitating, making
it easier to understand why even galaxies with a large enhancement to their
Newtonian baryonic rotation curve commonly have bars and spiral arms (McGaugh,
Schombert & Bothun, 1995; McGaugh & de Blok, 1998a, b).

MOND has also been applied on galaxy cluster and cosmological scales with the
neutrino hot dark matter ($\nu$HDM) paradigm developed by Angus (2009) and the
aether scalar tensor (AEST) model (Skordis & Złośnik, 2019, 2021, 2022) $-$
for a review, we refer the reader to section 9 of Banik & Zhao (2022). Here
again some of the successes relate to situations which are severely
problematic for $\Lambda$CDM, in particular the KBC void and Hubble tension
(Keenan, Barger & Cowie, 2013; Haslbauer, Banik & Kroupa, 2020) and the El
Gordo massive galaxy cluster collision at redshift $z=0.87$ (Asencio, Banik &
Kroupa, 2021, 2023). Thus, the successes of MOND extend well beyond its
initial motivation forty years ago from the rotation curves of high surface
brightness disc galaxies. At the same time, the motivation for considering
alternatives to $\Lambda$CDM extends well beyond a mere failure to detect the
hypothetical particles (for an extensive review, see Banik & Zhao, 2022).

###  1.1 The importance of wide binaries (WBs)

In this contribution, we seek to test the validity of MOND on very small
length scales by galactic standards, greatly diminishing the role of DM
(Acedo, 2020). One might think that the role of MOND would also be diminished,
but since MOND does not introduce a fundamental new length scale, it should
remain valid in systems with a much smaller mass and size than a typical
galaxy provided the accelerations are low enough. Equation 1 indicates that
Newtonian and Milgromian gravity diverge when $g\la a_{{}_{0}}$. Around an
isolated point mass $M$, this occurs beyond its MOND radius

$\displaystyle r_{{}_{M}}~{}\equiv~{}\sqrt{\frac{GM}{a_{{}_{0}}}}\,.$ (2)

The MOND radius of the Sun is 7000 AU (7 kAU) or 0.03 pc. This is much smaller
than the typical separation between stars in the Solar neighbourhood, where
the Newtonian Jacobi or tidal radius for a system with $M=2\,M_{\odot}$ is 350
kAU (equation 43 of Jiang & Tremaine, 2010). It is thus possible to test MOND
at kAU distances from a star provided suitable tracers can be found. The
properties of long-period comets that reach kAU distances have been used to
argue in favour of the Milgromian enhancement to the Solar potential (Penner,
2020), while its predicted anisotropy may explain the apparent orbital
alignment of Kuiper belt objects with a large pericentre (Brown & Mathur,
2023; Migaszewski, 2023). However, this evidence is rather indirect. More
direct evidence could be obtained using interstellar precursor missions
travelling at a few percent of the speed of light $c$ (Banik & Kroupa, 2019b),
but this remains well beyond current technology.

These difficulties can be overcome by using stars in WBs, defined in this work
as binary stars with orbital acceleration $\la a_{{}_{0}}$, implying kAU-scale
separations. The nearest star to the Sun is in such a WB (Beech, 2009, 2011).
MOND can thus be tested using the orbital acceleration of Proxima Centauri
around $\alpha$ Centauri A and B, which have a very small separation and can
be treated as one object (Kervella, Thévenin & Lovis, 2017). Unfortunately,
the $\approx 0.5\,\mu$as astrometric precision required is roughly two orders
of magnitude beyond the reach of current instruments (Banik & Kroupa, 2019a).
Instead of using the acceleration, it is much more promising to use the
relative velocity, but then the wide binary test (WBT) of gravity has to be
done statistically by considering a large number of WBs to average over
orbital phases and projection effects. This was first suggested by Hernandez,
Jiménez & Allen (2012), with Pittordis & Sutherland (2018) later exploring the
behaviour of WBs in several gravity theories. Since WBs come in a range of
masses and separations, the WBT is usually phrased as a statistical test
involving the distribution of the dimensionless parameter

$\displaystyle\widetilde{v}~{}\equiv~{}v_{\rm{rel}}\div\overbrace{\sqrt{\frac{GM}{r_{\rm{sky}}}}}^{\rm{Newtonian}\,v_{c}},$
(3)

where $M$ is the total mass of a WB with sky-projected separation
$r_{\rm{sky}}$ and relative velocity $\bm{v}_{\rm{rel}}$, with
$v\equiv\left|\bm{v}\right|$ for any vector $\bm{v}$. Thus, the relative
velocity is normalized to the Newtonian circular velocity at the projected
separation. Since the actual separation is larger, $\widetilde{v}$ is smaller
than if the Newtonian $v_{c}$ had been calculated using the 3D separation.
Moreover, we only consider the sky-projected part of $\bm{v}_{\rm{rel}}$
(which we denote $\bm{v}_{\rm{sky}}$) as the data is not yet accurate enough
for a 3D version of the WBT (though see Section 5.5). This is mainly due to
the parallaxes being too inaccurate to reliably constrain the 3D separation
(Section 2.3.2), though there is also some additional uncertainty from the
required correction for the gravitational redshift of each star (Loeb, 2022).
Considering only sky-projected quantities further reduces $\widetilde{v}$,
which in this contribution is based on using $\bm{v}_{\rm{sky}}$ in Equation 3
$-$ we do not use the full 3D $\widetilde{v}$. Note that even though we will
mostly be dealing with the magnitudes of $\bm{r}_{\rm{sky}}$ and
$\bm{v}_{\rm{sky}}$ in this work, these are actually 2D vectors within the sky
plane. This is important to the work of Hwang, Ting & Zakamska (2022), who
used the angle $\psi$ between $\bm{r}_{\rm{sky}}$ and $\bm{v}_{\rm{sky}}$ to
obtain interesting constraints on the WB orbital eccentricity distribution
(Section 3.1.2).

The first detailed MOND calculations of WB orbits showed that the Milgromian
circular velocity would be enhanced by 20% over the Newtonian expectation in
the asymptotic regime of large separations, a result which can be understood
analytically (Banik & Zhao, 2018c, hereafter BZ18).222$\widetilde{v}$ in this
contribution has the same meaning as $\widetilde{v}_{\rm{sky}}$ in BZ18.
Numerical MOND simulations of nearby star clusters also show a 20% enhanced
velocity dispersion over their Newtonian counterparts for essentially the same
reason (Kroupa et al., 2022). The enhancement is limited by the Galactic
gravity on the Solar neighbourhood, which is important in MOND because it is a
non-linear theory (Equation 1). As a result, the internal behaviour of a self-
gravitating system is influenced by an external source of gravity, even in the
absence of tidal effects. This non-standard external field effect (EFE) has
been an important aspect of MOND since the beginning, not only because it is
theoretically inevitable in a non-linear theory, but also because the EFE is
necessary to suppress the MOND enhancement to the velocity dispersions of the
Pleiades and Praesepe open star clusters (see section 3 of Milgrom, 1983).

Many further arguments for the EFE can be made nowadays (for a review, see
section 3.3 of Banik & Zhao, 2022). It significantly enhances the tidal
susceptibility of dwarf galaxies in a cluster environment by weakening their
self-gravity, thereby helping to explain the observed signs of tidal
disturbance in the Fornax Cluster’s dwarf galaxies and the lack of low surface
brightness dwarfs towards its centre (Asencio et al., 2022). As explained in
that work, tidal radii are necessarily in the EFE-dominated regime, making
their results particularly sensitive to the EFE. Its subtle imprint on the
outer parts of disc galaxy rotation curves has recently been detected by
comparing galaxies in isolated and more crowded environments (Chae et al.,
2020b, 2021). The EFE is also crucial when calculating the escape velocity
from a mass distribution, since the Milgromian potential of an isolated mass
is infinitely deep (Famaey, Bruneton & Zhao, 2007; Zhao, Li & Bienaymé, 2010;
Banik & Zhao, 2018b). This is not true once the EFE is included because at
large distances where the EFE dominates, the gravitational field of a point
mass returns to the Newtonian inverse square law, albeit with a higher
normalization and some angular dependence (see figure 1 of Banik & Zhao,
2022). This EFE-dominated regime is sometimes known as the quasi-Newtonian
regime, though in this contribution we simply refer to it as the asymptotic
regime because local WBs lack an extended region in which they are isolated
and at low acceleration. This is due to the Galactic external field $g_{e}$
being slightly stronger than $a_{{}_{0}}$, as evident from the Galactic
rotation curve (e.g. Zhu et al., 2023).

Given the importance of $g_{e}$ for the WBT, it is fortunate that we know
$\bm{g}_{e}$ from a direct measurement of the acceleration of the Solar System
relative to distant quasars using their changing aberration angle, which
directly tells us that $\bm{g}_{e}/c=5.05\pm 0.35\,\mu$as/yr in a direction
very close to that of the Galactic centre (Gaia Collaboration, 2021b). This
precludes substantial deviations of $g_{e}$ from the value deduced
kinematically from the Galactic rotation curve (Section 3.1). Even so, the
predicted 20% enhancement to the orbital velocities of local WBs is related to
the shape of the MOND interpolating function in the regime close to the
critical acceleration scale $a_{{}_{0}}$. It is thus somewhat model-dependent
as it does not directly relate to the weak-field asymptotic limit at the heart
of MOND. Nevertheless, recent detailed rotation curve studies have mostly
converged on the form of the interpolating function in this regime, though of
course slight variations are still possible. We will return to this point in
Section 5.3.

Prospects for the WBT were discussed extensively in BZ18, whose section 8
clarified the main systematic effects that would likely have to be considered.
The main issue was expected to be undetected close binary (CB) companions to
one or both of the stars in a WB. Indeed, the nearest star to the Sun is in
just such a triple system: Proxima Centauri is on a wide orbit about $\alpha$
Centauri A and B, whose mutual orbital semi-major axis of 23 AU is negligible
compared to the 13 kAU distance from their barycentre to Proxima Centauri
(Kervella et al., 2016, 2017). At larger distances from the Sun, it is
possible that a similar CB would not be resolved, but it would still have
significant effects on the observed kinematics that would not arise if Proxima
Centauri were orbiting a single star (Section 3.2.1). Pittordis & Sutherland
(2019, hereafter PS19) considered the WBT in light of actual data from _Gaia_
Data Release 2 (DR2; Gaia Collaboration, 2018), highlighting that there is
indeed an extended tail towards much higher $\widetilde{v}$ than could
plausibly arise from genuine WB orbital motion, which cannot realistically
yield $\widetilde{v}\ga 2$ for any reasonable modification to gravity. The
possible presence of CBs was strongly suggested by the analysis of Clarke
(2020) and later by Belokurov et al. (2020), who argued that some of the WBs
with unexpectedly large relative velocities also have a _Gaia_ astrometric
solution that poorly fits the observations. Since the astrometric solution
only includes parallax and proper motion, it was argued that the excess
‘noise’ could be due to astrometric acceleration induced by a CB, which was
previously shown to typically induce a much larger orbital acceleration than
the WB orbital motion (see section 8.2 of BZ18, ). Further follow-up
observations should be able to confirm the existence of CB companions around
some WBs (Manchanda, Sutherland & Pittordis, 2023), especially once the
astrometric time series are published in _Gaia_ DR4 (see Section 5.5).

Some line of sight (LOS) contamination is also inevitable, though this is
expected to be rather modest thanks to the excellent quality of the _Gaia_
data (see section 3 of PS19, ). Their work shows a clear main peak to the
$\widetilde{v}$ distribution due to WB orbital motion. Investigation of its
properties already decisively rejects theories where local WBs behave like
isolated MOND systems without the weakening of their self-gravity due to the
Galactic EFE (see their figure 11). This is because the gravity binding WBs
can be deduced from the main peak to the $\widetilde{v}$ distribution at
$\widetilde{v}\la 2$, with CB and LOS contamination leading to an extended low
amplitude tail out to much higher $\widetilde{v}$. If the orbital velocity
becomes asymptotically flat with increasing separation, we would expect the
mode of the $\widetilde{v}$ distribution to be significantly higher in systems
with large $r_{\rm{sky}}$. The predicted peak shift is in catastrophic
disagreement with the observations, which show at most only a modest peak
shift. This still leaves open the possibility that WBs obey MOND with its
inevitable EFE, but it falsifies the quantized inertia proposal (McCulloch &
Lucio, 2019) and repeated claims to have confirmed that WBs follow the
isolated MOND prediction without considering the full $\widetilde{v}$
distribution (Hernandez et al., 2019, 2022).

In this contribution, we conduct a detailed statistical analysis to test the
MOND prediction for local WBs. We use a rigorous grid solution to the MOND
field equation including the EFE (BZ18) and then add possible perturbations
from undetected CBs around one or both of the stars in each WB. We jointly
infer properties of the WB population, the undetected CB population, and the
extent of LOS contamination, which becomes important towards large
$r_{\rm{sky}}$ and $\widetilde{v}$. Both the model used and the exploration of
its parameter space are much more detailed than in any previous attempt at the
WBT. To mitigate possible biases, the plan was prepared in advance (Banik,
Pittordis & Sutherland, 2021) and has barely been modified to deal with the
actual _Gaia_ DR3 (Gaia Collaboration, 2023).

In the following, we explain how we reduce the _Gaia_ dataset to a form
suitable for our analysis and what quality cuts we employ (Section 2). We then
introduce our detailed model for the WB dataset and explain how we compare it
to observations (Section 3). Our results are presented in Section 4 and
discussed in Section 5, which includes a comparison with prior WBT results
(Section 5.2). We conclude in Section 6.

##  2 The observed WB distribution

The primary dataset for the WBT comes from the precise results obtained by the
_Gaia_ mission (Gaia Collaboration, 2016). Local WBs were extracted from the
_Gaia_ dataset using the methods described in Pittordis & Sutherland (2023,
hereafter PS23), who slightly adapted the methods in an earlier work (section
2 of PS19, ) based on _Gaia_ DR2 (Gaia Collaboration, 2018) due to the
improved quality in _Gaia_ Early Data Release 3 (EDR3; Gaia Collaboration,
2021a). We supplement EDR3 with radial velocities (RVs) from the full _Gaia_
DR3 (Gaia Collaboration, 2023).

In Section 2.1, we describe the basic quality cuts that we impose. We then
conduct a Monte Carlo propagation of the uncertainties (Section 2.3). This
leads us to impose further quality cuts, which we briefly summarize in Section
2.4 alongside our choice of restrictions on the parameter ranges.

###  2.1 The _Gaia_ DR3 sample & basic quality cuts

The details of the sample selection are given in PS23, whose main points we
summarize below. We begin by selecting all stars from _Gaia_ EDR3 at a
Galactic latitude $\lvert b\rvert>15^{\circ}$ with an apparent _Gaia_ -band
($G$-band) magnitude $m_{G}<17$ and measured parallax $\varpi>4$ mas
(estimated distance $<250$ pc) uncorrected for parallax bias (Lindegren et
al., 2021). From this sample of 2.1 million single stars, candidate WBs are
selected by requiring:

  1. 1.

Projected separation $r_{\rm{sky}}\leq 50$ kAU;

  2. 2.

Star distances consistent with each other within the lesser of 8 pc or
$4\times$ the combined distance uncertainty, i.e. $\lvert
d_{A}-d_{B}\rvert\leq\min\left(4\sigma_{d},8\,\rm{pc}\right)$, where $A$ and
$B$ label the stars in each candidate WB; and

  3. 3.

Projected velocity difference between the stars of $v_{\rm{sky}}<3$ km/s, as
inferred from the difference in proper motions but assuming both stars are at
the mean estimated distance $\bar{d}\equiv\left(d_{A}+d_{B}\right)/2$.

From this preliminary list of binaries (WB-EDR3), sky regions are removed
around four known open clusters (see table 1 of PS23, ).

To remove some probable triple systems or groups, we reject all WBs in which
either star is common to more than one candidate binary in WB-EDR3. We also
search for comoving companion stars to a fainter limit using a ‘faint star’
sample constructed from _Gaia_ EDR3 stars with $m_{G}<20$ and parallax
$\varpi>10/3$ mas. For each star in each candidate WB, we search for
companions in the faint star sample with the following criteria:

  1. 1.

Parallax distance consistent with the main star at $4\sigma$;

  2. 2.

Angular separation less than 2/3 of the main binary separation (since
hierarchical triples are expected to be unstable if the inner orbit separation
$\ga 0.4\times$ the outer separation);

  3. 3.

Angular separation $>0.5\arcsec$ to avoid barely-resolved companions; and

  4. 4.

Projected velocity difference from the main star $\leq 5$ km/s.

If any such ‘third star’ is found, the candidate binary is rejected. Note
however that this will not reject triples with an inner orbit closer than
about 100 AU or where the third star is so faint that $m_{G}>20$, so it
remains crucial to model the CB population. We return to this point in Section
3.2.3.

PS23 also applied a per-star quality cut based on the _Gaia_ parameters (see
equation 1 of Arenou et al., 2018). WBs were only considered if both stars
satisfy

$\displaystyle\sqrt{\frac{\chi^{2}}{\nu}}$ $\displaystyle\leq$ $\displaystyle
1.2\max\left(1,\,\exp\left[-0.2\left(m_{G}-19.5\right)\right]\right)\,,$ (4)
$\displaystyle\chi^{2}$ $\displaystyle\equiv$ $\displaystyle{\tt
astrometric\\_chi2\\_al}\,,$ $\displaystyle\nu$ $\displaystyle\equiv$
$\displaystyle{\tt astrometric\\_n\\_good\\_obs\\_al}-5\,.$

This cut is based on _Gaia_ DR2, so it should be fairly conservative in _Gaia_
DR3 given the greater number of observations per star and the longer time
baseline. Throughout this study, any quality cuts applied at the individual
star level are implemented by excluding the whole WB if either of its two
stars fails to pass the relevant quality cut.

Applying all the above quality cuts yields a sample of 73k candidate WBs. The
majority have $r_{\rm{sky}}\leq 1$ kAU, which is not very useful for the WBT.
We therefore publicly release a smaller version of this sample consisting of
19786 WBs that satisfy $1.5<r_{\rm{sky}}/\rm{kAU}<40$, which is somewhat wider
than the parameter range we use for the WBT.

From the _Gaia_ DR3 WBs that pass the above quality cuts, we restrict
ourselves to the range $r_{\rm{sky}}=2-30$ kAU and $\widetilde{v}\leq 5$.
These choices were fixed in advance to limit the possibility of biasing the
results (Banik et al., 2021). The upper limit to $\widetilde{v}$ is much
larger than can plausibly arise from WB orbital motion. This is necessary to
properly sample the extended tail to the $\widetilde{v}$ distribution, which
we need in order to properly constrain the CB and LOS contamination causing
the tail. In particular, LOS contamination is expected to become more
significant at high $\widetilde{v}$ because $\widetilde{\bm{v}}$ is a 2D
quantity, so we can get much better leverage on the extent of LOS
contamination if we have data at high $\widetilde{v}$. Even so, we do not want
to include data at very high $\widetilde{v}$ because the likelihood of
systematics increases the wider the ‘aperture’ used for the WBT.

We could increase the sample size substantially by going down to even lower
separations, but then there would be a risk of swamping the analysis with data
on Newtonian systems. In this case, small improvements to the fit in the
Newtonian regime would become more important to the statistical analysis than
checking whether the $\widetilde{v}$ distribution differs between systems
within their MOND radius and those with larger separations. On the other hand,
we need to go down to rather low $r_{\rm{sky}}$ to provide a ‘Newtonian
anchor’ population that cannot be affected by MOND, especially since the RAR
requires a fairly smooth MOND interpolating function such that significant
MOND effects persist even when the accelerations slightly exceed $a_{{}_{0}}$
(Section 5.3). Model parameters besides the gravity law can be constrained
much better if we have a significant number of WBs which cannot be much
affected by MOND, since then the parameters inferred here will not be
degenerate with the gravity law. In particular, a Newtonian anchor population
should be extremely valuable in constraining the CB population parameters $-$
WBs are expected to have similar CB populations regardless of whether their WB
orbital acceleration is below or above $a_{{}_{0}}$. Our analysis does not
explicitly consider some subset of our sample to be completely immune to MOND
(unlike Chae, 2023), but instead relies on the predicted MOND effects being
different at low and high $r_{\rm{sky}}$. We will see later that since our
sample goes up to $M=4\,M_{\odot}$ and thus has systems with MOND radii of up
to 14 kAU (Equation 2), we have plenty of WBs separated by a small fraction of
their MOND radius. This is especially true given the steeply declining
$r_{\rm{sky}}$ distribution of WBs (; ; ), which also implies that the 3D
separation is likely to only slightly exceed $r_{\rm{sky}}$ (Equation 8).

###  2.2 The mass-luminosity relation

We estimate the mass of each star in our WB sample from its absolute $G$-band
magnitude $M_{G}$ using a similar technique to PS19, which we briefly describe
below. We use the $M_{V}$ and $\left(V-I\right)$ colour at different stellar
masses $M$ as tabulated in Pecaut & Mamajek
(2013).333https://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt
[3.3.2021] We relate this to _Gaia_ photometry using the relation between
$\left(V-I\right)$ and $\left(G-V\right)$ colours from the first Johnson-
Cousins relation in table C2 of Riello et al. (2021), which states that

$\displaystyle\left(G-V\right)\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!-0.01597-0.02809\left(V-I\right)-0.2483\left(V-I\right)^{2}$
(5)
$\displaystyle\\!\\!\\!+0.03656\left(V-I\right)^{3}-0.002939\left(V-I\right)^{4}.$

Figure 1: The relation between the mass of a star and its absolute $G$-band
magnitude (red line). Results are based on the absolute $V$-band magnitude and
$\left(V-I\right)$ colour tabulated in Pecaut & Mamajek (2013) and the
relation between $\left(V-I\right)$ and $\left(G-V\right)$ colours from the
first Johnson-Cousins relation in table C2 of Riello et al. (2021). Our cubic
fit to the data (blue line; Equation 6) is used over the mass range between
the vertical solid grey lines. Outside this range, a linear relation is
assumed, as indicated by the solid purple lines. The dotted vertical line
shows the maximum mass of any star in our WB sample.

We now have a complete set of equations to obtain $M_{G}$ for any $M$.
However, this is difficult to invert because of small-scale irregularities in
the relation arising from the complexities of stellar astrophysics (Figure 1).
We therefore fit the relation using the cubic

$\displaystyle M_{G}~{}=~{}4.887-5.693x+0.4164x^{2}+0.9611x^{3}\,,$ (6)

where $x\equiv\ln M/M_{\odot}$. Since each unit decrease in $M_{G}$
corresponds to slightly less than an $\mathrm{e}$-fold increase in luminosity
$L$, Equation 6 tells us that $L\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}M^{5}$ when $M=M_{\odot}$ ($x=0$),
representing a steep mass-luminosity relation. We use our cubic fit over the
range $0.6<M_{G}<11.1$. At the low (high) mass end, we assume the relation
becomes linear and has a slope of $-3.3$ ($-2.9$). The fitted relation between
$M_{G}$ and $x$ is shown as the solid blue line, confirming that the above
cubic provides a fairly good fit over the range we use it, i.e. between the
vertical solid grey lines. The purple lines outside this range show that a
linear relation here also provides a reasonably good fit, at least up to the
dotted vertical grey line showing the maximum mass of any star in our sample.
A cubic is the lowest order polynomial which seems to provide a good match to
the tabulated data in Pecaut & Mamajek (2013), mainly because the slope is
similar at low and high masses but there is a steepening at intermediate
masses.

We invert our piece-wise analytic relation between $\ln M/M_{\odot}$ and
$M_{G}$ to obtain the mass of any star from its $M_{G}$, which we obtain from
its apparent magnitude using its trigonometric parallax. The linear relations
assumed at very low and very high masses can be readily inverted. We invert
our cubic fit over most of the mass and $M_{G}$ range using the Newton-Raphson
method. To speed this up, we obtain an initial guess for the mass using a
linear fit of the form

$\displaystyle M_{G}~{}\approx~{}5.023-5.102x\,.$ (7)

This can be inverted analytically without iterative techniques. Our estimate
that $-\frac{d\ln M}{dM_{G}}=0.20$ is similar to previous estimates of
$0.074\ln 10=0.17$ (equation 3 of PS23, ) and $0.0725\ln 10=0.17$ (equation 3
of PS19, ).

We consider uncertainties in $M_{G}$ arising from those in the trigonometric
parallax and the _Gaia_ photometric uncertainty in the apparent magnitude,
which we convert from a fractional uncertainty to a magnitude uncertainty. We
also allow a small uncertainty in the conversion of luminosity to mass. This
must play some role, but the $\approx 42\%$ asymptotic enhancement to the
radial gravity predicted in MOND combined with our estimate that
$\frac{dM_{G}}{d\ln M}=-5.1$ (Equation 7) implies that the MOND signal
corresponds to a difference in $M_{G}$ of about $5.1\ln 1.4=1.7$ magnitudes,
which is very large compared to the slight fluctuations about our cubic mass-
luminosity relation evident in Figure 1. In other words, the steep mass-
luminosity relation implies that mass uncertainties should have a very small
impact compared to uncertainties in the relative velocity arising from
uncertain proper motions, especially when the key quantity entering
$\widetilde{\bm{v}}$ is essentially the relative proper motion between the two
stars in each WB. The mass of a star with known absolute magnitude can be
constrained to within $\approx 6\%$ (Eker et al., 2015; Mann et al., 2015),
but due to the square root term in Equation 3, this would have only a 3%
impact on $\widetilde{v}$. The impact on the WBT is further reduced because it
is mostly concerned with the width of the $\widetilde{v}$ distribution, which
is only affected at second order by measurement uncertainties. Even so, we do
assign a modest 0.024 dex or 5.5% uncertainty to our stellar mass estimates at
fixed luminosity for reasons discussed in Section 2.3.3. We will see later
that our main conclusion is the same as a different study which used a linear
relation between $M_{G}$ and $x$ over the full mass range relevant for the
WBT, thus leading to rather different errors in the conversion of absolute
magnitudes to masses (PS23).

###  2.3 Monte Carlo error propagation

We use the full $5\times 5$ _Gaia_ covariance matrix to propagate
uncertainties in the sky position, parallax, and proper motion of each star in
our WB sample into an uncertainty in $\widetilde{v}$. We will see later that
this uncertainty is very small for a large number of WBs, so our approach is
to consider only those systems and obtain their
$\left(r_{\rm{sky}},\widetilde{v}\right)$ directly from the raw observables.
The resulting $\left(r_{\rm{sky}},\widetilde{v}\right)$ distribution is then
used in our main analysis, which therefore does not directly include
measurement uncertainties $-$ these are included indirectly in the sample
selection. Our procedure avoids adding measurement errors to values which
already contain measurement errors, as that would lead to adding measurement
errors twice to the latent values.

For our Monte Carlo error propagation, we diagonalize the correlation matrix
and generate Gaussian random numbers to propagate uncertainties along each of
the eigenvector directions, bearing in mind also the eigenvalues. The
correlation matrix is dimensionless, so we relate it to the covariance matrix
using the _Gaia_ catalogue uncertainty of each parameter. Uncertainties in the
parallax and proper motion are important, but we do not propagate the
uncertainty in the sky position of each star at the reference epoch, instead
taking the published position at face value. This simplifies our analysis
greatly because it enables us to work with a fixed definition of the sky plane
for each WB, which we take to be the plane whose normal vector is the angular
bisector of the directions toward each of the stars in the WB. Since the sky
positions of the stars in our WB catalogue are expected to be particularly
precise, there should be negligible uncertainty associated with this
definition. Apart from this minor change, we follow the same method as that
described in Banik (2019) using mock parallax and proper motion data which
samples the _Gaia_ error distribution as described above. We assume that the
projected separation of each WB has only a negligible uncertainty, so we
obtain $r_{\rm{sky}}$ from the raw data.

####  2.3.1 The systemic RV

Compared to the detailed plan we set up earlier (Banik et al., 2021), a major
simplification has become possible with regards to the systemic RV of each WB
because we can generally obtain this directly from _Gaia_ DR3 (Gaia
Collaboration, 2023) without having to guess it from the proper motion and
distance to the Galactic disc plane. The RV of each star enters the
calculation of $\bm{v}_{\rm{sky}}$ through perspective effects (Shaya &
Olling, 2011; Banik, 2019; El-Badry, 2019) $-$ a receding WB with no internal
motion will have a shrinking apparent separation on the sky. The impact is
suppressed by the angular separation of the WB on our sky plane, but the
systemic RV still cannot be completely neglected: for a WB with
$r_{\rm{sky}}=10$ kAU, a typical systemic RV of 20 km/s translates into an
apparent proper motion of 10 m/s at a heliocentric distance of 100 pc. Since
the Newtonian $v_{c}$ of two Sun-like stars with a 10 kAU separation is 420
m/s, the impact on $\widetilde{v}$ would be $\la 0.03$ depending on the
geometry. While this is small, we opt to try and include the impact of
recessional motion because WBs can be even closer and/or have a larger RV.

We can see from the above that the at most km/s level difference in RV between
the stars in each WB would have only a negligible impact on $\widetilde{v}$.
This is because relative motions along the LOS and within the sky plane should
be comparable, but the impact of the former are suppressed by the angular WB
separation on the sky. Consequently, we limit ourselves to including only the
systemic RV of each WB. We find the systemic RV by taking an inverse variance
weighted mean of each star’s RV when this is available for both stars, leading
to an uncertainty smaller than either star’s RV uncertainty.444We combine the
RV measurements when these are known for both stars because the typical RV
uncertainty is around 1 km/s, which is much more than the expected WB relative
velocity. Thus, differences in the RV are best understood as due to random
errors rather than intrinsic WB orbital motion. If the RV is only available
for one star, we adopt this as the systemic value. We then assign both stars a
recessional velocity equal to the systemic RV, errors in which are propagated
as part of the Monte Carlo error propagation.

####  2.3.2 Reducing uncertainty in the relative distance

Another source of perspective effects (discussed further in Shaya & Olling,
2011) is the difference in heliocentric distances to the stars in each WB. The
relative LOS separation is important because proper motions must be multiplied
by distances to convert them from angular to physical velocities. The _Gaia_
parallaxes typically do not provide very tight constraints on the LOS
separation: while the astrometry is spectacular by historical standards,
trigonometric distances to stars 100 pc away are not accurate at the kAU level
(the typical uncertainty is about 80 kAU; see section 6.2 of BZ18, ). Using
the same arguments as in Section 2.3.1, a system with a tangential motion of
20 km/s but whose LOS separation is uncertain by 80 kAU would have an almost
80 m/s uncertainty in $\bm{v}_{\rm{sky}}$. This translates to an uncertainty
in $\widetilde{v}$ of up to 0.2 depending on the geometry, which could hamper
the WBT.

We intuitively expect that a WB’s 3D separation is not much larger than its
$r_{\rm{sky}}$, as suggested by the referee to Banik (2019). We exploit this
insight using the technique discussed in its section 2.3, which we briefly
summarize here. The mean of the heliocentric distances to the two stars is
taken to be the inverse variance weighted mean of the parallax distances to
each of the stars in the WB, causing the mean distance to vary slightly
between trials and thereby affect the distance modulus. Assuming a power-law
prior on the 3D separation $r$ of the form $P\left(r\right)\propto r^{-1.6}$
(Lépine & Bongiorno, 2007; Andrews et al., 2017; El-Badry & Rix, 2018), we
obtain a posterior inference on the ratio $x\equiv r/r_{\rm{sky}}$. Following
equation 18 of Banik (2019), we obtain that

$\displaystyle
P\left(x\right)~{}\propto~{}\frac{x^{-2.6}}{\sqrt{x^{2}-1}}\,,\quad x\geq 1.$
(8)

As part of the Monte Carlo error propagation, we then randomly sample this
distribution to obtain the separation of each WB along the LOS. We then push
one of the stars away from us and bring the other star towards us by half this
amount, with the displacement of each star being purely along its sky
direction at the reference epoch. We use a random coin toss in each Monte
Carlo trial to decide which star is to be pushed away from us. The revised
distance to each star alters its estimated absolute magnitude and thus its
mass.

This technique drastically reduces the uncertainty on the relative LOS
distance in systems with a very low $r_{\rm{sky}}$ (compare the black and blue
points in figure 3 of Banik, 2019). One disadvantage is that it implicitly
assumes that the WB is genuine, so our results are less robust against LOS
contamination. However, we mitigate this in our statistical analysis by
including an allowance for LOS contamination (Section 3.3).

Our approach is similar to that used in Nelson et al. (2021), though they used
a prior of $r^{-1.5}$ instead of $r^{-1.6}$. A more important difference is
that those authors were only interested in the 3D separation between the
stars, whereas we need to find the heliocentric velocity of each star, which
requires knowledge of the individual LOS distances. This causes our calculated
$\widetilde{v}$ parameter to depend on which star is closer to us. Since we do
not have this information, we use Monte Carlo trials to quantify the resulting
uncertainty. We stress that in our main analysis, the adopted value of
$\widetilde{v}$ is based on assigning the same LOS distance to both stars. The
purpose of our deprojection is to quantify the uncertainty and combine this
with other sources of uncertainty so that we can remove systems with
insufficiently precise $\widetilde{v}$.

The technique discussed above assigns the relative LOS distance a value of
zero and a non-Gaussian uncertainty of $\approx r_{\rm{sky}}$. This would only
reduce the uncertainties if the parallaxes were not accurate enough to
directly constrain the relative distance. Thus, we only apply the above
technique if $r_{\rm{sky}}$ is smaller than the uncertainty in the relative
LOS distance given by the _Gaia_ trigonometric parallaxes. While we generally
expect this to be the case, we do not enforce it. The 3D geometry of nearby
WBs may well be clear from existing _Gaia_ data, with the definition of
‘nearby’ extending to larger distances as the data improves. Our approach is
thus to use whichever method is expected to be more accurate. Note that if
_Gaia_ parallaxes are used directly, then the relative distance along the LOS
may have a non-zero mean value when averaged across different Monte Carlo
trials as these would propagate uncertainties in each parallax. In this case,
the nominal values adopted for our main analysis do not rely on the
deprojection algorithm outlined above because each star is assigned its
observed parallax distance.

####  2.3.3 Improved mass estimates

It is possible to slightly improve the accuracy of our mass estimates by using
the Final Luminosity, Age and Mass Estimator (FLAME; Pichon, 2007) work
package in _Gaia_ DR3. Masses obtained in this way are not available for all
stars, but when they are available, they should be much more precise because
they involve a detailed analysis of the spectrum. Figure 3 of Hernandez (2023)
shows that FLAME masses can differ by $\approx 0.05\,M_{\odot}$ from estimates
using a linear relation between absolute magnitude and the logarithmic mass.
Unfortunately, restricting to only WBs where both stars have a FLAME mass
would reduce the sample size too much for the WBT to be feasible. We therefore
apply a small correction to the masses estimated using our cubic fit to the
Pecaut & Mamajek (2013) mass-luminosity relation (Equation 6). Denoting these
masses by $M_{PM}$ and the revised masses with the FLAME calibration as
$M_{F}$, a good fit is given by

$\displaystyle M_{F}$ $\displaystyle=$ $\displaystyle
M_{PM}-0.07\tanh\left(\frac{M_{PM}-0.75\,M_{\odot}}{0.16\,M_{\odot}}\right)f\,,$
(9) $\displaystyle f$ $\displaystyle=$
$\displaystyle\begin{cases}\exp\left(2-\frac{M_{\odot}}{M_{PM}}\right),\quad\quad\textrm{
if}~{}M_{PM}\leq 0.5\,M_{\odot}\,,\\\ 1,\qquad\quad\quad\,\textrm{
if}~{}0.5\,M_{\odot}\leq M_{PM}\leq 0.75\,M_{\odot}\,,\\\
\exp\left(\frac{0.75\,M_{\odot}-M_{PM}}{2\,M_{\odot}}\right),\textrm{
if}~{}M_{PM}\geq 0.75\,M_{\odot}\,.\end{cases}$

The adjustment is tapered using the factor $f$ so it rapidly decays for stars
with mass $M<0.5\,M_{\odot}$ as there are no FLAME masses below this, so we
assume that $M_{PM}$ becomes very reliable at the low-mass end. We will see
later that our results are not much affected by restricting our WB sample to
systems where both stars have $M>0.5\,M_{\odot}$, which has the advantage that
FLAME results can serve as a calibration. $M_{PM}$ and FLAME masses gradually
converge when $M\gg 0.75\,M_{\odot}$.

Our small adjustment to $M_{PM}$ is designed to eliminate a small systematic
discrepancy with the FLAME mass in the stars for which the latter is
available. Even so, our calculated $M_{F}$ does not exactly coincide with the
FLAME mass. Based on the typical difference between the two mass estimates and
assuming that the FLAME masses are very accurate, we assign our estimated
masses a random Gaussian uncertainty of 0.024 dex or 5.5% at fixed absolute
magnitude. A 5.5% uncertainty is in line with the fact that masses estimated
directly from luminosities were expected to have a 6% uncertainty (Eker et
al., 2015; Mann et al., 2015), as argued in section 6.3 of BZ18. Note that our
mass estimates also include uncertainty in the absolute magnitude due to that
in the apparent magnitude and the trigonometric parallax.

###  2.4 Refined quality cuts for the WBT

The quality cuts discussed in Section 2.1 provide a reasonably carefully
prepared sample of WBs that can be used in further analyses. Unfortunately,
some of these WBs are unsuitable for our highly precise analysis. The
additional quality cuts in this work beyond those mentioned there are to
exclude:

  1. 1.

Stars in sky directions with a high total Galactic extinction towards an
extragalactic source (Section 2.4.1);

  2. 2.

Stars which are substantially below the main sequence on a colour-luminosity
diagram, which can make the mass estimate inaccurate (Section 2.4.2);

  3. 3.

Stars which sometimes appear as multiple peaks in the _Gaia_ images,
suggesting the presence of a CB companion that can sometimes be resolved
(Section 2.4.3);

  4. 4.

WBs where the RV is unknown for both stars (Section 2.4.4);

  5. 5.

WBs where the RV is known for both stars and there is a mismatch at $>3\sigma$
confidence (Section 2.4.5); and

  6. 6.

WBs whose $\widetilde{v}$ is too uncertain (Section 2.4.6).

We also restrict ourselves to WBs with $2<r_{\rm{sky}}/\rm{kAU}<30$ and
$\widetilde{v}<5$ for reasons discussed earlier. To reduce the computational
cost, we only consider WBs with a total mass of $0.464-4.31\,M_{\odot}$. These
cuts are now described further.

####  2.4.1 Dust extinction

As discussed in Section 2.2, our estimated mass for each star relies entirely
on its absolute $G$-band magnitude, which in turn is inferred from its
apparent magnitude and trigonometric parallax. The difference between apparent
and absolute magnitudes is assumed to arise purely from the distance modulus,
which implicitly assumes negligible dust extinction. This should be a good
assumption for stars within 250 pc and $>15^{\circ}$ from the Galactic disc,
especially as known star clusters are excluded (table 1 of PS23, ). Even so,
we can further reduce the impact of dust by considering only WBs where both
stars have a $V$-band extinction $A_{V}<0.5$ based on the dust maps of
Schlegel, Finkbeiner & Davis (1998), which reduces our sample size by 8.5%. We
note that their published values refer to the total Galactic extinction along
the LOS towards an extragalactic source. Since the WBs in our sample are all
within 250 pc, the extinction towards any of the stars in our sample should be
far smaller, making our cut on $A_{V}$ quite conservative.

####  2.4.2 Colour-magnitude diagram

Our mass estimates also rely on a mass-luminosity relation designed for main
sequence stars (Pecaut & Mamajek, 2013). To check if this is a good
assumption, we use Figure 2 to plot the relation between the absolute $G$-band
magnitude and the colour, which we define as the difference in apparent
magnitude between the _Gaia_ blue and red passbands. As expected, the vast
majority of the stars in our sample are on the main sequence. However, a small
number of stars are very far below the main sequence $-$ these are probably
white dwarfs (WDs). We remove these by excluding stars below the solid blue
line, which achieves a good separation between the main sequence and the WD
track. Only a very small number of WDs are removed by this cut and these are
all well separated from the main sequence, so any remaining WDs should have a
negligible impact on our results.

Figure 2: The colour-magnitude diagram of our sample just prior to removing 66
WBs where one star is below the solid blue line, which is designed to remove
the white dwarfs evident at the lower left. The two dashed grey lines have an
offset corresponding to a factor of two in luminosity. The lower line passes
through the main sequence, while the upper line passes through the secondary
track evident at higher luminosity. We attribute this double main sequence to
unresolved binaries where the stars have an exactly equal mass (Appendix A).

An interesting feature of our main WB sample’s colour-magnitude diagram is a
concentration of stars parallel to but brighter than the main sequence (Figure
2). This ‘double main sequence’ is caused by unresolved CBs where both stars
have an exactly equal mass, as demonstrated by the parallel dashed grey lines
with an offset corresponding to a factor of two in luminosity. We do not
exclude stars belonging to this double main sequence because our analysis
already accounts for it (Appendix A).

A major complication to the WBT is the presence of undetected CB companions
(Section 1). We follow a ‘mitigate + simulate’ strategy to minimize the
uncertainties introduced by CBs. This entails reducing the presence of CB
companions and modelling the distribution of those that remain. Since we do
not aim to have a sample that is completely free of stars with a CB companion,
we benefit from a much larger sample size than the study of Hernandez (2023).
We discuss below how we attempt to reduce the proportion of contaminated WBs.
Details of our CB model will be presented later (Section 3.2).

####  2.4.3 Multiple peaks in _Gaia_ images

Figure 3: Distribution of the _Gaia_ DR3 parameter
$\tt{ipd\\_frac\\_multi\\_peak}$ for the WBs in our sample prior to imposing
the quality cut related to this parameter. For each WB, we take the star with
the higher value. Although $\tt{ipd\\_frac\\_multi\\_peak}$ can be as high as
100, it is typically 0, as expected for a single source.

If there is a clearly detected CB companion to a star in a WB and it is
reasonably likely that the third star is bound to the WB, then the system is
excluded from our analysis (Section 2.1). However, CBs which are fainter
and/or on a tighter orbit might not be clearly resolved. The contaminated star
might then appear resolved into two sources in only some _Gaia_ focal plane
transits. We exploit this by using the parameter
$\tt{ipd\\_frac\\_multi\\_peak}$ in the _Gaia_ DR3 catalogue (see also section
5 of El-Badry, Rix & Heintz, 2021). This is the percentage of successful focal
plane transits in which a source is detected as multiple peaks in the image. A
genuinely isolated star would ideally give $\tt{ipd\\_frac\\_multi\\_peak}=0$,
but since it might spuriously show up as multiple peaks in the _Gaia_ images
due to issues like cosmic ray hits, we allow
$\tt{ipd\\_frac\\_multi\\_peak}\leq 2$ (a threshold of 2 was also used in
Pace, Erkal & Li, 2022). We reject any star that is detected as multiple peaks
more frequently.

Figure 3 shows the distribution of $\tt{ipd\\_frac\\_multi\\_peak}$ for the
WBs in our sample prior to imposing this quality cut, with the higher value
shown for the two stars in each WB. It is apparent that our quality cut does
not reduce the sample size very much. At the same time, it seems unlikely that
a star really has a marginally detected companion if it appears to be isolated
98% of the time. However, if a star frequently appears as multiple peaks in
the _Gaia_ images, then we might reasonably suspect it to be part of a CB.

####  2.4.4 Having a systemic RV

Even though we are implementing a 2D version of the WBT using only motions
within the sky plane, the systemic RV of each WB does still enter into our
analysis, albeit scaled down substantially by the angular separation of the WB
(Section 2.3.1). For our technique to work, we are forced to reject WBs where
the RV is unknown for both stars, which loses roughly 1/3 of our sample. We
accept systems where the RV is known for only one star because the systemic RV
is not needed very precisely and the RVs of the stars in a WB would differ by
at most a few km/s, especially for the more widely separated systems where
perspective effects are more important. Requiring both stars to have a
measured RV reduces the sample size from 8611 to 5504, but the proportion of
systems with $\widetilde{v}>2.5$ (which is a rough proxy for the level of
contamination) only drops from 10.2% to 9.5%.

####  2.4.5 Similar RVs when both are available

While we do not require both stars in a WB to have a known RV, this is often
the case. When available, we use this information to reject WBs whose two
stars have RVs that differ from each other by more than triple the quadrature
sum of the RV uncertainties or triple the Newtonian $v_{c}$, whichever is
larger.555Even if we had perfect data, the RVs of the stars in a WB could
still differ by order their Newtonian $v_{c}$. This reduces the sample size by
3.5% (we lose 308 WBs out of 8919 without this cut). Since almost 2/3 of the
systems in our WB sample have a reported RV for both stars, the fact that only
a small proportion of WBs consist of stars with discrepant RVs suggests that
our selected WBs are mostly genuine. The use of systemic RVs constrained to
km/s precision by _Gaia_ spectroscopy substantially reduces the scope for
perspective effects to alter the results obtained in this contribution.

####  2.4.6 The $\widetilde{v}$ uncertainty

The above quality cuts are designed to reduce systematic uncertainties in the
WB parameters. These also have random uncertainties, which we need to
quantify. As discussed in Section 2.3, we use Monte Carlo error propagation to
obtain the uncertainty in the all-important $\widetilde{v}$ using $2^{12}$
trials. We find that the uncertainty is very small for a large proportion of
the WBs in our sample (as expected from a simpler analysis of _Gaia_ DR2
results; Banik, 2019). We therefore apply a quality cut such that

$\displaystyle\text{Error in
}\,\widetilde{v}~{}\leq~{}0.1\max\left(1,\frac{\widetilde{v}}{2}\right)\,.$
(10)

We found that this achieves a good balance between the quality and quantity of
data. The uncertainty is allowed to be somewhat larger at higher
$\widetilde{v}$ because we use wider bins in $\widetilde{v}$ here when
comparing to theory (Section 3.5). Moreover, such systems are less crucial for
the WBT as they are not part of the main peak at $\widetilde{v}\la 2$,
implying significant contamination from some source.

Figure 4: The frequency distribution of the $\widetilde{v}$ uncertainty as
estimated from Monte Carlo error propagation of the $5\times 5$ _Gaia_ DR3
covariance matrix (Section 2.3). Different curves show results for different
ranges of $\widetilde{v}$, as indicated in the legend. We only consider WBs
whose $\widetilde{v}$ uncertainty is at most
$0.1\,\max\left(1,\widetilde{v}/2\right)$, truncating each error distribution
at the right. The truncation at the left is caused by the 0.024 dex mass
uncertainty at fixed $M_{G}$ imposing a minimum fractional uncertainty in
$\widetilde{v}$.

The distribution of $\widetilde{v}$ uncertainties is shown in Figure 4 for
different $\widetilde{v}$ ranges of our WB sample. As the WBT is mostly
concerned with the main peak to the $\widetilde{v}$ distribution at
$\widetilde{v}\la 2$, our cut implies a maximum $\widetilde{v}$ uncertainty of
only 0.1 in the most crucial parameter range. Since the mode of the main peak
is at $\widetilde{v}\approx 0.5$, measurement uncertainties that add an extra
0.1 in quadrature would only boost the root mean square value by about 2%,
which is negligible in comparison to the predicted MOND enhancement of 20%.
Moreover, the typical $\widetilde{v}$ uncertainty for such systems is much
smaller than the maximum permitted by our selection. Indeed, PS19 found the
$\widetilde{v}$ distribution in _Gaia_ DR2 changed little when restricting to
systems where the uncertainty is below 0.25 (compare their figures 12 and 13).
We therefore neglect $\widetilde{v}$ uncertainties in the rest of our
analysis, bearing in mind that our WB sample has already been chosen to have
an accurate $\widetilde{v}$.

We need to integrate WB orbits for a range of masses covering the mass range
of our sample. This is because larger mass systems have a larger MOND radius
(Equation 2), reducing the impact of MOND at fixed separation (BZ18). To limit
the computational cost, we restrict our analysis to WBs with a total mass of
$0.464-4.31\,M_{\odot}$. This reduces our sample size very slightly. Including
the remaining WBs would require a considerable broadening of the mass range
covered by our WB orbit library or a substantial worsening of its mass
resolution, but there would be very little gain in sample size. It would also
leave our results more vulnerable to systematic issues with very low or very
high mass stars.

###  2.5 The WB sample

After applying the quality cuts described so far, we are left with a final
sample of 8611 WBs. Their mass distribution is shown in Figure 5, with blue
bars showing individual stars and red bars showing binary total masses. The
latter are more relevant physically, so we will use $M$ to mean the binary
total mass and $M_{\star}$ to denote the mass of a single star $-$ unless it
is clear that the discussion refers to stars considered individually.

Figure 5: The mass distribution of all the stars in our WB sample (blue bars).
The red bars show the distribution of the total mass of each binary. The outer
edges of the outermost bins show the range of the data.

As with previous detailed investigations into the WBT (; ), we perform it by
decomposing the data into pixels in the space of $r_{\rm{sky}}$ and
$\widetilde{v}$ (Figure 6). It is clear that there are many more systems at
low $r_{\rm{sky}}$ than at high $r_{\rm{sky}}$, in line with earlier results
(e.g. Andrews et al., 2017). The main peak to the $\widetilde{v}$ distribution
is quite prominent, with the mode located at $\widetilde{v}\approx 0.5$. There
is also an extended tail going out to much larger $\widetilde{v}$, which was
previously reported (PS19) and attributed to undetected CBs (Belokurov et al.,
2020; Clarke, 2020). Interestingly, a gap appears in the distribution at
$\widetilde{v}\approx 2$, but only at high $r_{\rm{sky}}$. We attribute this
to the minimum in the $\widetilde{v}$ distribution at $\widetilde{v}\approx
1.5$, which is evident also at lower $r_{\rm{sky}}$. This minimum is caused by
WBs having a sharp peak to their $\widetilde{v}$ distribution at
$\widetilde{v}\approx 0.5$ followed by a rapid decline as a consequence of
orbital mechanics and projection effects (e.g. ; , ). There is an extended
tail due to CB contamination, but this too would generally lead to a declining
$\widetilde{v}$ distribution because a slower CB orbital velocity corresponds
to a wider CB separation, which is less likely (a more detailed explanation is
given in Banik et al., 2021). However, we expect an increasing contribution
from chance alignments at larger $\widetilde{v}$ because their distribution is
expected to be linear in $\widetilde{v}$ (see Equation 27). The minimum in the
$\widetilde{v}$ distribution is not very apparent at $r_{\rm{sky}}\la 5$ kAU,
where WBs are very common and chance alignments are relatively unimportant.
The minimum becomes readily apparent at larger $r_{\rm{sky}}$ where WBs are
quite rare, increasing the relative contribution of chance alignments because
they should have a flat distribution with respect to $r_{\rm{sky}}$ (Equation
27). The fact that WBs are much less common at high $r_{\rm{sky}}$ naturally
leads to empty pixels around the minimum of the $\widetilde{v}$ distribution,
with the $\widetilde{v}$ range covered by these empty pixels expected to
gradually widen with increasing $r_{\rm{sky}}$.

Figure 6: The number of WBs in different pixels of
$\left(r_{\rm{sky}},\widetilde{v}\right)$. Results are shown binned according
to our pixellation scheme, which was fixed prior to any detailed analysis (the
pixel widths increase gradually towards higher values). The comparison with
theory is based on the likelihood of matching the number counts shown here,
which we find by multiplying together the binomial likelihood for each pixel
(Section 3.5). Notice the lack of WBs at $\widetilde{v}\approx 2$ for the
highest $r_{\rm{sky}}$ (empty pixels are shown in white). We interpret this as
due to dilution of the CB contamination over a wide $\widetilde{v}$ range due
to the low Newtonian $v_{c}$ of the WB. There is also little LOS contamination
here: it only becomes important when $r_{\rm{sky}}$ and $\widetilde{v}$ are
both large (Section 3.3).

It is clear that we need a detailed model including CB and LOS contamination
if we are to infer the law of gravity at low accelerations and search for the
predicted MOND signal. This is the subject of the next section.

##  3 Modelling the WB dataset

The $\left(r_{\rm{sky}},\widetilde{v}\right)$ distribution of our WB sample is
a consequence of several simultaneously occurring physical processes and
projection effects, which we try to model in this section to best reproduce
Figure 6. In addition to WB orbital motion (Section 3.1), we also consider how
an undetected CB population affects the observables (Section 3.2) and take
into account the possibility of some LOS contamination (Section 3.3). This
might be field stars, but it could also be from stars that formed in the same
star cluster. Given the results of PS19 and the excellent quality of _Gaia_
DR3, we expect LOS contamination to only slightly affect our WB sample and to
become important only at large separation and $\widetilde{v}$.

LOS contamination can be included in our model fairly easily. Moreover,
detailed WB orbit modelling has previously been conducted in Newtonian and
Milgromian gravity in a manner that can readily be repurposed for the WBT
(BZ18). As a result, the main complication that we need to handle is the
possibility of an undetected CB companion around 0, 1, or 2 of the stars in a
WB. Since recoil motion and additional mass created by an undetected CB
physically affect the WB rather than merely contribute an extra population to
the statistics, including CB contamination greatly increases the complexity
and computational cost of the WBT. We were eventually able to devise a
detailed plan that captures the essential physics in a computationally
feasible way (Banik et al., 2021). We refer the reader to that work for a more
detailed look at the computational techniques needed to implement the WBT in a
reasonable timeframe. This allows us to prepare the first detailed model of
the WB population whose parameter space is thoroughly explored using standard
statistical techniques.

###  3.1 The WB population

We set up a library of WB orbits using the prior work of BZ18 with minimal
adjustments, which we briefly describe. We need to run calculations over a
wider mass range in order to cover the range of WB masses in our sample (red
bars in Figure 5). We also found that the computational cost of including CB
contamination can be greatly reduced if instead of recording results with
respect to the projected separation $r_{\rm{sky}}$, we use

$\displaystyle\widetilde{r}_{\rm{sky}}~{}\equiv~{}\frac{r_{\rm{sky}}}{a}\,,$
(11)

where $a$ is the semi-major axis of the WB (this is motivated in Section
3.2.3). An incidental benefit of this approach is that at low $a$, recording
probabilities into different bins in $\widetilde{r}_{\rm{sky}}$ rather than
$r_{\rm{sky}}$ improves the resolution in $r_{\rm{sky}}$. Since we record
results separately for different $a$ and only marginalize over $a$ at a much
later stage, we can recover $r_{\rm{sky}}\equiv a\widetilde{r}_{\rm{sky}}$ as
and when needed, which is necessary for a comparison with observations because
$a$ is generally not known.

In MOND, the Galactic EFE plays a crucial role (BZ18). This can be estimated
as $g_{e}=v_{c,\odot}^{2}/R_{0}$, where $R_{0}$ is the Galactocentric distance
of the Sun and $v_{c,\odot}$ is the Galactic rotation curve amplitude at the
Solar circle (sometimes called the Local Standard of Rest or LSR). As in
section 3.6 of BZ18, we assume that $v_{c,\odot}=232.8$ km/s and $R_{0}=8.2$
kpc (McMillan, 2017), which implies that $g_{e}=1.785\,a_{{}_{0}}$. Subsequent
studies show that $R_{0}$ must be very close to 8.2 kpc based on combining
astrometry and spectroscopy of the star S2 near the Galactic centre black hole
(Gravity Collaboration, 2019). Since its proper motion is precisely known
(Reid & Brunthaler, 2004; Gordon, de Witt & Jacobs, 2023), the Galactocentric
tangential velocity of the Sun now has very little uncertainty, constraining
the LSR velocity to very close to 233 km/s (McGaugh, 2018; Zhou et al., 2023)
if the non-circular motion of the Sun is taken to be precisely known from
Schönrich, Binney & Dehnen (2010). The latter is somewhat more uncertain than
the other parameters entering into the LSR velocity (Schönrich, 2012; Francis
& Anderson, 2014; Bovy et al., 2015), but it is still clear that $v_{c,\odot}$
is very unlikely to differ from our assumed value by more than 10 km/s, which
implies that uncertainty in $v_{c,\odot}$ has little effect on the predicted
MOND signal in local WBs (Section 5.3). Moreover, $g_{e}$ has recently been
directly determined by the measurement of the Solar System’s acceleration
relative to distant quasars (Gaia Collaboration, 2021b). Their work discovered
the changing aberration angle of distant quasars due to the changing velocity
of the Solar System barycentre. Its directly measured acceleration points
directly towards the Galactic centre within a small uncertainty (see their
figure 10) and has a magnitude of $g_{e}=1.94\pm 0.13\,a_{{}_{0}}$. Since
gravitational fields from objects beyond the Galaxy are thought to be
negligible in comparison to the Galactic gravitational field, this provides a
completely independent and very direct confirmation of the kinematic estimate
of $\bm{g}_{e}$ used in BZ18 and in this contribution.

####  3.1.1 Orbit modelling

WBs are bound together by their mutual gravity. This is easy to calculate in
the Newtonian case but is somewhat more complicated in MOND, as discussed
next. Our calculation of the Milgromian gravitational field and integration of
WB orbits in this field use the same approach as BZ18, to which we refer the
reader for further details.

In an isolated spherically symmetric mass distribution, the asymptotic MOND
behaviour is given by Equation 1. To interpolate between the Newtonian and
deep-MOND regimes, we have to use an interpolating function $\nu$ with
argument $g_{{}_{N}}$ such that $\bm{g}=\nu\bm{g}_{{}_{N}}$. In the
quasilinear formulation of MOND (QUMOND; Milgrom, 2010), this can be
generalized to a more complicated geometry by taking the divergence of both
sides.

$\displaystyle\nabla\cdot\bm{g}~{}=~{}\nabla\cdot\left(\nu\bm{g}_{{}_{N}}\right).$
(12)

In this contribution, we adopt the simple interpolating function (Famaey &
Binney, 2005) for reasons discussed in section 7.1 of BZ18.

$\displaystyle\nu~{}=~{}\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{a_{{}_{0}}}{g_{{}_{N}}}}\,.$
(13)

The WBT depends significantly on the adopted interpolating function because
local WBs have an intermediate total acceleration. Using a very sharp
transition can cause WBs to be completely Newtonian even if MOND is correct.
However, this would not be consistent with the RAR (we discuss this further in
Section 5.3).

Since local WBs are significantly affected by the Galactic EFE, we need to
include the Galactic $\bm{g}_{e}$. This is assumed to point directly towards
the Galactic centre and to have a magnitude of $1.785\,a_{{}_{0}}$ based on
the Galactic rotation curve, as discussed above. In QUMOND, we need the
Newtonian-equivalent external field $\bm{g}_{{}_{N,e}}$, which is what the
Galactic gravity on the Solar neighbourhood would have been in Newtonian
gravity without DM. As discussed in section 9.3.1 of BZ18, we can use the
spherically symmetric relation between $\bm{g}_{e}$ and $\bm{g}_{{}_{N,e}}$
because the Sun is many disc scale lengths from the Galactic centre, reducing
the importance of the Galaxy’s disc geometry.

$\displaystyle\bm{g}_{e}\,\nu\left(g_{e}\right)~{}=~{}\bm{g}_{{}_{N,e}}\,.$
(14)

Solving this implicitly gives $g_{{}_{N,e}}=1.144\,a_{{}_{0}}$ in the Solar
neighbourhood. The Newtonian gravity in a WB is then found by adding
$\bm{g}_{{}_{N,e}}$ to the contribution from the WB itself, which is called
$\bm{g}_{{}_{N,i}}$ in some works as it is internal to the WB.

To limit the computational cost, we treat each WB as a test particle orbiting
a point mass containing the binary total mass $M$. This approach is valid in
Newtonian gravity, but it was shown in section 7.3 of BZ18 that it remains a
very good approximation in MOND. This is because the fairly strong Galactic
EFE causes local WBs to lack an extended regime in which they are isolated and
at low accelerations. For all intents and purposes, WBs are either in the
Newtonian regime or they are dominated by the EFE. In either case, MOND
becomes linear in the mass distribution, allowing the potentials of the two
stars to be superposed.

With this approximation, the problem becomes axisymmetric about $\bm{g}_{e}$.
We then solve Equation 12 using the ‘ring library’ protocol discussed in BZ18.
This involves finding the source term
$\nabla\cdot\left(\nu\bm{g}_{{}_{N}}\right)$ at every point in a 2D grid and
then using direct summation to find $\bm{g}$. The calculations are accelerated
at off-axis locations by scaling the previously calculated Newtonian gravity
of a unit mass and unit radius ring, avoiding the need to sum contributions to
$\bm{g}$ from different positions around each ring. An analytic allowance is
made for $\nabla\cdot\left(\nu\bm{g}_{{}_{N}}\right)$ being non-zero outside
the rather large region used in the integration (the extra contribution to the
potential is given in equation 9 of BZ18; see their appendix A for a
derivation).

Figure 7: The azimuthally averaged boost factor $\eta$ to the radially inward
gravity compared to the Newtonian expectation for WBs with different
separation relative to their MOND radius, assuming the simple interpolating
function (Equation 13). The solid red (blue) horizontal line shows the
asymptotic analytic expectation in QUMOND (AQUAL), which we give in Equation
37 (38). The numerical QUMOND result central to our study (solid black line;
BZ18, ) has the expected asymptotic behaviour. The dashed blue line shows
equation 15 of Chae & Milgrom (2022), which has the wrong asymptotic limit
(mismatch between blue lines towards the right) and is based on numerical
calculations that only deal with a weak external field rather than a dominant
one (see the text). The red histogram shows the distribution of
$r_{\rm{sky}}/r_{{}_{M}}$ for the WBs in our sample. The outer edges of the
most extreme bins show the full range of the dataset.

The Milgromian potential of a point mass is not spherically symmetric in the
presence of the EFE (Milgrom, 1986a, 2010; Banik & Zhao, 2018a). This has some
interesting consequences (Candlish et al., 2018; Thomas et al., 2018; Banik &
Kroupa, 2019b; Banik et al., 2020; Kroupa et al., 2022). Since the WB
separation $\bm{r}$ is expected to sample a range of directions relative to
the EFE, the WBT is mostly sensitive to the angle-averaged inward radial
gravity, which we denote $\langle g_{r}\rangle$. Since the use of
$\widetilde{v}$ rather than actual velocities effectively divides out the
Newtonian prediction $g_{{}_{N}}$, the WBT is really about the parameter

$\displaystyle\eta~{}\equiv~{}\frac{\langle g_{r}\rangle}{g_{{}_{N}}}\,,$ (15)

whose square root is probed by velocity measurements. We use Figure 7 to show
$\eta$ as a function of $r/r_{{}_{M}}$ (solid black line).666Zonoozi et al.
(2021) provide an analytic fit to these results for arbitrary $g_{e}$ with the
correct asymptotic behaviour (see their equations 23 and 24). There is no
boost to the Newtonian expectation when $r\ll r_{{}_{M}}$, but this changes
rapidly until $r\approx r_{{}_{M}}$, beyond which the Milgromian enhancement
to $\langle g_{r}\rangle$ saturates at just over 40% due to the EFE. Without
it, the boost factor would continue growing linearly with $r$. In the
asymptotic or quasi-Newtonian regime (discussed further in section 3.3 of
Asencio et al., 2022) where the enhancement factor has reached its maximum,
the WB gravity can be treated as a small perturbation about the Galactic
gravity, making the potential analytic (equation 30 of Banik et al., 2020).
The analytic expectation is shown using a horizontal red line, revealing an
excellent match to the numerical result when $r>r_{{}_{M}}$.

In this contribution, we focus on the QUMOND formulation because it is more
computer-friendly than the earlier aquadratic Lagrangian formulation (AQUAL;
Bekenstein & Milgrom, 1984). The asymptotic AQUAL result is shown in Figure 7
using a horizontal blue line (equation 35 of BZ18, ). As discussed there, this
result is very similar to the QUMOND result for the same interpolating
function, i.e. the same relation between $\bm{g}$ and $\bm{g}_{{}_{N}}$ in
spherical symmetry. Numerical AQUAL results are not available in the
transition zone where the EFE is neither negligible nor dominant, but we can
get a general idea using equation 15 of Chae & Milgrom (2022), which is a fit
to numerical results for a weak EFE as relevant to galaxy rotation curves. We
show this using a dashed blue line. It is clear that this equation has the
wrong asymptotic limit, as can also be seen from the equation: when the EFE
dominates, it takes the argument $1.1\,g_{{}_{N,e}}$ rather than
$g_{{}_{N,e}}$. Even if this issue is fixed, the functional form is not
compatible with the EFE-dominated point mass potential in AQUAL (equation 66
of Milgrom, 2010) or the resulting radial gravity (equation 32 of BZ18, ).
Despite these issues when the EFE is strong, equation 15 of Chae & Milgrom
(2022) should still give a reasonable idea of the AQUAL enhancement to
$\langle g_{r}\rangle$ when $r\ll r_{{}_{M}}$ because it is a fit to numerical
results in this regime. We see that when $r\la 0.5\,r_{{}_{M}}$, the AQUAL and
QUMOND results are almost identical. Therefore, results for the WBT are not
sensitive to which formulation of MOND we use. We discuss this further in
Section 5.3.

Another aspect of Figure 7 is the distribution of $r_{\rm{sky}}/r_{{}_{M}}$,
which we show using red bars (the left edge of the first bar shows the minimum
value). Our sample covers down to $r_{\rm{sky}}/r_{{}_{M}}=0.14$, so we
consider it to be reasonably sensitive down to $r/r_{{}_{M}}=0.2$, leaving
some allowance for 3D separations exceeding $r_{\rm{sky}}$. A WB separated by
only $r_{{}_{M}}/5$ has an acceleration of $g_{{}_{N}}=25\,a_{{}_{0}}$ and can
thus be considered isolated. In this case, $\nu=1/2+\sqrt{1/4+1/25}=1.04$,
which is very small compared to the predicted enhancement to $g_{{}_{N}}$ for
more widely separated WBs. This is even truer for interpolating functions that
have a somewhat faster return to Newtonian behaviour at high accelerations in
order to better match Solar System ephemerides (Hees et al., 2014, 2016). At
the other end, our sample probes well into the asymptotic regime where we
would expect almost the maximum possible enhancement to $g_{{}_{N}}$ given the
Galactic EFE. Our sample thus covers a wide enough range of internal WB
accelerations to allow a sensitive search for the predicted MOND effect. While
our WBs experience a narrower range of total gravitational field strengths
because the Galactic EFE essentially imposes a lower limit, we argue later
that this is not a major limitation given rotation curve constraints on the
MOND interpolating function in the transition zone (Section 5.3).

We use the assumed gravity law to integrate WB orbits with a range of mass
$M$, semi-major axis $a$, and orbital eccentricity $e$, defining $a$ and $e$
for a generalized gravity theory as in section 2.3.1 of BZ18. Following their
approach, each orbit is started at a separation of $a$ and integrated for 20
revolutions. This should be more than sufficient in the Newtonian case, but it
is also adequate in MOND because the resulting statistical distribution is
similar to that yielded by a much longer 5 Gyr integration (see their figure
15). We also consider a range of inclinations $i$ between the orbital pole and
the Galactic centre direction defining the EFE. Results for different
inclinations are marginalized assuming a prior of $\sin i$, which is
appropriate for an isotropic distribution of orbital poles. Orbits are
truncated if the stars get to within 50 AU or become more distant than 100 kAU
as we consider such WBs to have crashed together or to have likely become
unbound at some point by perturbations from passing stars, respectively
(section 2.3.2 of BZ18, ). At each timestep, we consider a dense 2D grid of
viewing directions and compute the
$\left(\widetilde{r}_{\rm{sky}},\widetilde{v}\right)$ inferred by a distant
observer along the chosen LOS, with results statistically weighted by the
solid angle covered by each LOS. We also weight the results at different
timesteps by the duration of the timestep, which varies because we use an
adaptive timestep to reduce the computational cost. In this way, we build up a
simulated $\left(\widetilde{r}_{\rm{sky}},\widetilde{v}\right)$ distribution
for different $\left(M,a,e\right)$.

In the Newtonian case, it is not necessary to consider different $M$ or $i$,
but we use a $5\times$ higher resolution in $e$ because the $\widetilde{v}$
distribution at fixed $e$ has sharp peaks corresponding to pericentre and
apocentre and $e$ does not change in a Newtonian orbit integration (unlike in
MOND; see Paučo & Klačka, 2016). To get a smooth distribution when
marginalizing over different $e$, BZ18 found that a rather high resolution is
required in $e$ for the Newtonian case. This is not needed in MOND because the
peaks are blurred out by changes in the shape of the orbit over time:
typically about 8 orbits are needed to go from nearly circular to very
eccentric back to nearly circular (see their figure 20).

####  3.1.2 Orbital parameter distribution

The mass distribution of WBs is obtained directly from that of our WB sample
(Figure 5). We count up how many WBs are in the mass range covered by each
value of $M$ at which we perform an orbit integration. This directly
determines the relative statistical weights of different $M$ values in our WB
orbit library.

Following Andrews et al. (2017), we assume that the semi-major axes of the WBs
follow a broken power law in $a$. The logarithmic slope at low $a$ is fixed to
$-1.6$, which is known to be valid down to at least 0.5 kAU (El-Badry & Rix,
2018). While this can in principle extend up to arbitrarily large $a$, we
assume that it breaks down at $a=a_{\rm{break}}$, beyond which the logarithmic
slope becomes $\beta<-1$ to ensure a convergent integral. Our assumed
distribution of $a$ is therefore given by

$\displaystyle
P\left(a\right)~{}\propto~{}\begin{cases}a^{-1.6}\,,&\textrm{if}~{}a<a_{\rm{break}}\,,\\\
a^{\beta}\,,&\textrm{if}~{}a>a_{\rm{break}}\,,\end{cases}$ (16)

where $a_{\rm{break}}$ and $\beta$ are free parameters that our analysis
infers from the data. Note that since we are modelling the distribution of our
sample of WBs rather than the full WB population, the inferred values of these
parameters may be somewhat biased. In particular, the requirement to have
accurate $\widetilde{v}$ measurements is quite difficult to satisfy at large
separations because the Newtonian $v_{c}$ is low. This is likely to reduce
$\beta$.

Following Hwang et al. (2022), we parametrize the distribution of orbital
eccentricities using a power law with index $\gamma$.

$\displaystyle P\left(e\right)~{}=~{}\left(\gamma+1\right)e^{\gamma}\,.$ (17)

$\gamma$ is left as a free parameter that is allowed to vary over the range
$0-4$ in all our analyses. We expect a nearly thermal distribution of
eccentricities ($\gamma=1$; Jeans, 1919; Ambartsumian, 1937) $-$ this is
discussed in more detail in section 4.2 of Kroupa (2008). We use a flat prior
on $\gamma$ for our nominal analysis. However, Hwang et al. (2022) found that
a slightly superthermal distribution ($\gamma>1$) better fits the observed
distribution of the angle $\psi$ between $\bm{r}_{\rm{sky}}$ and
$\bm{v}_{\rm{sky}}$ (this is related to the WB orthogonality test outlined in
section 8.2.1 of BZ18, ). We therefore run a variation of our nominal analysis
in which we assume a Gaussian prior of $\gamma=1.32\pm 0.09$, which should
shrink the uncertainties somewhat. We note that since CB contamination would
tend to flatten out the distribution of $\psi$ and thereby drive the inferred
eccentricity distribution towards thermal, it is quite likely that $\gamma$ is
slightly higher than inferred by Hwang et al. (2022) given that those authors
do not consider undetected CBs. On the other hand, their figure 10 shows that
$\gamma$ can be overestimated slightly at large separations for an underlying
slightly superthermal distribution.

####  3.1.3 Interpolating the gravity law

An important aspect of our analysis is that we consider an arbitrary gravity
law on a sliding scale between Newtonian and Milgromian. For this purpose, we
introduce the gravity law parameter $\alpha_{\rm{grav}}$, which is the most
important parameter for the WBT. We define this such that
$\alpha_{\rm{grav}}=0$ represents Newtonian gravity and $\alpha_{\rm{grav}}=1$
represents QUMOND with the simple interpolating function. We allow
$\alpha_{\rm{grav}}$ values somewhat outside this range to capture the
possibility of behaviour different to either theory. The most plausible
outcomes to the WBT are a strong preference for $\alpha_{\rm{grav}}$ values
very close to 0 or 1, with the large sample size hopefully tightening the
error budget enough that both possibilities are not simultaneously consistent
with the observations (a few thousand WBs should be sufficient; see BZ18, ).
Since the WBT will have some uncertainties, we allow $\alpha_{\rm{grav}}$ to
lie in the range $\left(-2,3.6\right)$. Negative values imply that the gravity
in the asymptotic regime has a lower than Newtonian normalization, while
$\alpha_{\rm{grav}}>1$ implies a normalization higher than expected in MOND.

To set up a parametrized gravity law, we obtain the simulated
$\left(\widetilde{r}_{\rm{sky}},\widetilde{v}\right)$ distribution in both
gravity theories for each $\left(M,a,\gamma\right)$, with different orbital
eccentricities marginalized over for the adopted $\gamma$ using the
eccentricity distribution it defines (Equation 17). The main idea of our
parametrization is to transform the Newtonian
$\left(\widetilde{r}_{\rm{sky}},\widetilde{v}\right)$ distribution into the
Milgromian one and then interpolate the coefficients involved in the
transformation so that arbitrary gravity laws can be considered.777We
marginalize over $M$ and $a$ at a later stage. To simplify our analysis, we
consider a simple 2D stretch in which the probability distribution is
stretched along the $\widetilde{r}_{\rm{sky}}$ axis by a factor $S_{r}$ and
then along the $\widetilde{v}$ axis by a factor $S_{v}$, with $S_{r}\neq
S_{v}$ in general. These are found by running a 2D gradient descent in
$\left(S_{r},S_{v}\right)$ to minimize the sum of squared residuals between
the stretched Newtonian array and the MOND array (Fletcher & Powell, 1963).
While a simple 2D stretch applied to the former does provide a good
approximation to the latter, there are still some residuals, which we
calculate for each pixel. This array of residuals is itself scaled by a factor
of $\alpha_{\rm{grav}}$ and stretched in 2D by some factors depending on
$\alpha_{\rm{grav}}$ and $\left(S_{r},S_{v}\right)$. The idea behind applying
this ‘corrections array’ is to recover the Newtonian (Milgromian)
$\left(\widetilde{r}_{\rm{sky}},\widetilde{v}\right)$ distribution exactly if
$\alpha_{\rm{grav}}=0$ (1). If we consider an intermediate case where
$\alpha_{\rm{grav}}=0.5$, we would expect the $\widetilde{v}$ distribution to
stop at a value about halfway between the maximum which arises in the
Newtonian and MOND distributions. To prevent the corrections array from
creating non-zero contributions all the way out to the maximum $\widetilde{v}$
which arises in MOND, the residuals array itself needs to be stretched
somewhat, though in this case we expect $S_{v}<1$. Since no corrections are
required if $\alpha_{\rm{grav}}=0$, we need to avoid applying any corrections
in this case, which is achieved by our scaling factor of $\alpha_{\rm{grav}}$.
We use linear interpolation in all the stretch and scaling factors discussed
above with respect to $\alpha_{\rm{grav}}$. To avoid negative probabilities
due to negative values in the corrections array, any negative entries in the
resulting $\left(\widetilde{r}_{\rm{sky}},\widetilde{v}\right)$ distribution
are converted to zero and the array renormalized.

###  3.2 Undetected companions

As discussed in Section 1.1, undetected companions are the main source of
contamination to the WBT given the low measurement uncertainties (Figure 4)
and the low LOS contamination fraction evident from the low number counts in
the rather wide pixels at high $r_{\rm{sky}}$ and $\widetilde{v}$ (Figure 6).
In particular, the declining $\widetilde{v}$ distribution it reveals is a
clear sign of some process that broadens the WB distribution, which we
otherwise expect to be concentrated at $\widetilde{v}\la 2$ for any plausible
gravity theory (PS19).

We therefore need to include the possibility that a star in a WB is actually
itself part of a CB, even though the CB companion is undetected. Including CBs
proved to be by far the most challenging aspect of the WBT, so we had to make
some simplifying assumptions to keep the complexity and computational cost
manageable. Due to their much smaller separation, we assume that the CB is
completely Newtonian internally, even if MOND affects the motion of its
barycentre around the other star in the WB.888Bekenstein & Milgrom (1984)
showed that a body with high internal accelerations can still exhibit MONDian
behaviour for its barycentre. It will become clear later that the undetected
CBs relevant to the WBT indeed have only a small separation compared to the
WB, justifying this approximation.

####  3.2.1 Impact on WB observables

Undetected companions have two main effects on a WB: increased mass and
induced recoil velocity. Since the CB companion is expected to tightly orbit
the star it contaminates, we neglect the impact of the CB on the WB’s
$r_{\rm{sky}}$. For this reason, we also assume that the light emitted by the
‘undetected’ star is blended with that of the contaminated star, i.e. they are
unresolved by _Gaia_. This is reasonable given that resolved third star
companions would cause the WB to be removed from our sample at an early stage
(Section 2.1). Indeed, we will see later that our analysis prefers to avoid
CBs whose separation is significant compared to that of the WB.

The actual mass of the CB exceeds its estimated mass by a factor

$\displaystyle\Delta\widetilde{M}~{}\equiv~{}\frac{\text{Total CB
mass}}{\text{Estimated CB mass}}\,-\,1\,.$ (18)

This is positive because the mass-luminosity relation is very steep (Figure
1), so blended light from the ‘undetected’ companion raises the estimated CB
mass by less than the companion’s mass. As a result, two unresolved Sun-like
stars would appear to have a total mass of only $1.19\,M_{\odot}$ even though
their total mass is actually $2\,M_{\odot}$, yielding
$\Delta\widetilde{M}=0.68$. Since the mass-luminosity relation is not an exact
power law, $\Delta\widetilde{M}$ depends on the masses of the stars involved.
We handle this complication by considering as a prior the mass distribution of
all the stars in our WB sample (blue bars in Figure 5). For a fixed fraction
$\widetilde{q}\leq 1/2$ of the CB total mass in the undetected companion, we
consider a dense grid of primary star masses covering the range shown
there.999‘Primary’ here refers to the more massive star in the CB, which may
itself be a sub-dominant component of the total WB mass. In each case, we use
the mass-luminosity relation from Pecaut & Mamajek (2013) to obtain the
fraction $\widetilde{L}$ of the CB’s total luminosity contributed by the
undetected companion $-$ this will be important later.101010For simplicity, we
do not implement the FLAME calibration discussed in Section 2.3.3 when
converting masses to luminosities or vice versa in this section. We then
numerically invert Equation 6 using the Newton-Raphson algorithm to find the
estimated mass of the CB given the total _Gaia_ -band luminosity of the two
stars in it, thus following the same approach as observers would use without
knowing about the CB. In this way, we obtain the relation between
$\Delta\widetilde{M}$ and $\widetilde{q}$, which we show as the solid red line
in Figure 8. We also use a dotted red line to show the results if we neglect
the blended light contributed by the undetected companion. The results are
similar either way unless the CB consists of two stars with an almost equal
mass: in this limit, blended light reduces $\Delta\widetilde{M}$ from 1 to
0.72.

The impact of blended light is more significant when it comes to the recoil
velocity induced by an undetected companion. In the limit that the CB has two
equal mass stars such that $\widetilde{q}=1/2$, the photocentre and barycentre
would coincide and move only due to the WB orbit. This is true regardless of
the relative velocity $\bm{v}$ of the CB orbit. In general, the velocity of
the primary star relative to the CB barycentre would be $\widetilde{q}\bm{v}$,
while the undetected star would move at $-\left(1-\widetilde{q}\right)\bm{v}$,
which is in the opposite direction. Once these velocities are averaged
weighted by the mass fractions, the CB barycentre is of course static. But if
we instead weight the stars using their luminosities, we get that the
photocentre velocity is

$\displaystyle\bm{v}_{\rm{phot}}~{}=~{}\left[\widetilde{q}\left(1-\widetilde{L}\right)-\widetilde{L}\left(1-\widetilde{q}\right)\right]\bm{v}~{}=~{}\left(\widetilde{q}-\widetilde{L}\right)\bm{v}\,.$
(19)

Since we are only interested in the magnitude of $\bm{v}_{\rm{phot}}$, the
impact on the observed WB relative velocity (based on photocentres) is
$\propto\left|\widetilde{q}-\widetilde{L}\right|$ (as also found by Penoyre,
Belokurov & Evans, 2022). Moreover, $\bm{v}$ is itself proportional to the
Keplerian velocity of the CB, which depends on its total mass and separation.
Finally, the impact on $\widetilde{\bm{v}}$ also needs to take into account
that blended light increases the apparent mass of the WB and thus its
Newtonian $v_{c}$ (Equation 3), but the extra mass is underestimated thanks to
the steep mass-luminosity relation. We assume for simplicity that both stars
in a WB are affected by CBs with a similar mass ratio to their primary star.

Figure 8: How the fraction of a CB’s total mass in its undetected component
affects the offset between its photocentre and barycentre (black lines) and
the fraction by which the total mass exceeds the mass estimated by observers
based on the luminosity (red lines). The scale for the left $y$-axis is
arbitrary: the important aspect is the functional form. The dot-dashed lines
neglect the increase to the estimated mass caused by the luminosity of the
undetected component. The solid lines include this effect, which helps to
somewhat mitigate the impact of CB contamination. The solid black line is
proportional to the CB’s impact on the WB $\widetilde{v}$ under the assumption
that an undetected CB increases the estimated mass of both WB components by a
similar fraction. Only the solid lines shown here enter into our model as
$\Delta\widetilde{M}$ (red; Equation 18) and $f_{\rm{pb}}$ (black; Equation
20).

Putting all this together, we get that the impact on $\widetilde{\bm{v}}$
scales with the parameter

$\displaystyle
f_{\rm{pb}}~{}=~{}\left|\widetilde{q}-\widetilde{L}\right|\overbrace{\sqrt{1+\Delta\widetilde{M}}}^{f_{M}}\,.$
(20)

This is shown as the solid black line in Figure 8. The dotted black line shows
the result if $\Delta\widetilde{M}$ is calculated neglecting blended light. It
is clear that $f_{\rm{pb}}$ is not sensitive to details of how
$\Delta\widetilde{M}$ is calculated for the purposes of the above equation.
This is because if $\widetilde{q}\ll 1/2$ and the CB has an extreme mass
ratio, then the undetected star contributes almost no blended light, making
$\Delta\widetilde{M}$ insensitive to whether we consider this light. In the
opposite limit that $\widetilde{q}\approx 1/2$, the stars in the CB are very
similar such that $\widetilde{q}\approx\widetilde{L}$, making $f_{\rm{pb}}$
very small and thus rendering irrelevant small differences in $f_{M}$.

####  3.2.2 Mass ratio distribution

The impact of an undetected CB on the WB observables depends on the fraction
of the CB mass in the undetected star (Figure 8). This requires us to assume
some distribution for $\widetilde{q}\leq 1/2$, or equivalently for the mass
ratio $q\equiv\widetilde{q}/\left(1-\widetilde{q}\right)\leq 1$. Following
Korntreff, Kaczmarek & Pfalzner (2012), we assume that

$\displaystyle P\left(q\right)\,\propto\,q^{0.4}\,,\quad q\leq 1\,.$ (21)

We expect this to broadly match the distribution of WB mass ratios. We
therefore use Figure 9 to show the mass fraction in the less massive star,
which is the quantity more relevant for our analysis (Equation 24). The
$q^{0.4}$ assumption is shown as a solid black line, while the distribution
for our WB sample is shown as a solid blue line. It is apparent that our
adopted $q^{0.4}$ distribution provides a good match to the WB distribution.
The WB sample contains a deficit of systems with a very extreme mass ratio,
but this is almost certainly due to selection effects: low mass stars are very
faint (Figure 1). Even if they are not present in our WB sample, we still need
to consider extreme mass ratio CBs in our analysis because they can affect the
WB observables if the CB has a small separation.

Figure 9: The mass ratio distribution of the WBs in our _Gaia_ DR3 sample
(thick solid blue line). We consider three different forms for the undetected
CB population. The nominal assumption of a $q^{0.4}$ dependency (solid black
line) is based on Korntreff et al. (2012), with $q\leq 1$ being the mass ratio
between the stars in each binary. This model fits the WB distribution fairly
well, especially given that selection effects would make it hard to find
binaries with $q\ll 1$. We also consider the case of a flat distribution in
$q$ (dot-dashed red line) and a linear distribution (dotted magenta line).
These should bracket the range of possible $q$ distributions.

To bracket the possible uncertainties in the CB mass ratio distribution, we
also consider a flat distribution of $q$ and a linear distribution, which we
show on Figure 9 as the dot-dashed red and dotted magenta line, respectively.
A flat distribution is a rough approximation to the three-part power law
recommended in figure 7 of El-Badry et al. (2019). We will see later that
using instead a flat or linear prior on $q$ gives similar results to our our
nominal $q^{0.4}$ assumption, so our conclusions are similar for any power-law
distribution of $q$ with exponent in the range $0-1$ (Section 5.1.2).

####  3.2.3 Semi-major axis distribution

The Keplerian orbital velocity of the CB is affected by its semi-major axis
$a_{\rm{int}}$. Following PS23, we impose an upper limit of $ak_{\rm{CB}}$,
where $a$ is the semi-major axis of the WB. Thus, we require that

$\displaystyle\widetilde{a}_{\rm{int}}~{}\equiv~{}\frac{a_{\rm{int}}}{a}<k_{\rm{CB}}\,.$
(22)

We work with $\widetilde{a}_{\rm{int}}$ rather than $a_{\rm{int}}$ because the
ratio between the CB and WB orbital velocities depends on the ratio of their
semi-major axes. The maximum allowed value of this ratio is $k_{\rm{CB}}$,
which is a free parameter in our analysis.

Since we remove systems with an obvious third star, $k_{\rm{CB}}$ refers to
the maximum $\widetilde{a}_{\rm{int}}$ for CB companions to the WBs in our
sample, bearing in mind that we do not consider stars which have a discernible
companion within $0.5\arcsec$ with a reasonably consistent parallax and proper
motion (Section 2.1). To convert this to a separation in AU, we note that as
our sample goes out to 250 pc, we can expect a reasonable number of systems
out to 200 pc. Thus, we expect WBs in our sample to typically not have CB
companions with projected separation beyond 100 AU. The upper limit to the
actual separation would be larger thanks to projection effects, but their
impact on $a_{\rm{int}}$ would be counteracted by the high likelihood of the
orbit being quite eccentric (Hwang et al., 2022) and the fact that orbital
phases near apocentre are more frequent. It is possible for a sufficiently
faint star to go undetected at an even larger separation, but it would have to
evade our faint star catalogue used to reduce CB contamination (Section 2.1).
This catalogue goes down to an apparent magnitude of $m_{G}=20$, which at a
distance of 200 pc corresponds to $M_{G}=20-5\log_{10}20=13.5$ and thus a mass
of about $0.1\,M_{\odot}$ (Figure 1). As a result, CB companions beyond
$0.5\arcsec$ from the contaminated star can remain undetected only in a narrow
range of mass at the very bottom of the main sequence. Such a low mass would
also reduce the perturbation to the velocity of the contaminated star (low
$f_{\rm{pb}}$ in Equation 20). Combined with the large separation, this makes
it unlikely that a distant low-mass companion which evades our faint star
catalogue would actually be relevant to our analysis. We thus assume that the
maximum relevant $a_{\rm{int}}\approx 100$ AU. To convert this to
$k_{\rm{CB}}$, we note that a reasonable $a$ for the WBs might be 40 kAU
because our WB sample goes up to $r_{\rm{sky}}=30$ kAU and we need to consider
all viewing angles and orbital phases. For this reason, our WB orbit library
goes up to $a=57$ kAU, though we expect our analysis to assign only very small
statistical weights at such high $a$ due to its declining distribution
(Equation 16) and the fact that we truncate orbits when they go beyond 100 kAU
to mimic disruption caused by Galactic tides and encounters with other stars.
Taking the ratio of the above estimates for $a_{\rm{int}}$ and $a$ suggests
that $k_{\rm{CB}}\approx 2.5$‰. We do not impose this in our analysis, but we
do set a prior that it should be at least 1‰ to avoid numerical difficulties.
This is very conservative because it implies a maximum $a_{\rm{int}}$ of about
40 AU, which corresponds to just $0.2\arcsec$ for a WB 200 pc away. One can
easily envisage that a CB with a larger separation would not be detected at
this distance.

At the opposite end, the lowest $a_{\rm{int}}$ that we need to consider is set
by the fact that CBs on very small orbits have a short orbital period, leading
to rapid but low amplitude astrometric oscillations that have little impact on
the inferred space velocity of the contaminated star. This phenomenon is not
directly included in our analysis, which for simplicity assumes that the stars
in a CB move at constant velocity over the whole _Gaia_ observing baseline. In
Appendix B, we estimate that this approximation breaks down when
$a_{\rm{int}}\la 3.2$ AU, below which the impact of the CB rapidly becomes
much smaller than our estimate based on uniform linear motion.

Based on the above considerations, we allow CBs to have a 1.5 dex range of
$\widetilde{a}_{\rm{int}}$ with a maximum at $k_{\rm{CB}}>1$‰. We assume that
$\widetilde{a}_{\rm{int}}$ has a flat distribution in log-space, as in the
classical Öpik law (Öpik, 1924). The logarithmic distribution of
$a_{\rm{int}}$ does indeed seem to be rather flat between 10 and 100 AU
(figure 2 of Offner et al., 2023), so this ought to be a reasonable
assumption. We will see later that our results are not much affected by using
a somewhat different shape or width to the $a_{\rm{int}}$ distribution
(Section 5.1.3).

In addition to properties of the CB, how it affects $\widetilde{\bm{v}}$ of
the WB also depends on the latter’s projected separation. At fixed $a$, the
impact on $\widetilde{\bm{v}}$ is reduced at low $r_{\rm{sky}}$ because this
raises the predicted Newtonian $v_{c}$ of the WB (Equation 3). Considering
also that the CB orbital velocity $\propto 1/\sqrt{\widetilde{a}_{\rm{int}}}$,
the perturbation to the WB $\widetilde{\bm{v}}$ scales as

$\displaystyle\Delta\widetilde{\bm{v}}~{}\propto\sqrt{\frac{\widetilde{r}_{\rm{sky}}}{\widetilde{a}_{\rm{int}}}}\,.$
(23)

The lack of any explicit dependence on $a$ greatly reduces the computational
cost of the WBT, but this is only possible because the WB libraries are stored
with respect to $\widetilde{r}_{\rm{sky}}$ rather than $r_{\rm{sky}}$. The
only slight downside is that we need the latter for a comparison with
observations, so we subsequently need to find $r_{\rm{sky}}\equiv
a\widetilde{r}_{\rm{sky}}$.

####  3.2.4 Convolving the CB and WB libraries

The fractional increase to the total mass of a WB is found by combining
$\Delta\widetilde{M}$ for each of the stars in it, which we denote using A and
B subscripts.

$\displaystyle\Delta\widetilde{M}~{}=~{}\widetilde{q}_{A}\Delta\widetilde{M}_{A}\,+\,\widetilde{q}_{B}\Delta\widetilde{M}_{B}\,,$
(24)

where $\widetilde{q}_{A}\equiv 1-\widetilde{q}_{B}$ is the fraction of the WB
total mass in star A and $\Delta\widetilde{M}_{A}$ is the fraction by which
its actual mass exceeds its apparent mass (Equation 18). If one of the WB
stars is not in a CB, then $\Delta\widetilde{M}=0$ for that star. We assume
that $\widetilde{q}_{A}$ and $\widetilde{q}_{B}$ have an identical
distribution, which we obtain from our WB sample (solid blue line in Figure
9). This is not exactly correct because undetected CBs can alter the mass
ratio. However, we expect this to be only a very minor issue because even in
the extreme case of having two equal Solar mass stars in a WB where one of the
stars has an undetected equal mass companion, observers would infer a mass
fraction in the less luminous component of 0.46 rather than the correct value
of 1/3. The error would be smaller in more typical cases, especially as both
stars in the WB could actually be CBs. We therefore expect the WB mass ratio
inferred from luminosities to be quite accurate.

To obtain the relative velocity between the photocentres and barycentres of
the CBs making up a WB, we again need to combine the contributions from both
CBs. In this case, the scaling factors are $\sqrt{\widetilde{q}_{A}}$ and
$\sqrt{\widetilde{q}_{B}}$ because the Keplerian orbital velocity of a CB
scales with the square root of its mass. Allowing also for the increased WB
orbital velocity due to its mass being higher than estimated, we get that the
WB $\widetilde{\bm{v}}$ should be revised as follows:

$\displaystyle\widetilde{\bm{v}}~{}\to~{}\widetilde{\bm{v}}f_{M}\,+\,\Delta\widetilde{\bm{v}}_{A}\sqrt{\widetilde{q}_{A}}\,+\,\Delta\widetilde{\bm{v}}_{B}\sqrt{\widetilde{q}_{B}}\,,$
(25)

where $f_{M}$ was defined in Equation 20 and $\Delta\widetilde{\bm{v}}_{A}$ is
the change to $\widetilde{\bm{v}}$ of the WB if all its mass was in star A
such that $\widetilde{q}_{A}=1$ ($\Delta\widetilde{\bm{v}}_{B}$ is defined
analogously). Since $\Delta\widetilde{\bm{v}}$ induced by each CB is a 2D
vector, we also need to allow for partial cancellation between the
photocentre-barycentre offsets arising from each CB. This leads to 4 possible
outcomes for whether the contributions along each direction are parallel or
opposite.

Following El-Badry et al. (2019), we implement a $\delta$-function in the
probability distribution of the CB mass ratio. Such an equal mass ratio
population is also evident in our colour-magnitude diagram (dashed grey lines
in Figure 2), with our analysis suggesting that the equal mass likelihood
$P_{\rm{eqm}}=0.04$ (see Appendix A). The origin of such a twin population is
unclear, but it could be related to gas accretion onto a forming binary star
tending to equalize the masses (Tokovinin, 2017; Tokovinin & Moe, 2020;
Tokovinin, 2023). To include this, we combine the CB perturbations calculated
assuming a smooth distribution of the CB mass ratio with similar calculations
assuming an exactly equal mass ratio for the CB, weighting the former at
$\left(1-P_{\rm{eqm}}\right)$ and the latter at $P_{\rm{eqm}}$ (further
details are provided in Banik et al., 2021). We demonstrate later that
changing the assumed CB mass ratio distribution has very little impact on our
results (Section 5.1.2).

The overall impact of CBs on the WB population is modulated by the likelihood
$f_{\rm{CB}}$ that a star in our WB sample has an undetected CB companion. We
assume that CB contamination of the two stars in each WB occurs independently,
so the likelihood of both stars being contaminated is $f_{\rm{CB}}^{2}$, which
is not negligible. We found that one of the most computationally expensive
parts of the WBT is calculating the CB perturbations to the WB parameters in
this double contaminated case. This is due to the need to consider a range of
CB semi-major axes, velocity perturbations, and mass contributions from both
CBs, which almost doubles the already large number of nested loops.

Table 1: The $r_{\rm{sky}}$ and $\widetilde{v}$ range used for the WBT and the binning scheme, which we fixed in advance to limit the possibility of bias (Banik et al., 2021). $r_{\rm{sky}}$ (kAU) | $\widetilde{v}$
---|---
$2-5$ (6 bins) |  $0-1.6$ (20 bins)
$5-10$ (5 bins) |  $1.6-2.6$ (5 bins)
$10-16$ (3 bins) |  $2.6-5$ (5 bins)
$16-22$ (2 bins) |
$22-30$ (2 bins) |

Combining the CB perturbations with the WB library is non-trivial because CBs
cannot be considered as adding an extra population to the statistics $-$ they
physically affect the WB. As summarized in Equation 25, we need to scale up
the WB $\widetilde{\bm{v}}$ to account for the fractional mass increase
$\Delta\widetilde{M}$ induced by the CB(s) and then add the velocity
perturbations from each CB in the directions parallel and orthogonal to
$\widetilde{\bm{v}}$ of the WB, which is a 2D vector in the sky plane. We then
use the Pythagoras Rule to calculate the resultant $\widetilde{v}$, or more
precisely the appropriate $\widetilde{v}$ pixel (Table 1). Since there is a
reasonable chance that the resultant $\widetilde{v}>5$, CBs can reduce the
overall normalization of the $\left(r_{\rm{sky}},\widetilde{v}\right)$
distribution in our WB library. This may well have happened in reality, but
since our _Gaia_ DR3 sample only contains WBs with $\widetilde{v}<5$, we
renormalize the $\left(r_{\rm{sky}},\widetilde{v}\right)$ distribution that we
obtain after the CB-WB convolution step. The computational cost of convolving
the CB and WB libraries is kept low using the techniques detailed in Banik et
al. (2021).

###  3.3 Chance alignments

So far we have only considered bound systems. It is inevitable that some of
the WBs in our sample will be unbound. Observations are unlikely to catch a WB
in the process of disruption due to a passing third star or a molecular cloud
(see section 8.1 of BZ18, ). Even so, chance alignments of stars can still
arise, especially when we get to large $r_{\rm{sky}}$ and $\widetilde{v}$ (see
figure 7 of PS19, ). The likelihood is enhanced due to the possibility of the
stars being born in the same star cluster which later dissolved (Oh et al.,
2017; Dinnbier et al., 2022), giving them a much lower relative velocity than
the Galactic velocity dispersion. This raises the chance of observing the
stars close to each other on the sky, which could cause them to be
misidentified as a WB given the _Gaia_ parallax uncertainty. LOS contamination
should be reduced somewhat by excluding regions of the sky with an enhanced
density of stars, e.g. known star clusters (table 1 of PS23, ).

To limit the complexity of our analysis, we set up a very simple model for LOS
contamination under the assumption that the separation and relative velocity
of stars born in the same cluster can greatly exceed that of the WBs in our
sample. Since $\bm{r}_{\rm{sky}}$ and $\bm{v}_{\rm{sky}}$ are both 2D vectors
within the sky plane, we assume that the number density of chance-aligned WBs
scales as

$\displaystyle\frac{dN_{\rm{LOS}}}{dr_{\rm{sky}}\,dv_{\rm{sky}}}~{}\propto~{}r_{\rm{sky}}v_{\rm{sky}}\,.$
(26)

Bearing in mind that $v_{\rm{sky}}\propto\widetilde{v}/\sqrt{r_{\rm{sky}}}$
(Equation 3), the population distribution of LOS contamination is

$\displaystyle\frac{dN_{\rm{LOS}}}{dr_{\rm{sky}}\,d\widetilde{v}}~{}\propto~{}\widetilde{v}\,.$
(27)

The dependence on $r_{\rm{sky}}$ cancels out because raising $r_{\rm{sky}}$
reduces the range of $\bm{v}_{\rm{sky}}$ for the same range in
$\widetilde{v}$.

We integrate this distribution across each of the $r_{\rm{sky}}$ and
$\widetilde{v}$ pixels shown in Table 1 to obtain the relative number of
chance-aligned systems that we expect to find in each pixel. We then normalize
the results to a sum of $f_{\rm{LOS}}$ and combine with the result of our CB-
WB convolution (Section 3.2.4) weighted by $\left(1-f_{\rm{LOS}}\right)$,
where $f_{\rm{LOS}}$ is the LOS contamination fraction. We leave
$f_{\rm{LOS}}$ as a free parameter in our analysis. Since there are only a few
WBs in our catalogue at high $r_{\rm{sky}}$ and $\widetilde{v}$ where LOS
contamination would be most significant (Figure 6), we expect $f_{\rm{LOS}}$
to be at most a few percent.

###  3.4 The simulated distribution of WBs

Figure 10: How we build up the simulated $\widetilde{v}$ distribution in our
best MOND model (Section 4.1), with results shown in four different
$r_{\rm{sky}}$ ranges (different panels). In each case, the distribution due
to WBs alone (blue) is significantly altered by including CBs (red). Further
adding LOS contamination yields our final simulated distribution (black),
which is used to normalize the results shown here. This helps to highlight
that including LOS contamination slightly raises the predicted number of WBs
at high $r_{\rm{sky}}$.

To simulate the distribution of WBs, we tie together a model for the WB
population (Section 3.1) with a model for an undetected population of CB
companions (Section 3.2) and then allow for chance alignments (Section 3.3).
We use Figure 10 to illustrate the separate stages for the particular example
of our best-fitting MOND model (Section 4.1). The simulated $\widetilde{v}$
distribution is shown in four different $r_{\rm{sky}}$ ranges considering only
WBs (blue curves), including CBs (red curves), and finally including also the
LOS contamination (black curves). Since only a small proportion of the WB
sample is expected to lie in any single $r_{\rm{sky}}$ interval used here, we
normalize the $\widetilde{v}$ distributions using the final simulated
distribution, so only the black curve in each panel is guaranteed to have an
integral of 1. WBs alone almost never yield $\widetilde{v}\ga 1.5$, even in
MOND $-$ the theoretical limit of 1.7 only refers to systems on highly
eccentric orbits seen close to pericentre from a direction nearly orthogonal
to the relative velocity (this is also apparent in figure 3 of BZ18, ). CB
companions create an extended declining tail that is easily able to reach much
higher $\widetilde{v}$. We capture much of this extended tail by going up to
$\widetilde{v}=5$, a choice we discussed in Section 2.1 and fixed in advance
of the WBT to mitigate possible biases (Banik et al., 2021). This extra
leverage allows our analysis to constrain the CB population, which dominates
the behaviour once $\widetilde{v}\ga 2$. A flat log-space distribution of
$a_{\rm{int}}$ (Öpik, 1924) yields a $\widetilde{v}$ distribution that is also
flat in log-space in the CB-dominated region, helping to explain why the
distribution here is $\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}1/\widetilde{v}$. We expect this to flatten
out and ultimately start rising at very high $\widetilde{v}$ due to LOS
contamination. This is most evident in the highest $r_{\rm{sky}}$ interval,
where there is less contribution from only the gravitationally bound systems
because WBs have a declining distribution of $r_{\rm{sky}}$ (Equation 16). Our
WB sample reaches high enough $\widetilde{v}$ for its distribution to become
flat, but a clearly rising trend is not apparent because we have limited the
considered $\widetilde{v}$ range to avoid being swamped by chance alignments.
Other studies which consider a wider parameter range do show a rising trend
(e.g. see figure 12 of El-Badry et al., 2021).

###  3.5 Comparison with observations

We compare the observed number of WBs in each
$\left(r_{\rm{sky}},\widetilde{v}\right)$ pixel with the expectation of our
model for that pixel. We use binomial statistics because we are dealing with
integer statistics and the total sample size $N$ is finite. Thus, the binomial
likelihood of observing $k$ systems in a pixel is

$\displaystyle
P_{\rm{pixel}}~{}=~{}\frac{N!}{\left(N-k\right)!k!}p^{k}\left(1-p\right)^{N-k}\,,$
(28)

where $p$ is the fraction of the total number of WBs in the comparison region
($r_{\rm{sky}}=2-30$ kAU, $\widetilde{v}<5$) which should be located within
the pixel under consideration. For computational reasons, we use the
logarithmic version of this equation. The log-likelihood $\ln P$ of the model
is then given by considering all 540 pixels summarized in Table 1.

$\displaystyle\ln P~{}=~{}\sum_{\rm{Pixels}}\ln P_{\rm{pixel}}\,.$ (29)

This is used to infer the optimal values and confidence intervals for our
model parameters associated with the CBs ($f_{\rm{CB}}$, $k_{\rm{CB}}$, and
$\gamma$), LOS contamination ($f_{\rm{LOS}}$), and the WBs ($a_{\rm{break}}$,
$\beta$, and $\alpha_{\rm{grav}}$). Note that $\gamma$ affects both the CB and
WB populations because they are assumed to follow the same eccentricity
distribution for simplicity. The CBs are assumed to follow an Öpik law
distribution of semi-major axes (Öpik, 1924) over a 1.5 dex range (Section
3.2.3), while the WBs follow a broken power law that reduces to $a^{-1.6}$ at
small separations (Equation 16).

Notice that regardless of how many WBs a pixel contains, $\ln
P_{\rm{pixel}}<0$ in every pixel as its predicted $p$ is always non-zero once
we include CB and LOS contamination. Since we have quite a large number of WBs
spread across a large number of pixels, we expect $\ln P\ll 0$ (a similar
situation is evident in table 2 of Asencio et al., 2022). As a crude estimate,
we can suppose that our model predicts $8611/540=16$ WBs in each of our 540
pixels, so the actual number of WBs in any pixel is expected to nearly follow
a Poisson distribution with a standard deviation of $\sqrt{16}=4$. If the
model works well, then the actual number of WBs would be approximately
uniformly distributed over the range $12-20$, with all of these outcomes
having a likelihood of $P_{\rm{pixel}}=1/9$. The combined likelihood would
then be $\ln P=-540\,\ln 9=-1200$. We therefore expect to get a combined log-
likelihood of about this much, but this is of course only a very crude
estimate. In reality, not all pixels would be expected to have the same number
of WBs, and even if they were, outcomes in the tail of the Poisson
distribution are less likely than the inverse of the Poisson noise. Moreover,
it is inevitable that our model does not capture all the subtleties of the
actual WB population because of the simplifying assumptions made to keep the
complexity and computational cost manageable. Nonetheless, we will see later
that the above estimate is fairly accurate.

##  4 Results

Before presenting the results of our detailed statistical analysis discussed
in the previous section, we first look for trends in the median
$\widetilde{v}$ of our WB population with $r_{\rm{sky}}/r_{{}_{M}}$, which is
a proxy for the internal acceleration of each WB. For this purpose, we sort
our WBs in order of $r_{\rm{sky}}/r_{{}_{M}}$ and use this to create ten
equally sized subsamples with no overlap. The thin magenta line in Figure 11
shows the median $\widetilde{v}$ as a function of the median
$r_{\rm{sky}}/r_{{}_{M}}$ of each subsample. A clear rising trend is evident,
with the increase being close to the 20% expected in MOND. However, LOS
contamination would be more important at high $r_{\rm{sky}}$ and high
$\widetilde{v}$, potentially inflating the median $\widetilde{v}$ in a manner
which appears acceleration-dependent. To mitigate against this possibility, we
note that MOND would cause an acceleration-dependent broadening to the main
peak of the $\widetilde{v}$ distribution at $\widetilde{v}\la 2$ (Figure 6;
see also figure 7 of PS19, ). In MOND, we expect almost no WBs with
$\widetilde{v}>1.5$ in the absence of contaminating effects (see figure 3 of
BZ18, ). We therefore consider the median $\widetilde{v}$ of only those WBs in
each $r_{\rm{sky}}/r_{{}_{M}}$ bin with $\widetilde{v}<1.5$, 2, or 2.5, which
we show in Figure 11 using the red, black, and blue line, respectively. In all
three cases, the rising trend disappears: the median $\widetilde{v}$ becomes
flat with respect to our proxy for the internal acceleration. There is of
course some scatter, but we would expect a random scatter of about
$0.6/\sqrt{861}=0.02$ due to the finite number of WBs, which is roughly in
line with the actual scatter.

Figure 11: The median $\widetilde{v}$ of our WB sample as a function of the
median $r_{\rm{sky}}/r_{{}_{M}}$, shown for ten bins in the latter which each
contain 861 WBs except the last bin which has 862. The solid lines show
results with four different upper limits to $\widetilde{v}$ of 1.5 (red), 2
(black), 2.5 (blue), and 5 (thin magenta). Qualitatively, the MOND expectation
is that the median $\widetilde{v}$ rises for larger $r_{\rm{sky}}/r_{{}_{M}}$
and then becomes flat once the Galactic gravity dominates. We indicate this by
two dashed grey lines with different normalizations. Their shape is based on
Figure 7, which implies that the main peak in the $\widetilde{v}$ distribution
for $\widetilde{v}\la 2$ should broaden at low accelerations. Neglecting
projection effects and assuming the current separation is representative of
the whole orbit, this would lead to a rise in the median $\widetilde{v}$ by a
factor of $\sqrt{\eta}=1.2$ that occurs roughly linearly over the range
$r_{\rm{sky}}/r_{{}_{M}}=0.2-1$ before rapidly flattening out. The Newtonian
prediction is a flat relation, which is consistent with the data up to
$\widetilde{v}<2.5$, within which lies 89.8% of the sample. The full sample of
WBs (with $\widetilde{v}<5$) does show a rising median that looks
qualitatively similar to the MOND expectation, but since this rise is not seen
for the subsamples restricted to the main peak region of the $\widetilde{v}$
distribution where $\widetilde{v}$ is low enough to be physically plausible
for an uncontaminated system, we attribute this behaviour to LOS contamination
becoming more important at larger separations (Section 3.3).

Since the upper limit to $\widetilde{v}$ should substantially alleviate the
impact of CB and LOS contamination without much affecting a genuine MOND
signal, it is clear that the apparent MOND signal discussed above is almost
certainly a consequence of contaminating effects, highlighting how easily one
can reach erroneous conclusions from the WBT using simple population
statistics without a good understanding of the astrophysical systematics
(McCulloch & Lucio, 2019; Hernandez et al., 2019, 2022). Focusing on WBs with
$\widetilde{v}\la 2$, the very weak dependence of the median $\widetilde{v}$
on $r_{\rm{sky}}/r_{{}_{M}}$ despite the significant range covered by our
sample strongly suggests that the WBT will return a Newtonian result. We
highlight this using the dashed grey lines in Figure 11, which approximately
show the MOND prediction that the median $\widetilde{v}$ should rise by a
factor of $\sqrt{\eta}=1.2$ over the range $r_{\rm{sky}}/r_{{}_{M}}=0.2-1$.
This is motivated by Figure 7, which shows that the enhancement factor to the
radial Newtonian gravity rises from 1 up to its asymptotic value roughly
linearly over this range of $r/r_{{}_{M}}$ and is almost flat afterwards.
While this is a rather crude way of considering the situation because it
neglects projection and orbital phase effects, it should approximately capture
the expected broadening of the $\widetilde{v}$ distribution in the main peak
region as one considers WB subsamples with a lower internal acceleration. The
expected Newtonian behaviour is not shown as there is already a gridline at
0.57. This provides a much better fit to the median $\widetilde{v}$ for the
$\widetilde{v}<2$ subsample than the approximate MOND expectation mentioned
above.

This strong hint that the WBT will disfavour MOND is in line with the detailed
analysis of PS23, who prepared a model for the WB population similar to ours
and found a strong preference for Newtonian gravity over MOND. Their Newtonian
model fits the observations surprisingly well given their limited exploration
of the parameter space (see their figure 12). Moreover, the right panel of
their figure 17 shows that WBs have a median $v_{\rm{sky}}\propto
1/\sqrt{r_{\rm{sky}}}$, especially when restricting to only those WBs with
$\widetilde{v}\la 2$. This is the expected behaviour if Kepler’s Third Law
works even in the low-acceleration regime, as illustrated by the dashed lines
on this figure.

The lack of a trend in the median $\widetilde{v}$ is unlikely to be caused by
observational uncertainties, which would generally lead to a larger
$\widetilde{v}$ uncertainty at high $r_{\rm{sky}}$ due to the lower Newtonian
$v_{c}$ (Equation 3). This would if anything create the appearance of a MOND-
like signal in actually Newtonian data rather than precisely cancel out a
genuine MOND signal (c.f. Section 5.2). Moreover, the $\widetilde{v}$
uncertainties in our WB sample are very small (Figure 4).

Although our results favour the Newtonian model, the distribution of
$\widetilde{v}$ contains much more information than just its median value.
Before drawing strong conclusions about the gravity law, we try to exploit
this information using a forward model that considers the most relevant
factors (Section 3). We now turn to the results of a comparison between this
model and the WB dataset.

###  4.1 Gradient ascent with a fixed gravity law

To check if Newtonian gravity indeed provides a better fit to the data than
MOND, we fix the gravity law to Newtonian ($\alpha_{\rm{grav}}=0$) and run a
gradient ascent in the remaining model parameters to try and maximize $\ln P$