Robust Evidence for the Breakdown of Standard Gravity at Low Acceleration from Statistically Pure Binaries Free of Hidden Companions

It is found that Gaia DR3 binary stars selected with stringent requirements on astrometric measurements and radial velocities naturally satisfy Newtonian dynamics without hidden close companions when projected separation $s \lesssim 2$ kau, showing that pure binaries can be selected. It is then found that pure binaries selected with the same criteria show a systematic deviation from the Newtonian expectation when $s \gtrsim 2$ kau. When both proper motions and parallaxes are required to have precision better than 0.005 and radial velocities better than 0.2, I obtain 2,463 statistically pure binaries within a `clean' $G$-band absolute magnitude range. From this sample, I obtain an observed to Newtonian predicted kinematic acceleration ratio of $\gamma_g=g_{\rm{obs}}/g_{\rm{pred}}=1.49^{+0.21}_{-0.19}$ for acceleration $\lesssim 10^{-10}$ m s$^{-2}$, in excellent agreement with $1.49\pm 0.07$ for a much larger general sample with the amount of hidden close companions self-calibrated. I also investigate the radial profile of stacked sky-projected relative velocities without a deprojection to the 3D space. The observed profile matches the Newtonian predicted profile for $s \lesssim 2$ kau without any free parameters but shows a clear deviation at a larger separation with a significance of $\approx 5.0\sigma$. The projected velocity boost factor for $s\gtrsim 5$ kau is measured to be $\gamma_{v_p} = 1.20\pm 0.06$ (stat) $\pm 0.05$ (sys) matching $\sqrt{\gamma_g}$. Finally, for a small sample of 40 binaries with exceptionally precise radial velocities (fractional error $<0.005$) the directly measured relative velocities in the 3D space also show a boost at larger separations. These results robustly confirm the recently reported gravitational anomaly at low acceleration for a general sample.


INTRODUCTION
Wide binaries (widely-separated, long-period, gravitationally-bound binary stars) provide crucial testbeds for probing gravitational dynamics in the lowacceleration regime (e.g., Hernandez et al. 2012;Banik & Zhao 2018;Pittordis & Sutherland 2018;Banik & Kroupa 2019;Pittordis & Sutherland 2019;Hernandez et al. 2022).A couple of recent studies by Chae (2023) and Hernandez (2023) of wide binary stars selected from Gaia data release 3 (DR3; Vallenari et al. 2023) have reported a gravitational anomaly at low acceleration ≲ 10 −9 m s −2 , or for a sky-projected separation s ≳ 2 kau (kilo astronomical unit) for typical binaries with total masses of ∼ (0.5 − 2)M ⊙ .This gravitational Corresponding author: Kyu-Hyun Chae chae@sejong.ac.kr, kyuhyunchae@gmail.comanomaly implies a low-acceleration breakdown of both Newtonian dynamics and general relativity and so has immense implications for astrophysics, cosmology, and fundamental physics.Thus, one cannot overemphasize the importance of confirming the claimed anomaly from as many independent studies as possible.
Chae (2023) considered wide binaries selected from El-Badry et al. ( 2021) that are statistically free of both chance-alignment cases and resolved (> 1 ′′ ) triples and higher-order multiples.Because of the initial selection, additional quality cuts, and the availability of dust extinction correction, Chae (2023) used only up to 26, 615 wide binaries within 200 pc.Chae (2023) then selfcalibrated the occurrence rate (f multi ) of triples and higher-order multiples hosting hidden (i.e.unresolved) close companions by requiring that binaries must satisfy Newtonian dynamics at a close enough separation, or at a high enough acceleration ≳ 10 −8 m s −2 , as predicted by all currently available plausible theories including modified Newtonian dynamics (MOND; Milgrom 1983).Chae (2023) also paid a particular attention to projection effects and employed a Monte Carlo (MC) method to deproject measured sky-projected relative velocities v p to the three-dimensional (3D) space physical velocities v, and compared a kinematic acceleration v 2 /r with the corresponding Newtonian prediction.Chae (2023) obtained up to a 10σ significance for the gravitational anomaly based on MC analyses.Moreover, the magnitude and trend of the anomaly matched well the prediction of MOND-type modified gravity such as AQUAL (Bekenstein & Milgrom 1984) and QUMOND (Milgrom 2010) under the external field effect (EFE) of the Milky Way.Hernandez (2023) took a different approach that tried to remove all cases of both chance alignments and hierarchical multiples.Because Hernandez (2023) applied various strict cuts, his final sample includes only 450 pure binaries in the distance range d < 125 pc or 125 < d < 170 pc.Hernandez (2023) calculated the dispersion of one-dimensional velocity components on the plane of the sky and compared it with the Newtonian prediction by Jiang & Tremaine (2010).Hernandez (2023) checked that small-separation (s ≲ 2 kau) systems matched the Newtonian prediction indicating that kinematic contaminants are indeed negligible and Newtonian dynamics holds in the high-acceleration regime.Then, he found that the observed sky-projected velocities systematically deviated from the Newtonian expectation at large-separation (s ≳ 2 kau) systems.Because the sample size was small and the observed kinematics of the binaries was compared with simulations for other binaries, Hernandez (2023) did not quantify a statistical significance of the anomaly seen at s ≳ 2 kau.Nevertheless, the final result from Hernandez (2023) agreed with that of Chae (2023).
Unlike Chae (2023) and Hernandez (2023), another recent study by Pittordis & Sutherland (2023) based on a Gaia database considered only low-acceleration binaries and thus could not calibrate f multi among their wide binaries.Their analysis was also compounded by their inclusion of chance-alignment cases.Pittordis & Sutherland (2023) erroneously concluded that Newtonian dynamics matched their low-acceleration data with their "fitted" f multi (without a proper calibration).As Chae (2023) demonstrated, Pittordis & Sutherland (2023) conclusion is not surprising because the dynamics uncovered by Chae (2023) is pseudo-Newtonian with a rescaling of Newton's constant G → 1.4G and this boost can be canceled by a higher f multi (without a proper calibration).
Just recently, Banik et al. (2023) argued that Newtonian dynamics was preferred over MOND (using specifically the QUMOND model) based on a method similar to Pittordis & Sutherland (2023).They did not include the Newtonian regime (≲ 2 kau) data to calibrate f multi but claimed that gravity and f multi could be simultaneously constrained based only on data from the Newton-MOND transition and MOND regimes (see below).They obtained a high value of f multi ≈ 0.70 for their preferred Newtonian model although their sample included only binaries passing a strict cut on relative velocities.Their value of f multi ≈ 0.70 is unlikely high compared with the observed range 0.3 ≲ f multi ≲ 0.5 (e.g., Raghavan et al. 2010;Riddle et al. 2015;Moe & Stefano 2017) even for general binary samples without kinematic cuts.Moreover, I note that binaries selected by strict kinematic cuts have significantly lower f multi than that for a general sample.I will further discuss their sample selection, analyses, and results at relevant places (see in particular appendices).
Figure 1 shows an AQUAL numerical (Chae & Milgrom 2022) prediction on the radial acceleration for circular orbits under the EFE of the Milky Way that approximately matches the AQUAL analytic asymptotic limit at an average inclination of the external field.Due to the strong external field of the Milky Way, internal dynamics is expected to switch from the Newtonian regime to the pseudo-Newtonian regime with a boosted Newton's constant.As shown in Figure 1, most of the transition is expected to occur abruptly in the narrow acceleration range of −9.6 ≲ log 10 (g N /m s −2 ) ≲ −8.8.As shown in Figure 2, for typical Gaia wide binaries with a total mass of ≈ 1.4M ⊙ the transition acceleration range corresponds approximately to the sky-projected separation range of 2 ≲ s ≲ 5 kau (kilo astronomical units).This MOND prediction was supported by the two published analyses (Chae 2023;Hernandez 2023) and is intended to be further tested in this study.
Here I consider a new analysis that is complementary to Chae (2023) and Hernandez (2023)) and can provide a robust test of gravity.The analysis of Chae (2023) involved a complex chain of steps with various observational inputs for a general sample and obtained a maximal statistical power.The complexity arises largely from modeling the kinematic effects of hidden close companions with currently available observational inputs.As f multi → 0, the complexity is gone and thus any possibilities of systematic errors involving close companions can be removed at the cost of losing statistical power.The question is whether a statistically significant sample of f multi = 0 can be selected in a systematic and veri- fiable way so that gravitational anomaly can be tested robustsly and with a sufficient statistical power.Hernandez (2023) obtained a sample that was supposed to be largely free of hierarchical systems but his analysis was limited in two ways.Hernandez (2023) did not carry out a Monte Carlo simulation or any other statistical procedure for his sample to quantify a statistical significance of the seen anomaly (although a comparison with an independent simulation by Jiang & Tremaine (2010) was made).Also, each sample defined by Hernandez (2023) seems too small with just 450 binaries.
Here I obtain a much larger sample of 2,463 wide binaries of f multi = 0 in a systematic and verifiable way.Then, I test gravity in an acceleration plane with the algorithm developed in Chae (2023) with f multi = 0.More importantly, I investigate stacked velocity profiles with a MC simulation to do a quantitative statistical test.In this way, the present work will complement both Chae (2023) and Hernandez (2023).I will show that the results from pure binaries agree excellently with those of Chae (2023) reaffirming the validity of the procedure and conclusion of Chae (2023).
The structure of this paper is as follows.Section 2 describes how a sample of pure binaries can be selected in a systematic way.Section 3 describes a Monte Carlo modeling of pure binary stars.Section 4 presents the results.Section 5 discusses any possible source of systematic errors.In Section 6, I offer the conclusion and discuss future works.In Appendix A, I describe a correction to Chae (2023) and revise the representative results.In Appendix B, I consider some kinematic quality cuts on general samples of wide binaries and their effects on inference on gravity.Python scripts used for this work and the sample of pure binaries can be accessed at Zenodo: doi:10.5281/zenodo.8416435.

A SYSTEMATIC SELECTION OF PURE BINARIES
Following Chae (2023), I work with the catalog of one million candidate binaries derived by El-Badry et al. (2021) from Gaia DR3 astrometric measurements.This catalog provides estimated values of chance-alignment probability (R) so that chance-alignment cases can be effectively excluded.The catalog also excludes triples and higher-order multiples whose components are all resolved by > 1 ′′ .Thus, by requiring R < 0.01 (or something similar), one can choose binaries that may include only unresolved close companions.
Specifically, I consider 26,615 binaries with R < 0.01 within 200 pc from the Sun whose components have sky-projected separation (s) in the range 0.2 < s < 30 kau and have absolute magnitudes in the 'clean range' 4 < M G < 14 defined by Chae (2023) where M G is the dust-extinction-corrected absolute magnitude in the Gaia G band.As Chae (2023) showed, samples of binaries selected by a precision threshold of 0.01 or 0.005 imposed on the measured proper motions (PMs) require f multi ≳ 0.3 to statistically satisfy Newtonian dynamics at high enough accelerations ≳ 10 −8 m s −2 .
Because binaries must satisfy Newtonian dynamics at high enough accelerations, binaries selected in a systematic way may be regarded as a sample of pure binaries statistically free of hidden close companions if they require f multi = 0 to match Newtonian dynamics at accelerations ≳ 10 −8 m s −2 .For this procedure to be valid masses of individual stars must be reliably known.Fortunately, Chae (2023) provides a couple of reliable magnitude-mass relations (see figure 7 and table 1 of Chae (2023)) in the Gaia G band.In search of a sample of pure binaries satisfying f multi = 0, various selection criteria have been tried guided by observational and simulation studies (e.g., Belokurov et al. 2020;Penoyre et al. 2022).
If a close undetected companion is present, it can have various effects.First, the image may not be well modeled so that Gaia derived ruwe values may be significantly larger than 1.Second, because the close companion induces additional motions, the measured uncertainties of parallaxes and proper motions (PMs) will become larger.Third, the additional motions will also increase the measurement uncertanties of radial velocities.
It turns out that a sample of pure binaries can be obtained by requirements on both astrometric measurements and radial velocities.The requirements can be specified as follows.Throughout the brighter (more massive) star is referred to as component A while the fainter star component B.
1. Reported values of ruwe for both components are smaller than 1.2.
3. Distances of two components agree within 2σ uncertainties and a maximum possible difference from the elliptical orbit ∆d max orbit , i.e., where ∆d max orbit = 6s from a 99% statistical limit from random inclinations and orbital phases for elliptical orbits of observational eccentricities.
4. Radial velocities of both components have relative measurement errors smaller than 0.2.In other words, we take only binaries whose components have radial velocities with S/N > 5.
5. Radial velocities of two components agree within 2σ uncertainties and a maximum possible difference from the elliptical orbit ∆v max r,orbit , i.e., where M tot is the total mass of the binary system in units of solar mass (M ⊙ ) and s is the sky-projected separation between the two components.The factor 1.3 represents a maximum possible value (see Section 3 below) arising from random inclinations and orbital phases for elliptical orbits of observational eccentricities.The last factor 1.2 allows for a possible boost of velocity in MOND-type modified gravity theories so as not to preclude such theories, though it is practically not important because other uncertainties are larger.
In the above, the third and fifth requirements gaurantee that two stars form a true binary system.When ε = 0.005 is used, the total number of pure binary systems is N tot = 2463.The distribution of masses of the selected binaries can be found in Figure 3.The total masses are in the range 0.5 ≲ M tot /M ⊙ ≲ 2.1 with a mass ratio M B /M A ≥ 0.5 for 88% of the systems.The mean (median) mass for the entire sample is 1.35 (1.36) M ⊙ , and the binned mean varies only mildly with s as the right panel of Figure 3 shows.I also consider a subsample within a narrow range of total mass 1.1 < M tot /M ⊙ < 1.8.This subsample has a mean (median) mass of 1.44 (1.43) M ⊙ and its binned mean does not vary with s.
Since the mean mass varies little or not at all with s in the entire sample or the subsample, it is possible to investigate a stacked radial velocity profile.This is important because a statistical analysis of velocity profiles will be a main part of this work.
For the observed right ascension (α) and declination (δ) components of the PMs in a binary, (µ * α,A , µ δ,A ) and (µ * α,B , µ δ,B ),1 along with the accurately and precisely measured distances d A and d B , the magnitude of of the plane-of-sky relative velocity v p is given by (4) For Equation (4) to be used for actual data, the precision of d A and d B must be extremely good to prevent spurious boost of v p in some systems caused by random measurement errors.Thus, in practice it is more accurate to use where d is a representative distance in pc to the binary system, and with all PM values given units of mas yr −1 .For d I take an error-weighted mean (d M ) of d A and d B .Tests with the Gaia sample show that velocities estimated with Equations ( 4) and ( 5) are statistically equivalent only when the precision of distances is better than ε ≈ 0.002.In this work Equation ( 5) will be used because up to ε ≈ 0.005 is considered.Note also that because the distance range is 9 ≲ d/pc ≲ 200 (Figure 4) and 8 × 10 −6 ≲ s/d ≲ 2.8 × 10 −3 with a median of 5 × 10 −5 , it is sufficient to assume a plane geometry for the sky region of a binary system.
The uncertainty of the PM magnitude (Equations ( 6)) is estimated following El-Badry et al. (2021) as where The uncertainty of the sky-projected velocity is given by The normalized velocity parameter ṽ (Banik & Zhao 2018) and its uncertainty are given by ṽ ≡ v p /v c , σ ṽ ≡ ṽ where v c ≡ GM tot /s is the Newtonian circular velocity defined at the projected separation s (Banik & Zhao 2018) and I take σ vc /v c = 0.05 assuming a total mass uncertainty of 10%.The uncertainties of v p and ṽ are introduced here as additional means to check/control data quality.
Figure 5 shows the distributions of the estimated uncertainties of v p and ṽ for the pure binaries of the main sample shown in Figure 3. Virtually all pure binaries selected with strict astrometric and kinematic criteria have good signal-to-noise (S/N ≳ 3) with a median of about 40 for v p .Banik et al. (2023) advocate a cut based on σ ṽ < 0.1 max(1, ṽ/2).Figure 5 shows that this cut is already satisfied by virtually all pure binaries.Thus, for the pure binary sample I will not consider any artificial cut using either σ vp /v p or σ ṽ .In Appendix B, I will discuss the effects of cuts in general samples.
Gaia DR3 radial velocities typically have much less precision than PMs.Thus, most radial velocities cannot be used to measure the relative radial velocity v r between the two components because large random errors in individual radial velocities can create spurious boosts in many systems.However, for exceptionally precise radial velocities with the measurement precision comparable to that of ∆µ (Equation ( 6)), two relative velocities v p and v r can be combined to reliably estimate the relative physical velocity between the two stars Considering that all PM components have relative errors < ε, the relative error of v p,i (i = α, δ)2 is < √ 2ε for ε = 0.005.To require that the relative error of v r is comparable with that of v p,i , I require relative errors of individual radial velocities < 0.005.Finally, I note that gravitational redshifts (El-Badry 2022) from the surface gravities of the stars are irrelevant for the stars used in this work because stellar mass-to-radius ratio M/R varies little for the mass range 0.1 ≲ M/M ⊙ ≲ 1 (Demory et al. 2009).

A MONTE CARLO METHOD OF TESTING GRAVITY WITH STACKED VELOCITY PROFILES OF PURE BINARIES
Testing gravity in an acceleration plane with pure binaries will be done using the algorithm of Chae (2023) with f multi = 0.Here I describe a Monte Carlo method of testing gravity with stacked velocity profiles of pure binaries.The description of elliptical orbits will be the same as that of Chae (2023).The pertinent question is how to predict the sky-projected relative velocity v p (s) and the physical relative velocity v(s) between the two orbital plane  2023)) The left panel shows a one-particle equivalent description of orbital motions of the two stars in a binary system.The right panel defines the observer's viewpoint at an inclination i.
stars with a sky-projected separation s that can be compared with the observed velocities, Equations ( 5) and (11).
Figure 6 shows an equivalent one-body description of the elliptical orbit of binary dynamics taken from Chae (2023).The orbit is described in the plane polar coordinates (r, ϕ) by the equation where e is the eccentricity, a is the semi-major axis, and ϕ 0 is the longitude of the periastron.
In Newtonian dynamics, the magnitude of the relative physical velocity between the two stars is given by Physical separation r is related to the sky-projected separation s by where i is the inclination and ϕ is the azimuthal angle of the physical separation vector on the orbital plane as shown in Figure 6.
Combining Equations ( 12), (13), and ( 14), we can express the magnitude of the relative physical velocity as a function s, The magnitude of the sky-projected velocity to the observer is given by where (I note that the factor π is physically irrelevant and added to match the definition given in Figure 6 exactly.)For the observed set of (M tot , s), one MC realization of Newtonian velocities of Equations ( 15) and ( 16) follow from MC realizations of ϕ 0 , ϕ, i and e.Because possible ranges of these parameters are broad, the predictions for one binary system cannot be meaningfully compared with the observed velocities to test gravity.However, if a number of binary systems are considered simultaneously, the individual random fluctuations are averaged out and thus the mean of the predictions can be meaningfully compared with the mean of the observed velocities.Moreover, if MC realizations are repeated many times, on can derive the probability distribution of the mean in a sample and thus estimate its statistical uncertainty.This procedure allows one to test gravity in a quantitative way.
MC realizations of ϕ 0 , ϕ, i and e follow those described in Chae (2023).They can be summarized as follows: (1) ϕ 0 is drawn randomly from the range (0, 2π); (2) ϕ comes from the time along the orbit randomly drawn from (0, T ) where T is the period; (3) i is randomly drawn from (0, π/2) with a probability density p(i) = sin(i); finally, (4) e is drawn from the individual ranges provided by Hwang et al. (2022) as shown in figure 8 of Chae (2023).Note particularly that the eccentricity distribution for each binary is specified by three values: the most likely value (e m ), a lower-bound value (e l ), and an upper-bound value (e u ).In an MC, eccentricity is drawn using a combination of two truncated Gaussian functions: e m is taken as the median and each side is assumed to follow a truncated Gaussian function with a "σ" of e u −e m or e m −e l with the total range bounded by the limit 0.001 < e < 0.999.I also consider a power-law distribution p(e; α) = (1 + α)e α with a systematically varying with x ≡ log 10 (s/au) based on figure 7 and table 1 of Hwang et al. (2022).Equation ( 18) is valid for the range 2 ≲ x ≲ 4.5.For further details, the reader is referred to Chae (2023).

RESULTS
The sample of 2,463 pure binaries defined in Section 2 was obtained from a systematic investigation using the code developed in Chae (2023) and revised as described in Appendix A. Various samples defined in Chae (2023) always had f multi > 0 when binary motions were fitted to the Newtonian expectation in the high-acceleration regime ≳ 10 −8 m s −2 .This means that generally defined samples always include undetected close companions as widely appreciated.Also, the values of f multi ≈ 0.2−0.5 obtained in Chae (2023) with the revision of Appendix A agree well with the results from various surveys (Raghavan et al. 2010;Tokovinin 2014;Riddle et al. 2015;Moe & Stefano 2017) indicating that the calibration procedure of Chae ( 2023) is reliable.
In Section 4.1, I show the results in the acceleration plane for the sample of pure binaries obtained with the code of Chae (2023).The results clearly show that binary motions match the Newtonian expectation with f multi = 0 in the high-acceleration regime.The results will then provide a new test of gravity in the lowacceleration regime.In Section 4.2, I present the main results of this work, i.e., the stacked velocity profiles compared with the MC predictions of Newtonian gravity.

Results on the acceleration plane
Figure 7 shows one MC result for the kinematic acceleration g ≡ v 2 /r defined in Chae (2023) against the Newtonian gravitational acceleration g N ≡ GM tot /r 2 for the main sample of 2,463 pure binaries.Here stellar masses are estimated through the standard magnitudemass (M G -M ) relation (the first choice in table 1 of Chae ( 2023)) that is based on the Pecaut & Mamajek (2013) V -band magnitude-mass relation.The Newtonian expectation of the g N -g relation is compared with that for the Gaia data.In particular, the median orthogonal deviations ⟨∆ ⊥ ⟩ in the acceleration bins are quantitatively compared, as shown in the bottom panels of Figure 7.
Another MC gives different ⟨∆ ⊥ ⟩ values, and probability distributions of ⟨∆ ⊥ ⟩ can be derived from a number of MC results as demonstrated in Chae (2023).Figure 8 shows the distributions from 400 MC results.Those shown in the left column are from the MC results with the standard M G -M relation while those in the right are with the J-band based M G -M relation (the second choice in table 1 of Chae ( 2023)).It is clearly seen that the Gaia result naturally matches the Newtonian expectation in the highest acceleration bin at x 0 ≈ −8 with f multi = 0 (i.e.without any calibration) whichever M G -M relation is used.The parameter δ obs−newt ≡ ⟨∆ ⊥ ⟩ observed − ⟨∆ ⊥ ⟩ Newton defined in paper I has values of δ obs−newt = −0.004± 0.019 and −0.009 ± 0.021 at x 0 ≈ −8, well consistent with zero.This agreement is remarkable considering that there are no free parameters.
To check that the above agreement with the Newtonian expectation at x 0 ≈ −8 in the main sample is true rather than a coincidence, I consider subsamples selected with more stringent astrometric requirements (along with the same requirement on radial velocities).When relative errors of PMs and parallaxes are required to be < ε, I consider ε = 0.004 and 0.0025.The MC results on the acceleration plane for these subsamples can be found in Figure 9.These results agree well with δ obs−newt = 0 at x 0 ≈ −8 though with larger statistical errors due to smaller sample sizes.
The above results are based on individual eccentricities estimated by Hwang et al. (2022) and thus represent currently most likely results.Figure 10 shows a result based on eccentricities drawn statistically from a powerlaw distribution with a varying index given by Equation (18).Because binary-specific eccentricities are replaced by statistical eccentricities, the deviation is somewhat diluted as already noted in Chae (2023).However, the result still favors the AQUAL model over Newton.
So far I have considered three bins so that each bin has a maximal number of MC points for bins of significantly different accelerations.Now I consider finer bins to test any dependence of the results with binning.Figure 11 shows the results.The result with the standard inputs including individual eccentricities follows the AQUAL curve remarkably well over the entire bins.
If the astrometric requirements or the requirement on radial velocities are relaxed progressively from PM and parallax relative errors < 0.005 or RV relative errors < 0.2, one can check that δ obs−newt at x 0 = −8 with f multi = 0 increases progressively from zero.Note that the astrometric and RV requirements can be less stringent than the presently chosen requirements depending on the tolerance about the value of δ obs−newt at x 0 = −8 with f multi = 0.For example, one could allow δ obs−newt to be consistent with zero only within the MC estimated 1σ.Here I am very conservative and require δ obs−newt to be consistent with zero within a small fraction of the MC estimated 1σ.
The above results verify that the main sample with the presently chosen astrometric and RV requirements is statistically free of hierarchical systems.Of course, there still can be a few individual systems that have minor undetected close companions.However, even if they are present, they are statistically negligible.The above results also reassure that the whole procedure, the Gaia data, and the empirical M G -M relations are all reliable.
Now the derived values of δ obs−newt at lower accelerations are not consistent with zero indicating that Newtonian gravity breaks down in the low-acceleration regime.Because the same astrometric and RV requirements are imposed on all binaries regardless of the separation s, it is unreasonbale to imagine that only more widely (s ≳ 2 kau) separated binaries preferentially have large amounts of undetected companions while the less widely (s ≲ 1 kau) separated binaries have none.Thus, these results provide robust evidence for two aspects of gravity: (1) Newtonian gravity holds for acceleration ≳ 10 −8 m s −2 , (2) Newtonian gravity (and thus general relativity) breaks down in the low acceleration regime ≲ 10 −9 m s −2 .
While the evidence from statistically pure binaries is robust, the statistical significance is much weaker than in Chae (2023) due to the much smaller sample size.From the results with the standard M G -M relation, δ obs−newt = 0.122 ± 0.041 at x 0 ≈ −10.3 with a significance of 3.0σ.At x 0 ≈ −9.1, δ obs−newt = 0.033 ± 0.022 with a significance of 1.5σ.Taken together these results rule out Newtonian gravity with a > 3σ confidence.At x 0 ≈ −10.3, the acceleration boost factor is where g obs = v 2 /r is the MC realized kinematic acceleration from the Gaia data, and g pred is the corresponding Newtonian prediction.This value is in excellent agreement with the values obtained in Chae (2023) for general samples with f multi self-calibrated.

Profiles of stacked velocities
Here I present profiles of stacked velocities for the pure binary sample and compare them with the Newtonian MC predictions as described in Section 3. I pay most attention to sky-projected velocity (v p ) profiles because radial velocities (v r ) as precise as v p are quite rare.However, I will also consider, for the first time, profiles of physical velocities v (Equation ( 11)) for a few dozens of binaries with exceptionally precise v r measurements.

Testing the general sample including impure binaries
Before presenting results on pure binaries I test the general sample of 26,615 binaries defined in Chae (2023) that includes "impure" binaries hosting unresolved hidden close companions.Figure 12 compares the observed sky-projected velocities with the Newton-predicted values without taking into account hidden close companions.As expected, the observed velocities are higher than the Newton-predicted values regardless of the separation between the two stars.This is a clear indication that the observed velocities are largely affected by additional masses from hidden companions.
However, Figure 12 also reveals that the observed binned medians do not follow the Keplerian scaling of ∝ s −0.5 .The measured scaling has a slope of −0.437 ± 0.003 exhibiting a 21σ deviation from −0.5.This indicates that f multi must increase sharply with s to be consistent with the Keplerian scaling if gravitational anomaly or other factors are not permitted.I note that any bias arising from distances cannot explain this slope because a subsample with a narrow distance range shows a similar slope.Thus, we are in a situation where "general" binaries at the same distances selected with the same criteria show stong differences in kinematics depending only on s.This motivates investigations of "pure" binaries in the following sections.

Main results
Figure 13 shows the profile of stacked v p values for all 2,463 pure binaries in the main sample and compares it with the Newtonian prediction with the standard M G -M relation.Small red dots represent individual observed values while the big red dots represent the median, the 16th percentile, and the 84th percentile in the bins of s.Small blue dots represent an output from one Newtonian MC run.The mean values of the median, the 16th percentile, and the 84th percentile in the bins of s are obtained from 400 MC runs.The big blue dots and errorbars shown in the left panel of Figure 13 represent the distributions from MC runs.Note that the Newtonian MC predicted velocities can have tangible uncertainties when numbers are small as in the relatively larger-s bins.Histograms in the right panels of Figure 13 show the distributions of log 10 (v p,obs /v p,N ) in the 6 bins.Note that the widths of the distributions are entirely determined by the distributions of vp,N .
Figure 13 shows that in the three smallest-s bins the observed median velocities agree well with the Newtonpredicted values with log 10 (v p,obs /v p,N ) = −0.008± 0.015, −0.001 ± 0.013, and 0.005 ± 0.016.Moreover, the 16th and 84th percentiles of the observed velocities agree well with the Newton-predicted values for the first two bins indicating that the observed distributions are fully consistent with the Newton-predicted distributions.Thus, for binary systems with 0.5 ≲ M tot /M ⊙ ≲ 2 Newtonian dynamics holds for s ≲ 2 kau as expected from Figures 1 and 2. This result is consistent with the results on the acceleration plane presented in Section 4.1.
The Newton-predicted median velocities follow the Keplerian profile ∝ s −1/2 as expected because the median mass varies little with s in the sample as shown in Figure 3.However, the observed median velocities do not follow the Newtonian profile for the probed entire s range.There is a jump in the vp,obs profile around s ≈ 2 kau.The deviation from the Newtonian predictions in the 4th bin is log 10 (v p,obs /v p,N ) = 0.045 ± 0.019 at log 10 (s/au) = 3.53.However, in the last two bins the deviations are 0.076 ± 0.027 and 0.085 ± 0.040 at log 10 (s/au) = 3.89 and 4.22.These deviations together represent a ≈ 5.0σ detection of gravitational anomaly.Hereafter the statistical significance is estimated as follows.For the last three bins that deviate from the solid blue line, the probabilities p(x < 0) (where x ≡ log 10 (v p,obs /v p,N )) are calculated based on the estimated µ and σ values and the product of the three probabilities is obtained.Then, the product value is used to estimate a Gaussian equivalent significance.
In the two largest-s bins of Figure 13 the boost factor for projected velocities is Because the kinematic acceleration was defined to be v 2 /r, γ vp is expected to match √ γ g if projection effects are averaged out in a statistical sample.This is indeed the case from Equations ( 19) and (20).
A qualitatively important aspect of the stacked velocity profile shown in Figure 13 is that the last two bins follow the Keplerian scaling ∝ s −1/2 indicating that they are pseudo-Newtonian with a boosted gravity as exactly predicted by MOND gravity under an EFE as shown in Figure 1.Thus, Figure 13 agrees well with the unique trait of MOND gravity.
Figure 14 shows the result with the J-band-based M G -M relation.This result is consistent with the result with the standard M G -M relation indicating that the results are robust within the currently available M G -M relations.The significance of the gravitational anomaly from the last three bins is 4.9σ.However, below I will further probe the effects of systematically varying the observed Gaia magnitudes.
The above results are for a mass range 0.5 ≲ M tot /M ⊙ ≲ 2.5.Now I consider subsamples with a limited mass range 1.1 < M tot /M ⊙ < 1.8. Figure 15 shows the results for the subsample with 1.1 < M tot /M ⊙ < 1.8 and the clean magnitude range 4 < M G < 14.
Figure 16 shows the results for the subsample with 1.1 < M tot /M ⊙ < 1.8 and a narrower magnitude range 4 < M G < 10.These results agree well with the results with the main sample.The statistical significance of the deviations is ≈ 4.0σ for both results.
It is interesting to consider subsamples obtained in the limiting cases of extreme precision of PMs and distances.Figure 17 shows the result for 1, 206 pure binaries with ε = 0.0025, i.e. twice better precision.The result is well consistent with the results for the main sample.The statistical significance of the deviations in the last three bins taken together is 3.9σ.
Finally, I consider statistical eccentricities based on the power-law distribution with the slope systematically varying with s as given by Equation ( 18), rather than individual ranges of eccentricities reported by Hwang et al. (2022).The overall trend of the stacked velocity profile agrees well with the results with individual eccentricities.The statistical significance of the deviations in the last 3 bins taken together is 4.1σ.The boost factor estimated based on the two largest-s bins is slightly lower than the value given in Equation ( 20): The main key results are summarized in Table 1.The projected velocity boost factor is in the range 1.15 ≤ γ vp ≤ 1.24.The statistical significance of the gravitational anomaly is in the range (3.9σ, 5.0σ).

Auxiliary analyses
In Section 4.2.2, the Newtonian analysis with f multi = 0 showed that the observed velocities in the smaller-s bins matched the Newton-predicted values while those in the larger-s bins were higher than the Newton-predicted values.It is interesting to explore any possibility of attributing the boosted velocities in the larger-s bins somehow to additional masses from hidden close companions or a systematically shifted magnitude-mass relation.Here I consider varying f multi to make the observed velocities in the large-s bins agree with the Newtonpredicted values.5)) while blue ones are Newtonian velocities (Equation ( 16)) in one MC realization.Big dots indicate the median, the 16th percentile, and the 84th percentile velocities in the bins.The error bars of the big blue dots are estimated from 400 MC realizations.The Newtonian median velocities vp,N follow the Keplerian scaling of vp,N ∝ s −1/2 .The observed median velocities vp,obs in the three lowest-s bins naturally match the Newtonian predictions.However, the observed median velocities deviate from the Newtonian predictions in the larger-s bins.The dotted line indicates the boosted velocity in the two largest-s bins.The right panels show the probability distributions of log 10 (v p,obs /vp,N) derived from the MC realizations.
I use the procedure of modeling close companions described in Chae (2023).Figure 19 shows the result with f multi = 0.5.In this case, the observed median velocities in the two largest-s bins match well the Newtonpredicted medians.However, the Newton-predicted distributions are broader as revealed by the 16th and 84th percentiles due to the scatters arising from added components.More importantly, the Newton-predicted velocities in the three smallest-s bins now severely deviate from the observed velocities.Thus, we would be in a Newtonian world where binaries in the three smallest-s bins require f multi ≈ 0 while those in the two largest-s bins require f multi ≈ 0.5 although two subsamples satisfy identical photometric, astrometric, and radial velocity requirements.Since this huge difference in f multi between the small-and large-s bins is ad-hoc and unlikely, this result reinforces the results of Section 4.2.2.
I also consider a pseudo-Newtonian analysis with a rescaled Newton's constant.This analysis is motivated because the results of Section 4.2.2 suggest that binaries      with s ≳ 5 kau follow pseudo-Newtonian dynamics.Figure 20 shows the result with G ′ = 1.44G.The observed median velocities in the two largest-s bins agree well with the pseudo-Newtonian predictions.The observed 84th and 16th percentiles also match well the pseudo-Newtonian predictions within the statistical uncertainties.These results suggest that binaries with s ≳ 5 kau truly obey pseudo-Newtonian dynamics.However, the velocities in the three smallest-s bins deviate severely from the pseudo-Newtonian predictions as well expected because they are in the Newtonian regime.

Physical velocity
Just 40 out of 2,463 binaries of the main sample have measured radial velocities satisfying the precision of < 0.005 required for proper motions in each dimension.These exceptional systems can be used to probe the profile of physical velocities in the 3D space.
Figure 21 shows the measured physical velocities along with the Newtonian predicted values.Only three bins of s are considered to obtain median velocities.As expected, the MC generated distributions of median velocities are quite broad.In the lowest-s bin the observed median velocity matches well the Newtonian prediction.In the middle bin the observed median velocity is consistent with the Newtonian prediction within 1.5σ.However, in the highest-s bin satisfying s > 3.5 kau there is an indication that vobs is higher than vN with log 10 (v obs /v N ) = 0.249 ± 0.089, which represents a 2.8σ deviation.This result is consistent with the results from the analyses of sky-projected velocity profiles.

DISCUSSION
All the above results except for those from the auxiliary analyses of Section 4.2.3 have been obtained without any free parameters.They are derived or calculated quite naturally from the Gaia measurements, the magnitude(M G )-stellar mass(M ) relations from Chae (2023), and the individual eccentricity ranges from Hwang et al. (2022).Are there any possibilities that any of the observational inputs or analyses are grossly in error?
Gaia data themselves cannot be a source of systematic error because only exceptionally precise data have been used in this study (see Figure 5) and the results remain consistent as the precision increases (see Figure 9 and 17).However, one possible concern may be that the above results are based on a composite of binaries at significantly different distances ranging from 9 pc up to 200 pc.Here I consider two subsamples in significantly narrower distance ranges of 50 < d M < 125 pc and 125 < d M < 200 pc as defined in Figure 4.The M G -M relation is also not like to be any source of systematic error.As shown in Chae (2023), the two M G -M relations reliably cover the likely range at least for the clean magnitude range 4 < M G < 14, and the results with the two relations give consistent results (see Figure 8 and compare Figure 13 with Figure 14).Moreover, the result for the subsample with a narrower magnitude range 4 < M G < 10 (Figure 16) is consistent with those for the full sample.(Note that the M G -M relation is particularly accurate in the range 4 < M G < 10 as shown in figure 7 of Chae (2023).)Nevertheless, I consider systematically shifting magnitudes by ±0.5 mag just as a wild possibility in the estimated G-band magnitudes.I note that masses of Pecaut & Mamajek (2013) agree well with the masses measured directly from confirmed close binaries at the same K-band magnitudes (Mann et al. 2019).Thus, it suffices to consider a systematic shift in G-band magnitudes only.Figure 24 shows the result with a shift of −0.5.Even in this case, the median velocities in the two largest-s bins deviate significantly with a velocity boost factor of ≈ 1.14 while the two smallest-s bins deviate in the opposite direction.
Kinematic contaminants such as hidden close companions cannot be a source of systematic error.Because binaries with small and large s values have similar total masses and satisfy the same selection criteria, if additional masses should be added to the binaries, similar masses must be added to the binaries regardless of their separations.Then, the fractional boost in velocity is similar for all binaries regardless of s because the boosted-to-initial velocity ratio v ′ /v = M ′ /M is independent of s.As shown in Figure 19, when additional masses are controlled by f multi , the two largest-s bins require f multi = 0.5 to match Newtonian dynamics while the three smallest-s bins require f multi = 0.This extremely unlikely difference means that it seems impossible to attribute the gravitational anomaly to differential kinematic contaminants of hidden close companions.This view is bolstered by the fitted values of f multi in subsamples with progressively stricter kinematic criteria as presented in Appendix B. When kinematic criteria get stricter and stricter, f multi gets lower and lower, approaching zero eventually.Thus, it is hardly expected to be as high as f multi = 0.5 in the statistically pure binary sample selected in this work.Because the sample  of 2,463 "pure" wide binaries is made public, this claim can be directly tested with observations (see Manchanda et al. 2023).
The individual eccentricity ranges from Hwang et al. (2022) are the best available empirical information on eccentricities at present.I have also considered the power-law distribution of eccentricities with s-dependent exponent (Equation 18).However, here I consider the uniform thermal probability distribution of eccentricity p(e) = 2e for all binaries to gauge possi-ble source of systematic error arising from eccentricities.Note that this choice is deliberately biased from the empirical information (see figure 24 of Chae ( 2023)).
Figure 25 shows the result with the thermal eccentricity distribution.In this case, the deviation is somewhat weakened but the combined statistical significance of the deviations in the two largest-s bins is 3.7σ.Thus, it is not possible to do away with the deviations by reasonably modifying eccentricities.

CONCLUSION AND FUTURE PROSPECTS
A sample of 2,463 statistically pure wide binaries in the mass range 0.5 ≲ M tot /M ⊙ ≲ 2 provides two crucial results on non-relativistic gravitational dynamics.
First, the observed orbital motions of binaries with relatively small sky-projected separations (s ≲ 2 kau) are statistically consistent with Newtonian dynamics, naturally and without any adjustment of Gaia data or other observational inputs.This is a nontrivial result that provides direct evidence that currently available gravity theories including Newtonian and Milgromian theories hold in the non-relativistic regime of accelerations ≳ 10 −8 m s −2 .
Second, the observed orbital motions of binaries with relatively larger s ≳ 2 kau are statistically inconsistent with Newtonian dynamics.The Gaia measured skyprojected velocities are boosted with a statistical significance of ≈ 5.0σ.In the bin of the investigated range 5 < s < 30 kau, the velocity boost factor is measured    to be γ vp = 1.20 ± 0.06 (stat) ±0.05 (sys) (see Table 1).When the pure binaries are analyzed in the acceleration plane defined in Chae (2023), the kinematic acceleration v 2 /r with MC-deprojected v and r is systematically higher in the low-acceleration regime ≲ 10 −9 m s −2 than the corresponding Newtonian prediction, with an acceleration boost factor of γ g = 1.49+0.21  −0.19 satisfying the expected relation γ g = γ 2 vp .In a small sample of 40 pure wide binaries with exceptionally precise Gaia radial velocities, the physical velocity (v) is directly measured (i.e.without any deprojection).Despite the small number statistics, it is seen that the mean measured velocity in the largest-s bin is boosted compared with the Newtonian prediction while that in the smallest-s bin naturally matches the Newtonian prediction.
The present results from analyses of statistically pure binaries provide a robust confirmation of the results of Chae (2023) for a much larger general sample that includes hierarchical systems with undetected companions.The present results also complement Hernandez (2023) who obtained a similar boost in projected velocities but with a lower or roughly quantified statistical significance.
Just recently and almost concurrently with this work, Hernandez et al. (2023) has carried out a statistical analysis of the Hernandez (2023) sample of 667 wide binaries within 125 pc.They obtained a boost factor of γ g = 1.512 ± 0.199 which is in good agreement with the results from this work and Chae (2023).
Unlike this work and other recent studies (Chae 2023;Hernandez 2023;Hernandez et al. 2023), Banik et al. (2023) claimed an opposite conclusion and argued particularly that the gravitational anomaly obtained by Chae (2023) was largely affected by kinematic contaminants.As shown in Figure 5, the sample of pure binaries used in this work already satisfies the Banik et al. (2023) kinematic cut and yet shows essentially the same gravitational anomaly reported by Chae (2023); Hernandez (2023); Hernandez et al. (2023).Moreover, as presented in Appendix B, new results with kinematic cuts imposed confirm the gravitational anomaly.
The evidence for the gravity boost in the low acceleration regime is now clear enough although the scientific community should keep gathering further evidence from future observations.What seems now more important is to precisely characterize the gravity boost to the point that the theoretical direction can be narrowed down.Given that precise radial velocities measured in just 40 binaries already show a mild evidence of the gravity boost in the low-acceleration regime (Figure 21), precise measurements of radial velocities in more pure binaries in the future may turn out to be quite fruitful in characterizing the gravity in the low acceleration regime.
In principle, theoretical interpretations of the gravitational anomaly obtained here and in Chae (2023) and Hernandez (2023) are wide open.However, the most straightforward interpretation at hand is that nonrelativistic gravitational dynamics is governed/described by MOND-type Lagrangian theories of gravity (Bekenstein & Milgrom 1984;Milgrom 2010Milgrom , 2023)).However, because MOND-type Lagrangian theories are nonrelativistic phenomenological theories, something like phenomenological quantization rules before the full development of quantum physics, it is unclear what will be the underlying fundamental theory that will explain the MOND phenomenology eventually.Because MOND breaks the strong equivalence principle (Chae et al. 2020(Chae et al. , 2021;;Chae 2022) while keeping the Einstein equivalence principle, even non-quantum gravity must be different from Einstein's general relativity (Einstein 1916) and reformulated encompassing the successful aspects of both MOND and general relativity perhaps in the spirit of Mach's principle.
Hypothetical dark matter was introduced as a solution to gravitational anomalies in the presently investigated low-acceleration regime in galaxies and galaxy clusters assuming that standard gravity holds in that regime.Now that standard gravity breaks down in the same low-acceleration regime regardless of hypothetical dark matter and in agreement with MOND, dark matter interpretation is seriously questioned as a valid solution.Thus, no direct detection of dark matter despite intensive worldwide campaigns can now be seen as a natural outcome.Because there has been no direct detection of dark matter, all circumstantial arguments and indirect "evidence" for dark matter assuming standard gravity can now be overridden by the present direct evidence for the breakdown of standard gravity.This means that the dark matter paradigm seems now doomed to be abandoned and we are entering an era of a paradigm shift.Implications of the gravitational anomaly in the low-acceleration regime for astrophysics, cosmology, and fundamental physics are truly far-reaching.
In particular, the standard cosmology based on general relativity seems no longer valid even in principle.Development of MOND-based cosmology and structure formation (e.g., Sanders 1998Sanders , 2001Sanders , 2008;;Wittenburg et al. 2023) is now well-motivated and much-needed in parallel with theoretical advancement of MOND (e.g.Thomas et al. 2023).

ACKNOWLEDGMENTS
The author thanks Kareem El-Badry for discussion on radial velocities of stars.The revised version was based on an insightful report for which the author thanks the referee, a plenary talk given at the Korean Astronomical Society Fall 2023 meeting, and an invited talk given at the Pacific Rim Conference on Stellar Astrophysics held at Sejong University in October 2023.In particular, the author spotted an error in the code described in Chae (2023), which is typographic in nature, while preparing for the plenary talk, and the correction was reflected in this revision.This work was supported by the National Research Foundation of Korea (grant No. NRF-2022R1A2C1092306).The right column is the same for a subsample selected with a tighter precision on PM and additional requirements on radial velocities.In the latter, I require that radial velocities of both components have S/N > 2 and match each other within 3σ.Note that the subsample largely satisfies the kinematic cuts indicated by vertical dashed lines.
Because wide binary test is intended to test whether the observed velocity or velocity squared over radius (i.e.kinematic acceleration) is boosted or not compared with the Newtonian prediction, I also consider a more direct cut based on the uncertainty of the observed sky-projected velocity v p .
The left column of Figure 28 shows the distributions of the two uncertainty parameters for the main sample of Chae (2023) with the PM error cut < 0.01.About 27% of binaries from the Chae (2023) main sample do not pass the cut defined by Equation (B2) (this is higher than 21% obtained by Banik et al. (2023) because my cut is somewhat stricter).For the Chae (2023) sample with the stricter PM error cut < 0.005, a smaller fraction of about 19% do not pass the cut.
But, even if the cut such as Equation (B2) can be fully justified, what effect would it have on inferring gravity from wide binaries?Banik et al. (2023) argued that a qualitative examination of the scaling of ṽ with sky-projected radius normalized by the MOND radius r M ≡ GM tot /a 0 could question the conclusion of Chae (2023).However, the question can only be answered through detailed quantitative analyses.Also, the Banik et al. (2023) cut itself needs to be examined.Thus, a thorough investigation could be done properly only through a heavy independent work.Here I present only the default results based on the cut given by Equation (B2) and the standard inputs of Chae (2023).
Figure 29 shows the results for the two main samples of Chae (2023).The deviations from the corresponding Newtonian predictions are only mildly lowered compared with those without the kinematic constraint (Figure 26).In particular, the result (the right panel) for the sample with PM error < 0.005 is barely affected by the cut.These results indicate that the Banik et al. (2023) argument based on the qualitative appearances of ṽ curves does not hold.However, these results are not surprising because much more precise data used in the main part of this work already agreed with the results of Chae (2023).
Figure 29 also shows that the fitted values of f multi are lowered if the kinematic constraint of Equation (B2) is imposed.This is easily understood by the fact that hidden close companions are likely to add some kinematic noise.The fitted values of f multi = 0.36 and 0.31 are quite reasonable considering that some noisy data were removed.Although the sample of Banik et al. (2023) is not the same as the samples used here, the extremely high value of f multi ≈ 0.70 for their preferred Newtonian model indicates that their preference of Newtonian gravity is likely to be due to the largely biased f multi .The bias may have occurred because they chose not to calibrate f multi using Newtonian regime data.
I also consider a more reliable subsample than the main samples shown in Figure 29.In this sample, it is required that S/N values of the reported radial velocities for both components are greater than 2 and two radial velocities match each other within 3σ (adjusting Equation ( 2)).This sample will be intermediate between the main samples of Chae (2023) and the most precise pure binary sample from this work.The right column of Figure 28 shows the distributions of the two uncertainty parameters for this sample.Interestingly, both kinematic constraints are largely satisfied by this sample.I carry out detailed tests of gravity with this sample.

Figure 1 .
Figure1.This figure shows an AQUAL numerical(Chae & Milgrom 2022) prediction on the internal gravitational field for circular orbits under the external field of the Milky Way as estimated inChae (2023).The numerical acceleration matches well the analytic asymptotic value when the external field is inclined at an average angle of 60 • from the orbital axis.Internal dynamics is expected to switch from the Newtonian regime to the MOND regime over the narrow acceleration range indicated by the cyan-colored band.

Figure 2 .
Figure 2.This figure shows the range of sky-projected separation s corresponding to the Newton-MOND transition range in acceleration shown in Figure 1.The inset shows the Keplerian prediction on the range of orbital periods for Gaia wide binaries used in Chae (2023) and this study.

Figure 3 .
Figure 3.The left panel shows distributions of individual masses and total masses in 2,463 pure binaries of the main sample defined in the text.The dashed black vertical lines define a narrow mass range of total masses 1.1 < Mtot/M⊙ < 1.8.The right panel shows the mean masses in the 6 bins defined by sky-projected separation s as indicated by vertical red dashed lines.These 6 bins will be used in the analyses of sky-projected velocities.

Figure 4 .
Figure 4.This figure shows the distribution of dM (the error-weighted mean distance) for the sample of statistically pure binaries shown in Figure 3. Two distance bins are indicated by the vertical dashed lines.These bins will be used to probe any possible systematic effect of distances.

Figure 5 .
Figure 5.The upper panel shows the distribution of σṽ (Equation (10)) while the bottom panel shows the distribution of the fractional uncertainty of the sky-projected velocity σv p /vp for the main sample of 2,463 pure binaries.The upper panel indicates the cut similar to that suggested by Banik et al. (2023).Note that virtually all pure binaries satisfy the Banik et al. (2023) cut and S/N > 2 for the skyprojected velocity.

Figure 6 .
Figure 6.(Adapted from Chae (2023)) The left panel shows a one-particle equivalent description of orbital motions of the two stars in a binary system.The right panel defines the observer's viewpoint at an inclination i.

Figure 7 .
Figure 7. MC realized distributions of 2,463 pure binaries in the acceleration plane defined in Chae (2023).The quantity gN ≡ GMtot/r 2 is the Newtonian gravitational acceleration between the two stars and g ≡ v 2 /r is an empirical kinematic acceleration, where r and v are deprojected 3D separation and relative velocity.The left panel shows a Newton-predicted distribution while the right panel shows a distribution from Gaia measurements from one MC realization.Big dots indicate the medians in the orthogonal bins indicated by magenta dotted lines.The orthogonal deviation ∆ ⊥ of a point from the diagonal line is indicated in the inset of the upper left panel.The bottom panels show the medians of ∆ ⊥ in the bins.In the right panels, the Gaia medians are compared with he Newtonian medians.

Figure 8 .
Figure 8.The upper left panel shows distributions of median orthogonal deviations ⟨∆ ⊥ ⟩ (defined in Figure 7) from 400 MC results with the standard MG-M relation for the main sample of 2,463 pure binaries.The upper right panel is with the J-band-based MG-M relation.The bottom panels show distributions of the difference δ obs−newt ≡ ⟨∆ ⊥ ⟩ obs − ⟨∆ ⊥ ⟩N.In the lowest acceleration bin at x0 ≈ −8.0, δ obs−newt = 0 is naturally satisfied without any adjustment of f multi .The magenta curve in the bottom panels represents the AQUAL prediction for circular orbits with the Milky Way external field (see Figure 1).

Figure 9 .
Figure 9.The same as the left panel of Figure 8 but for subsamples with stricter astrometric requirements than the main sample.The results have larger statistical uncertainties and are consistent with those for the main sample.

Figure 10 .
Figure10.The same as the left panel of Figure8but with statistical eccentricities drawn from a power-law distribution with a varying index given by Equation (18).

Figure 11 .
Figure 11.Similar to Figure 8 but with 7 bins.The left panel is with the standard inputs including individual eccentricities while the right panel is with statistical eccentricities replacing individual eccentricities.

Figure 12 .
Figure 12.The left panel shows sky-projected velocities with respect to sky-projected separation s for the general sample of 26,615 binaries with the standard MG-M relation.Small red dots are the observed velocities (Equation (5)) while blue ones are Newtonian velocities (Equation (16)) in one MC realization.Big dots indicate median, 16th, and 84th percentile velocities in the bins of s defined in Figure3.The Newtonian median velocities vp,N follow the Keplerian scaling of vp,N ∝ s −1/2 .The observed median velocities deviate from the Newtonian predictions in the entire range and the scaling as a slope clearly different from −1/2.The right panels show the probability distributions of log 10 (v p,obs /vp,N) derived from 400 MC realizations.

Figure 13 .
Figure 13.The left panel shows sky-projected velocities with respect to sky-projected separation s for the main sample of 2,463 pure binaries with the standard MG-M relation.Small red dots are the observed velocities (Equation (5)) while blue ones are Newtonian velocities (Equation (16)) in one MC realization.Big dots indicate the median, the 16th percentile, and the 84th percentile velocities in the bins.The error bars of the big blue dots are estimated from 400 MC realizations.The Newtonian median velocities vp,N follow the Keplerian scaling of vp,N ∝ s −1/2 .The observed median velocities vp,obs in the three lowest-s bins naturally match the Newtonian predictions.However, the observed median velocities deviate from the Newtonian predictions in the larger-s bins.The dotted line indicates the boosted velocity in the two largest-s bins.The right panels show the probability distributions of log 10 (v p,obs /vp,N) derived from the MC realizations.

Figure 14 .
Figure 14.The same as Figure 13 but with the J-band-based MG-M relation.

Figure 15 .
Figure 15.The same as Figure 13 but for the subsample with the limited mass range 1.1 < Mtot/M⊙ < 1.8.

Figure 16 .
Figure 16.The same as Figure 13 but for the subsample with the narrower magnitude range 4 < MG < 10 and the limited mass range 1.1 < Mtot/M⊙ < 1.8.

Figure 17 .
Figure17.The same as Figure13but for the subsample with a more stringent astrometric requirement of ε = 0.0025.

Figure 18 .
Figure18.The same as Figure13but with statistical eccentricities from an s-dependent power-law distribution.
show the results.Because sample sizes are smaller, these results have larger statistical uncertainties, but are consistent with each other and the main results.Note particularly that the 50 < d M < 125 pc sample has particularly larger statistical uncertainties because the two largest-s bins have relatively few binaries.

Figure 20 .
Figure 20.The same as Figure 13 but with a boosted Newton's constant G ′ = 1.44G.

Figure 21 .
Figure 21.The result for the physical velocity v = v 2 p + v 2 r for binaries with exceptionally precise radial velocities.The relative errors of individual radial velocities are required to be smaller than 0.005 as described in the text.

Figure 22 .
Figure 22.The same as Figure 13 but for a subsample within a narrow distance range of 50 < dM < 125 pc.

Figure 23 .
Figure 23.The same as Figure 13 but for a subsample within a narrow distance range of 125 < dM < 200 pc.

Figure 24 .
Figure 24.The same as Figure 13 but with G-band magnitudes systematically shifted by −0.5.

Figure 25 .
Figure25.The same as Figure13but with eccentricities from a 'thermal' probability distribution p(e) = 2e ignoring individual eccentricities.The result is deliberately biased toward a weakened boost of velocities.Nevertheless, the boost is still very significant.

Figure 26 .
Figure 26.Gravitational anomalies of Chae (2023) are revised for the main samples with a correction to the published code with Equation (A1).The upper panels show distributions of median orthogonal deviations ⟨∆ ⊥ ⟩ (defined in Figure 7) from 200 MC results.The upper right panel is with the J-band-based MG-M relation.The bottom panels show distributions of the difference δ obs−newt ≡ ⟨∆ ⊥ ⟩ obs − ⟨∆ ⊥ ⟩N.Parameter f multi was fitted through an iterative procedure so that δ obs−newt = 0 is satisfied at x0 ≈ −8.0.The magenta curve in the bottom panels represents the AQUAL prediction for circular orbits with the Milky Way external field.

Figure 27 .
Figure 27.Same as Figure 26 but with eccentricities drawn from the power-eccentricity distribution with the slope varying with s as given by Equation (18).

Figure 28 .
Figure28.The left column shows the distributions of two uncertainties σṽ and σv p /vp for the main sample ofChae (2023).The right column is the same for a subsample selected with a tighter precision on PM and additional requirements on radial velocities.In the latter, I require that radial velocities of both components have S/N > 2 and match each other within 3σ.Note that the subsample largely satisfies the kinematic cuts indicated by vertical dashed lines.

Figure 29 .
Figure 29.Same as Figure 26 but with the kinematic constraint of Equation (B2) imposed.
Figures 30 and 31 show the results from testing gravity with the subsamples shown in the right column of Figure 28.Two samples with different kinematic cuts imposed return similar results.The results with individual eccentricities match the AQUAL prediction remarkably well.It is also interesting to note that for these relatively more reliable samples (compared with the Chae (2023) main samples) there is a clear indication for the gravitational anomaly, though somewhat weakened, even with statistical eccentricities that are less informative than individual eccentricities.The fitted values of f multi = 0.15 or 0.19 are even lower than those shown in Figure 29, indicating further that the Banik et al. (2023) value of ≈ 0.70 is unlikely for samples satisfying kinematic cuts.Finally, Figures 32 and 33 show the results for seven bins from testing gravity with the subsamples used for Figures 30 and 31.The results with the more reliable individual eccentricities closely follow the AQUAL prediction for the entire range of acceleration covered by the data.

Figure 30 .
Figure 30.Similar to the right panel of Figure 29 but for a subsample with radial velocities shown in the upper-right panel of Figure 28.Note that statistical eccentricities are also considered to provide systematically biased results.

Figure 31 .
Figure 31.The same as Figure 30 but with a different kinematic cut shown in the lower-right panel of Figure 28.

Figure 32 .
Figure 32.Similar to Figure 30 but for finer bins.

Figure 33 .
Figure 33.The same as Figure 31 but for finer bins.