Lowacceleration gravitational anomaly is investigated with a new method of exploiting the normalized velocity profile inlineformula mmlmath overflowscrollmmlmover accenttruemmlmrowmmlmivmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmovermmlmommlmommlmsubmmlmrowmmlmivmmlmimmlmrowmmlmrowmmlmipmmlmimmlmrowmmlmsubmmlmrowmmlmo stretchytruemmlmommlmrowmmlmsubmmlmrowmmlmivmmlmimmlmrowmmlmrowmmlmicmmlmimmlmrowmmlmsubmmlmath inlineformula of wide binary stars as a function of the normalized skyprojected radius sr SUBMSUB where v SUB p SUB is the skyprojected relative velocity between the pair v SUB c SUB is the Newtonian circular velocity at the skyprojected separation s and r SUBMSUB is the MOND radius. With a Monte Carlo method Gaia observed binaries and their virtual Newtonian counterparts are probabilistically distributed on the sr SUBMSUB versus inlineformula mmlmath overflowscrollmmlmover accenttruemmlmrowmmlmivmmlmimmlmrowmmlmrowmmlmommlmommlmrowmmlmovermmlmath inlineformula plane and a logarithmic velocity ratio parameter  is measured in the bins of sr SUBMSUB. With three samples of binaries covering a broad range in size data quality and implied fraction of hierarchical systems including a new sample of 6389 binaries selected with accurate distances and radial velocities I find a unanimous systematic variation from the Newtonian flat line. With   0 at sr SUBMSUB  0.15 or s  1 kau I get   0.068  0.015 stat inlineformula mmlmath overflowscrollmmlmsubsupmmlmrowmmlmrowmmlmrowmmlmommlmommlmn0.015mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.024mmlmnmmlmrowmmlmsubsupmmlmath inlineformula syst for sr SUBMSUB  0.7 or s  5 kau. The gravitational anomaly i.e. acceleration boost factor given by  SUB g SUB  10SUP2SUP is measured to be inlineformula mmlmath overflowscrollmmlmsubmmlmrowmmlmimmlmimmlmrowmmlmrowmmlmigmmlmimmlmrowmmlmsubmmlmommlmommlmsubsupmmlmrowmmlmn1.37mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.09mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.10mmlmnmmlmrowmmlmsubsupmmlmath inlineformula stat inlineformula mmlmath overflowscrollmmlmsubsupmmlmrowmmlmrowmmlmrowmmlmommlmommlmn0.09mmlmnmmlmrowmmlmrowmmlmommlmommlmn0.16mmlmnmmlmrowmmlmsubsupmmlmath inlineformula syst. With a reduced  SUP2SUP test of Newtonian and Milgromian nonrelativistic theories I find that Newtonian gravity is ruled out at 5.8 inlineformula mmlmath overflowscrollmmlmsubsupmmlmrowmmlmimmlmimmlmrowmmlmrowmmlmimmlmimmlmrowmmlmrowmmlmn2mmlmnmmlmrowmmlmsubsupmmlmommlmommlmn9.4mmlmnmmlmath inlineformula by the new sample and 9.2 by the largest sample used. The Milgromian AQUAL theory is acceptable with inlineformula mmlmath overflowscrollmmlmn0.7mmlmnmmlmommlmommlmsubsupmmlmrowmmlmimmlmimmlmrowmmlmrowmmlmimmlmimmlmrowmmlmrowmmlmn2mmlmnmmlmrowmmlmsubsupmmlmommlmommlmn3.1mmlmnmmlmath inlineformula. These results agree well with earlier results with the accelerationplane analysis for a variety of samples and the stacked velocity profile analysis for a pure binary sample.
['Gravitation', 'Wide binary stars', 'Modified Newtonian dynamics', 'Non-standard theories of gravity', 'Binary stars', '661', '1801', '1069', '1118', '154', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Theory']
Measurements of the Lowacceleration Gravitational Anomaly from the Normalized Velocity Profile of Gaia Wide Binary Stars and Statistical Testing of Newtonian and Milgromian Theories
{'No Header': 'Measurements of the Low-Acceleration Gravitational Anomaly from the Normalized Velocity Profile of Gaia Wide Binary Stars and Statistical Testing of Newtonian and Milgromian Theories  \nKyu-Hyun Chae 1  \n1 Department of Physics and Astronomy, Sejong University, 209 Neungdong-ro Gwangjin-gu, Seoul 05006, Republic of Korea', 'ABSTRACT': "Low-acceleration gravitational anomaly is investigated with a new method of exploiting the normalized velocity profile ˜ v ≡ v p /v c of wide binary stars as a function of the normalized sky-projected radius s/r M where v p is the sky-projected relative velocity between the pair, v c is the Newtonian circular velocity at the sky-projected separation s , and r M is the MOND radius. With a Monte Carlo method Gaia observed binaries and their virtual Newtonian counterparts are probabilistically distributed on the s/r M versus ˜ v plane and a logarithmic velocity ratio parameter Γ is measured in the bins of s/r M . With three samples of binaries covering a broad range in size, data quality, and implied fraction of hierarchical systems including a new sample of 6389 binaries selected with accurate distances and radial velocities, I find a unanimous systematic variation from the Newtonian flat line. With Γ = 0 at s/r M ≲ 0 . 15 or s ≲ 1 kilo astronomical units (kau), I get Γ = 0 . 068 ± 0 . 015 (stat) +0 . 024 -0 . 015 (syst) for s/r M ≳ 0 . 7 or s ≳ 5 kau. The gravitational anomaly (i.e. acceleration boost) factor given by γ g = 10 2Γ is measured to be γ g = 1 . 37 +0 . 10 -0 . 09 (stat) +0 . 16 -0 . 09 (syst). With a reduced χ 2 test of Newtonian and Milgromian nonrelativistic theories, I find that Newtonian gravity is ruled out at 5 . 8 σ ( χ 2 ν = 9 . 4) by the new sample (and 9 . 2 σ by the largest sample used). The Milgromian AQUAL theory is acceptable with 0 . 7 ≲ χ 2 ν ≲ 3 . 1. These results agree well with earlier results with the 'acceleration-plane analysis' for a variety of samples and the 'stacked velocity profile analysis' for a pure binary sample.  \nKeywords: :Binary stars (154); Gravitation (661); Modified Newtonian dynamics (1069); Non-standard theories of gravity (1118); Wide binary stars (1801)", '1. INTRODUCTION': "The nature of gravity in the low acceleration limit can be directly probed by wide binaries (widely separated, long-period, gravitationally bound binary stars) since any conceivable dark matter density (even if it were detected) in the Milky Way cannot affect their internal dynamics (e.g., Hernandez et al. 2012; Banik & Zhao 2018; Pittordis & Sutherland 2018; Banik & Kroupa 2019; Pittordis & Sutherland 2019; Hernandez et al. 2019; El-Badry 2019; Clarke 2020; Hernandez et al. 2022; Pittordis & Sutherland 2023; Hernandez 2023).  \nFor this reason and thanks to the state-of-the-art data provided by the Gaia satellite several statistical analyses of wide binaries have been recently carried out aiming at testing gravity quantitatively and even definitely (Chae  \nchae@sejong.ac.kr, kyuhyunchae@gmail.com  \n2023a; Banik et al. 2024; Chae 2024a; Hernandez et al. 2024) based on the Gaia data release 3 (DR3; Vallenari et al. 2023).  \nChae (2023a) put forward an algorithm (Chae 2023b) that calculates the probability distribution of a kinematic acceleration g = v 2 /r with respect to the Newtonian gravitational acceleration g N between the two stars where v is the relative velocity and r is the separation in the three-dimensional real space, and compares it with the corresponding Newtonian prediction. One key aspect of this algorithm is that the occurrence rate ( f multi ) of multiplicity higher than two (i.e. harboring hidden additional components) can be self-calibrated at a high enough acceleration and checked at another high acceleration. Through this algorithm, to be referred to as the 'acceleration-plane analysis', Chae (2023a) found that the observed acceleration started to get boosted from the Newtonian prediction for g N ≲ 10 -9 m s -2 with a boost factor of ≈ 1 . 4 for g N ≲ 10 -10 m s -2 , at an extremely high ( > 5 σ ) statistical significance. Chae (2023a) con-  \nvarious samples by varying selection criteria and noted that the low-acceleration gravitational anomaly persisted.  \nChae (2024a) further considered a sample of statistically pure binaries (i.e. the limiting case of f multi = 0) for an independent test. For this sample Chae (2024a) employed another algorithm to be referred to as 'stacked velocity profile analysis' as well as the accelerationplane analysis. The stacked velocity profile analysis compares the observed distribution of the sky-projected relative velocities against the sky-projected separations with the corresponding Newtonian prediction. Chae (2024a) found that the results for the pure binary sample through the two independent algorithms agreed well with each other and the Chae (2023a) results for general or 'impure' samples. This means that the gravitational anomaly at low acceleration is robust with respect to a large variation of the sample selection between f multi ∼ 0 . 5 (Chae 2023a) and f multi = 0 (Chae 2024a). Hernandez et al. (2024) carried out a statistical analysis of the distribution of normalized velocities in a pure binary sample that was selected by Hernandez (2023). Here the normalized velocity on the sky plane is defined (Pittordis & Sutherland 2018; Banik & Zhao 2018) by  \n˜ v ≡ v p v c ( s ) , (1)  \nwhere v p is the observed sky-projected relative velocity between the pair and v c ( s ) is the theoretical Newtonian circular velocity at the sky-projected separation s . The Hernandez et al. (2024) algorithm and sample are completely independent of Chae (2023a) and Chae (2024a) though the sample was drawn from the same Gaia DR3 database. Hernandez et al. (2024) obtained a gravitational anomaly that was well consistent with those obtained by Chae (2023a) and Chae (2024a).  \nThese recent quantitative analyses (Chae 2023a, 2024a; Hernandez et al. 2024) consistently show that the gravitational acceleration is boosted by a factor from ≈ 1 . 3 -1 . 5 and the relative velocity is boosted by a factor of ≈ 1 . 2 for g N ≲ 10 -10 m s -2 or s ≳ 5 kau (kilo astronomical units). Despite these consistent results from various samples of wide binaries with various independent methods, Banik et al. (2024) claimed an opposite conclusion based on their own statistical method and raised a concern for data quality control in Chae (2023a) samples.  \nHowever, Hernandez et al. (2023) soon pointed out the problems of their methodology. Above all, they knowingly excluded the Newtonian-regime ( g N > 10 -9 m s -2 or s < 2 kau) binaries that are essential for an accu-  \nrate determination of f multi in any sample. 1 Then, they employed a statistical method of fitting the distribution (or number count) of ˜ v (Equation (1)) in cells in an attempt to simultaneously constrain gravity, f multi , and other parameters without the Newtonian-regime data. However, in such a approach they improperly used cells smaller than the errors of ˜ v although their approach relies on number counts in cells. The reader is referred to Hernandez et al. (2023) for further details.  \nThe concern for data quality control raised by Banik et al. (2024) based on the uncertainty of ˜ v is not relevant for the acceleration-plane analysis of Chae (2023a) because ˜ v is not used, but probability distributions of g and g N are directly derived on the acceleration plane with sufficiently precise projected relative velocities, and taking into account possible ranges of parameters such as eccentricity, inclination, and orbital phase in the deprojection to the three-dimensional space. At any rate, the pure binary sample of Chae (2024a,b) already meets high signal-to-noise ratios ( S/N ≳ 10) for ˜ v . Moreover, Appendix B of Chae (2024a) shows that even with an artificial error cut of ˜ v imposed on the Chae (2023a) main samples the acceleration-plane analysis returns similar results for the gravitational anomaly for g N ≲ 10 -9 m s -2 .  \nAlthough the Banik et al. (2024) criticism of the Chae (2023a) acceleration-plane results based on their error cut of ˜ v is baseless, it is interesting to investigate the profile of ˜ v with respect to the projected separation s normalized by the MOND radius  \nr M ≡ √ GM tot a 0 (2)  \nwhere G is Newton's constant, M tot is the total mass of the system, and a 0 is the critical acceleration whose value is taken to be 1 . 2 × 10 -10 ms -2 . Banik et al. (2024) showed that the observed ˜ v profiles for the Chae (2023a) main samples were significantly affected by a cut on the uncertainty of ˜ v , and used this fact to question the gravitational anomaly reported by Chae (2023a). However, they did not calculate how the Newtonian prediction of the ˜ v profile in a sample was affected by modifying the sample with such a cut and thus did not make a valid comparison.  \nHere I propose an analysis of the ˜ v ( s/r M ) profile, to be referred to as 'normalized velocity profile analysis', as a method different from the other methods employed  \nin Chae (2023a, 2024a). This analysis is different from the stacked velocity profile analysis in that both velocity and separation are normalized. It is different from the acceleration-plane analysis in that no deprojection to the 3D space is considered. This method is interesting in its own right as an intermediate method between the acceleration-plane and stacked velocity profile analyses. It can also be used to directly and precisely address the effect of a cut on the uncertainty of ˜ v .  \nIn Section 2, I describe the wide binary samples and the methodology. In Section 3, I present the results. In Section 4, I compare the results with relevant previous results and discuss possible systematic errors. I conclude in Section 5. In Appendix A, I present the accelerationplane analyses results from which relevant f multi values are derived and used in the main analyses. Appendix B addresses effects of the proper motion (PM) measurement uncertainties in the gravity test methods used in this work. Revised and new python scripts used for this work and the sample of binaries can be accessed at Chae (2023b). Most analyses of this work are based on the medians of ˜ v (Equation (1)) in various bins. Throughout this paper they are denoted by ⟨ ˜ v ⟩ . This work uses several samples of wide binaries with various implied fractions of hierarchical systems including the limiting sample with f multi → 0. The latter sample was obtained in a systematic way by Chae (2024a) and revised by Chae (2024b). It will be referred to as the 'pure binary' sample as in Chae (2024a). However, it should be understood that the sample may include a tiny fraction of individual exceptions that evaded the stringent restrictions that Chae (2024a,b) applied statistically.", '2.1. Binary Samples': "As in Chae (2023a, 2024a), this work is based on a sample of binaries within 200 pc from the Sun that are free of clusters, background pairs, and resolved (1 '' ) triples or higher-order multiples. The sample was drawn from the El-Badry et al. (2021) catalog of 1.8 million binary candidates within 1 kpc from the Sun. The sample 2 of 81,880 binary candidates defined by Chae (2023a) has the range 0 . 05 < s < 50 kau. It is not free of chancealignment (fly-by) cases or unresolved close companions. Thus, they need to be dealt with carefully in using the Chae (2023a) catalog of 81,880 binary candidates.  \nChance-alignment cases (i.e. false binaries) and binaries including undetected close companions (i.e. hierarchical systems) are two main issues in assembling a sam-  \nple and carrying out gravity tests with it. There can be an infinite number of ways of defining a sample with a varying degree of including chance-alignment cases and hierarchical systems, and in principle any sample is permissible as long as chance-alignment cases and hierarchical systems are consistently modeled for the sample.  \nThe El-Badry et al. (2021) catalog has unique features. They introduced and estimated a chance-alignment probability R for each system. For this they used abundant observational information including angular separation, distance, parallax difference uncertainty, local sky density, relative tangential velocity, parallax difference over error, and PM difference over error. In particular, for the relative tangential velocity, i.e. the skyprojected relative velocity between the pair v p , to be consistent with a gravitationally bound system, they required that v p be smaller than the Newtonian expectation for a 5 M ⊙ system (considering that most systems have masses well below this limit) with a tolerance of 3 σ or 6 σ combined uncertainties depending on the separation (their Equations (3) through (7)).  \nThe left panel of Figure 1 displays all 1.8 million candidates in the plane spanned by the sky-projected separation s and v p . Clearly, this base catalog includes chance-alignment cases above the red dashed line given by Equation (7) of El-Badry et al. (2021). The middle panel with R < 0 . 1 excludes chance-alignment cases but still includes many systems above the red dashed line. The right panel exhibits the base sample defined by Chae (2023a), which includes some chance-alignment cases as well as hierarchical systems. Here I note that the El-Badry et al. (2021) Equation (7) primarily plays the role of removing chance-alignment cases but may also remove some pronounced hierarchical systems.  \nConsidering the uncertainties of chance-alignment and hierarchical cases in any sample, it is necessary to consider a range of samples rather than relying on one specific sample and check whether gravity test results are robust against the sample variation. 3 Here taking advantage of the unique features of the El-Badry et al. (2021) sample, I remove chance-alignment cases but allow a varying fraction of hierarchical systems.  \nChae (2023a) considered a 'clean' range of 4 < M G < 14 4 where M G is the Gaia G -band absolute magnitude so that a relatively more reliable part of an empirical magnitude-mass ( M G -M ) relation could be used (see  \nFigure 1. The left panel exhibits all 1.8 million binary candidates assembled by El-Badry et al. (2021). The diagonal dashed red line indicates a kinematic criterion used (with tolerance of 3 -6 σ combined measurement uncertainties) by El-Badry et al. (2021) to identify chance-alignment cases. Indeed, most chance-alignment cases are above this line. The middle panel shows cases with low chance-alignment probabilities satisfying R < 0 . 1 where R is the parameter introduced by El-Badry et al. (2021). Note that a significant number of systems are still above the red dashed line. Those systems are more likely to be hierarchical systems hosting hidden close companion stars although hierarchical systems do exist even below the line. The right panel shows the sample within 200 pc from the Sun defined by Chae (2023a).  \n<!-- image -->  \nFigure 7 of Chae (2023a)). Chae (2023a) presented a couple of M G -M relations based on the Pecaut & Mamajek (2013) results. Since the two M G -M relations give similar results as shown in Chae (2023a, 2024a), here I consider the standard choice only (the first choice given in Table 1 of Chae (2023a)).  \nIn this work I consider two approaches to remove chance-alignment cases. The first approach is to use the cut R < 0 . 01 following (Chae 2023a). Specifically, I use a main sample of 19,716 binaries defined in Chae (2023a) by requiring that the PM relative (i.e. fractional) errors are < 0 . 005 for all components and 0 . 2 < s < 30 kau. 5 This sample will be referred to as the 'Chae (2023a) sample'. Table 1 gives a summary of the samples to be used in this work.  \nThe second approach is to use RVs as well as distances to remove chance-alignment cases for those with measured RVs. Because two stars orbiting each other must be at the same distance from the Sun up to the minor difference of the orbit size and have the same system RV up to the minor difference from the orbital relative motion 6 , a requirement of consistent distances and RVs can effectively remove chance-alignment cases detected  \non the sky. For this, I first select binaries with wellmeasured PMs, distances, and RVs. The required relative errors are respectively < 0 . 005, < 0 . 01, and < 0 . 5 for PMs, distances, and RVs. Then, I require that distances satisfy  \n| d A -d B | < √ 9( σ 2 d A + σ 2 d B ) + (6 s ) 2 , (3)  \nand RVs satisfy  \n| v r,A -v r,B | < √ 9( σ 2 v r,A + σ 2 v r,B ) + (∆ v max r, orbit ) 2 (4)  \nwith  \n∆ v max r, orbit = 0 . 9419 km s -1 √ M tot s × 1 . 3 × 1 . 2 , (5)  \nwhere M tot is the total mass of the system in units of solar mass. Equations (3) and (4) are similar to Equations (1) and (2) of Chae (2024a) differing only in that 3 σ (instead of 2 σ ) differences are allowed. A sample based on this selection will be referred to as the 'new sample' (Table 1).  \nThe new sample is interesting for a few reasons. First, the selection is not based on R that has been used from Chae (2023a). I note that this selection is considered not because the El-Badry et al. (2021) R has any problems in a statistical sense but to provide an independent sampling. It turns out that 99.5% of the binaries in the new sample satisfy R < 0 . 01 automatically, proving that two independent methods are reliable and consistent. Second, the new sample includes only binaries with measured RVs with S/N > 2 and thus system velocities are  \nTable 1. Samples of binaries used in this studyNote. (1) The Chae (2023a) limited sample is considered for the purpose of investigating/illustrating the effects of a limited dynamic range. (2) The pure binary sample needs to be understood as a statistical sample satisfying f multi = 0 in the Newtonian regime (it may include rare individual exceptions of hierarchical systems).  \nknown. Lastly, the new sample has the additional requirement on the precision of distances (i.e. parallaxes). Overall, this sample is intermediate between the Chae (2023a) main sample and the pure binary sample defined in Chae (2024a) that has more stringent requirements.  \nIn this work, I also use the pure binary sample of Chae (2024a,b) for the normalized velocity profile analysis. Finally, I consider a subsample from the 'Chae (2023a) sample' that has a limited separation range of 2 < s < 30 kau. This sample will be referred to as the 'Chae (2023a) limited sample'. This sample is considered only for the purpose of illustrating or investigating the effects of limiting the dynamic range because Banik et al. (2024) (see also Pittordis & Sutherland 2023) recently considered such a narrow dynamic range.", '2.2. Empirical Relations to Test Gravity: from Acceleration Relation to Normalized Velocity': 'For the observed R. A. ( α ) and decl. ( δ ) components of the PMs in a binary, ( µ ∗ α,A , µ δ,A ) and ( µ ∗ α,B , µ δ,B ) , 7 the plane-of-sky scalar relative velocity v p is given by  \nv p = 4 . 7404 × 10 -3 km s -1 × ∆ µ × d (6)  \nwhere d is the distance in pc to the binary system, and  \n∆ µ = [ ( µ ∗ α,A -µ ∗ α,B ) 2 +( µ δ,A -µ δ,B ) 2 ] 1 / 2 (7)  \nwith all PM components given in units of mas yr -1 . Once the system is known/determined to be gravitationally bound, d is given by an error-weighted mean of the measured distances of the two components.  \nThe uncertainty of ∆ µ (Equations (7)) is estimated as  \nσ ∆ µ = [ ( σ 2 µ ∗ α,A + σ 2 µ ∗ α,B )(∆ µ α ) 2 +( σ 2 µ δ,A + σ 2 µ δ,B )(∆ µ δ ) 2 ] 1 / 2 / ∆ µ, (8)  \nwhere  \n(∆ µ α ) 2 = ( µ ∗ α,A -µ ∗ α,B ) 2 , (∆ µ δ ) 2 = ( µ δ,A -µ δ,B ) 2 . (9)  \nThen, the uncertainty of v p (Equation (6)) follows from the uncertainty of ∆ µ .  \nThe normalized velocity ˜ v (Equation (1)) follows from the measured value of v p (Equations (6)) and the theoretical Newtonian circular velocity at s given by  \nv c ( s ) ≡ √ GM tot s , (10)  \nwhere M tot is the total mass of the system including any hidden additional component(s).  \nThe observed distribution of ˜ v with respect to s/r M can test the low-acceleration behavior of gravity. This can be understood as follows. Consider a binary with a scalar relative velocity v and a relative separation r in the real 3D space so that g = v 2 /r and g N = GM tot /r 2 = v 2 c ( r ) /r . We then have  \ng g N = ( v v c ( r ) ) 2 (11)  \nand  \ng N a 0 = ( r r M ) -2 . (12)  \nAs done in Chae (2023a), a direct test of gravity can be done by the distribution of g with respect to g N in the acceleration plane. This requires a deprojection from the observed quantities v p and s to v and r with a Monte Carlo (MC) method. This test is reminiscent of the radial acceleration relation (McGaugh et al. 2016) test in galaxies. Now one can alternatively consider the ratio g/g N with respect to g N , or equivalently v/v c ( r ) with respect to r/r M using Equations (11) and (12). Finally, without a deprojection one can use ˜ v as a proxy for v/v c ( r ) and s/r M for r/r M .', '2.3. Normalized Projected Velocities from Proper Motions: Measured versus Newtonian Simulation': "The observed normalized velocities ˜ v (= v p /v c ( s )) will be statistically compared with their Newtonian counterparts that are generated by a Monte Carlo method. Here  \n𝑦𝑦  \nFigure 2. Adapted from Chae (2023a). The left panel shows a one-particle equivalent description of the relative motion between the two stars in a binary system. The right panel defines the observer's viewpoint at an inclination i .  \n<!-- image -->", 'orbital plane (face-on view)': '𝑧𝑧  \n<!-- image -->', "observer's view (3D geometry)": 'I describe how normalized velocities are calculated for the Gaia observed PMs and mock Newtonian PMs. The methodology follows from Chae (2023a) with a minor correction described in Appendix A of Chae (2024a). I summarize the essential concepts and elements while referring the reader to the relevant further details in Chae (2023a, 2024a).', '2.3.1. Measured Proper Motions': 'The projected relative velocity v p (Equation (6)) follows from the observed PMs through Equation (7). The Newtonian circular velocity (Equation (10)) depends on the total mass that may include the mass(es) of any unresolved hidden component(s). In each MC, a randomly selected fraction of binaries with probability f multi are given additional components: 70% of the hierarchical occurrences are assumed to be triples while 30% are quadruples. The mass of a close companion is probabilistically assigned with the empirical distribution of magnitude difference between the host and the companion as described in Section 3.2 of Chae (2023a).', '2.3.2. Newtonian Proper Motions: without hidden companions': "Each wide binary was observed by Gaia only at one phase of the orbital motion (in its long period) of an unknown eccentricity at an unknown inclination. Lacking the detailed information on the orbit, it is not possible to obtain a Newtonian mock PM as a precise counterpart to the observed PM. One can obtain only a random realization out of broad possible ranges. Thus, an individual mock Newtonian PM cannot be individually  \ncompared with the observed counterpart PM. Only a sufficiently large collection of mock Newtonian PMs can be statistically compared with the counterpart collection of observed PMs.  \nKeeping the above statistical nature in mind, now I describe how a random Newtonian PM is obtained for the observed binary system with masses M A and M B and the sky-projected separation s . Consider the oneparticle equivalent motion described in Figure 2 taken from Chae (2023a). The one-particle equivalent motion can be conveniently used for the Newtonian simulation of the relative motion between the two stars.  \nThe Newtonian mock PMs are obtained as follows.  \n- 1. For each system, eccentricity ( e ) is assigned using the observational input as described in Section 2.4. The longitude of the periastron ( ϕ 0 ) is assigned from the range (0 , 2 π ). The phase angle ϕ is obtained by solving Equation (10) of Chae (2023a) for a time t randomly drawn from (0 , T ) where T is the period also determined from the same equation. The inclination angle i is drawn from (0 , π/ 2) with a probability density function p ( i ) = sin i .\n- 2. The 3D separation r is obtained by s/ √ cos 2 ϕ +cos 2 i sin 2 ϕ . The semi-major axis follows from a = r (1 + e cos( ϕ -ϕ 0 )) / (1 -e 2 ).\n- 3. The magnitude of the relative 3D velocity is given by  \nv ( r ) = √ GM tot r ( 2 -r a ) . (13)  \n𝑦𝑦  \nThe sky-projected relative velocity components are then given by  \nv p,x ' = v ( r ) cos ψ, v p,y ' = v ( r ) cos i sin ψ, (14)  \nwith  \nψ = tan -1 ( -cos ϕ + e cos ϕ 0 sin ϕ + e sin ϕ 0 ) + π. (15)  \n4. The mock PM components are given by  \nµ ∗ α,A = µ ∗ α,M ± ( M B /M tot ) v p,x ' /d A , µ ∗ α,B = µ ∗ α,M ∓ ( M A /M tot ) v p,x ' /d B , µ δ,A = µ δ,M ± ( M B /M tot ) v p,y ' /d A , µ δ,B = µ δ,M ∓ ( M A /M tot ) v p,y ' /d B , (16)  \nwhere d A and d B are distances to the components, which can be taken to be the same as the errorweighted mean of the observed distances, 8 and µ ∗ α,M and µ δ,M are physically irrelevant constants chosen to be the error-weighted means of the observed PM components.  \nOnce the mock PMs (Equation (16)) are obtained, v p and then ˜ v can be calculated in the same manner as used for the observed PMs.", '2.3.3. Newtonian Proper Motions: with hidden companions': "If the binary system has a hidden inner binary for one component (or two inner binaries for both components), the apparent motion of a photocenter with respect to the barycenter in each inner binary is calculated and added to the outer velocity components. The details can be found in step #7 of Section 3.4 of Chae (2023a). Here I give a brief description.  \nThe same companion masses randomly assigned for the observed binaries (as described in Section 2.3.1) are used for the simulation. Thus, mock and real systems share the same masses and the same projected separation. The semi-major axis of the inner orbit a in is sampled from 0.01 au to ( d /pc) au with a steep powerlaw distribution of probability density p ∝ ( a in /a out ) -0 . 8 (Figure 14 of Chae 2023a; a nearly uniform probability density in log space) that was derived from near-by surveys (e.g. Tokovinin 2008). The upper limit represents 1 '' (Lindegren et al. 2021) within which unresolved component may be hidden. Other parameters such as i , e ,  \nϕ 0 and ϕ are assigned in a similar way as for, but independently of, the main outer orbit (see Chae 2023a for the details).  \nOnce the sky-projected relative velocities for the inner orbit are obtained, they are multiplied by a dimensionless factor η phot that mainly accounts for the distance of the photocenter from the barycenter normalized by the 3D separation between the host and the companion. The value of η phot is given by  \nη phot = { 0 ( P in < 3 yr) , M h M c ( M α -1 h -M α -1 c ) ( M h + M c )( M α h + M α c ) ( θ in < 1 '' ) , (17)  \nwhere P in is the period, θ in is the projected angular separation, and M h and M c are the masses of the host and the companion, and luminosity ∝ (mass) α with α = 3 . 5. As indicated in Equation (17), if the period of the inner orbit is shorter than the time span of 3 years of Gaia DR3 observations, the motion does not make a fixed contribution to PM components but may have contributed to the reported uncertainties of PM components and parallax. 9 Finally, the PM components from the inner orbit follow from the adjusted velocities, and are component-by-component added to the PM components from the outer orbit.", '2.4. Observational Inputs': "The MC procedure requires various observational inputs including the magnitude-mass relation, orbital eccentricities, and the statistical properties of hidden close companions (see Chae (2023a) for the details). As noted in Section 2.1, it is sufficient to consider the standard magnitude-mass relation for the considered magnitude range. Chae (2023a) considered several variations in the statistical properties of hidden close companions and found no significant systematic concern as long as f multi is self-calibrated for each sample.  \nEccentricities can have relatively large impact on wide binary tests of gravity because the Newtonian prediction on the relative velocity and the 'kinematic acceleration' depend significantly on the orbital eccentricity for the observed mass and sky-projected separation of the binary system. Eccentricities for all binaries of the ElBadry et al. (2021) sample were investigated by Hwang et al. (2022) by analyzing the angle between the relative displacement vector and the relative velocity vector on  \nthe sky. They provide binary-specific eccentricity ranges for all binaries.  \nAs in Chae (2023a, 2024a) the Hwang et al. (2022) eccentricity ranges will be used as the standard input. From the ranges, eccentricities are drawn using a truncated asymmetric Gaussian function as described in Chae (2023a, 2024a). However, I update the code (Chae 2023b) so that the truncated part is handled more properly. In the previous version, when a sampled e falls outside the range (0 . 001 , 0 . 999), its value is replaced by the boundary value. This procedure has the caveat that some occurrences can accumulate at the boundaries. To remedy this caveat, in the updated version e is resampled from an initial pool of e values that are within the range (0 . 001 , 0 . 999). I note that this modification has a small effect. Unless specified otherwise, eccentricities from the updated code are used in this work as the standard/default choice.  \nWhile the Hwang et al. (2022) individual ranges are most informative at present, the ranges are broad and may not be so accurate in some cases. Thus, I also consider sampling eccentricities from a power-law probability distribution of the form p ( e ; α ) = (1 + α ) e α . I consider two cases: (1) an empirically determined α varying as a function of s as given by Equation (18) of Chae (2024a), (2) the 'thermal' eccentricity distribution with α = 1. In the former, α is close to thermal at relatively smaller separation (0 . 2 < s ≲ 1 kau) but becomes superthermal at larger separations. Eccentricities sampled from this systematically varying power-law distribution will be referred to as 'statistical eccentricities' while those sampled from the α = 1 distribution 'thermal eccentricities'. Note that based on the Hwang et al. (2022) measurements, thermal eccentricities are systematically biased.", '2.5. Profiles of Normalized Projected Velocity and Mass with a Normalized Projected Radius': 'As described in Section 2.2, the relation between ˜ v and s/r M can be used as a proxy for the relation between g/g N and g N without a deprojection of the observed 2D quantities to the 3D space. Before considering a statistical method of testing gravity with the relation between ˜ v and s/r M , in this subsection I investigate the observed properties of ˜ v and associated quantities with respect to s/r M . Here four quantities are involved: v p (Equation (6)), v c ( s ) (Equation (10)), s , and r M (Equation (2)). But, both v c ( s ) and r M depend on M tot , and thus the key quantities are v p , s , and M tot .  \nThe projected separation s has negligible uncertainty. Also, because the required fractional uncertainty for the scalar relative PM is < 0 . 01 (since the four PM compo-  \nnents have fractional errors < 0 . 005) from the sample selection, the measurement uncertainty of v p 10 is practically negligible compared with the uncertainty arising from the uncertainties of hidden companions. In other words, the mass uncertainty associated with hidden components dominates the uncertainties of both ˜ v and s/r M . In obtaining and using the relation between ˜ v and s/r M , i.e. the profile ˜ v ( s/r M ), the uncertainties need to be taken into account. When hierarchical systems are present as in the Chae (2023a) sample and the new sample, the nominal analyses will be carried out including the uncertainties associated with hierarchical systems but not the PM measurement uncertainties. However, for the pure binary sample the PM measurement uncertainties are always included. In Appendix B, I explicitly show that the PM measurement uncertainties can only have a minor effect when hierarchical systems are present.  \nI use the MC method described in Section 2.3 that is similar to those used in Chae (2023a, 2024a) to take into account the uncertainties. The mass of additional component(s) is assigned for a given value of f multi with the method presented in Section 3.2 of Chae (2023a). One MC draws one set of values from the respective probability distributions of ˜ v and s/r M . By drawing many MC sets, the uncertainty of the distribution of (˜ v, s/r M ) will be derived and used in the statistical analysis to be described in the next subsection.  \nFigure 3 shows distributions of (˜ v , s/r M ) for the four samples summarized in Table 1 from one MC set for each sample. Here I use f multi = 0 . 43 for the Chae (2023a) sample, and f multi = 0 . 18 for the new sample as determined from this work (see Appendix A). I consider six bins uniformly spaced in log 10 ( s/r M ) within the specified range bounded by the vertical dotted lines and calculate the medians of ˜ v in the bins. As will be shown below, the lower and upper limits are set to exclude data that give biased distributions of M tot at s/r M .  \nFigure 4 shows one MC distribution for a virtual Newtonian sample as a Newtonian counterpart of the Chae (2023a) sample. The virtual Newtonian sample is obtained by replacing the observed PMs with simulated PMs as described in Section 2.3.  \nThe derived binned medians ⟨ ˜ v ⟩ for the Gaia samples increase slowly with s/r M and an elevation is apparent for the bins of s/r M ≳ 0 . 5 compared with the lower s/r M bins. On the contrary, the virtual Newtonian sample is  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 3. Each panel shows a distribution of binaries in the plane spanned by logarithmic values of two dimensionless quantities s/r M (sky-projected separation over the MOND radius) and ˜ v (Equation (1)) from one MC realization of binary total masses including probabilistically assigned masses of hidden close binaries. Four samples summarized in Table 1 are shown. The bins and the range used for analyses in log 10 ( s/r M ) are indicated by vertical lines. For each sample the overall median value of ˜ v is indicated by the horizontal dashed line and the binned medians are indicated by big dots and the thick solid line. I emphasize that the scattered points should be regarded as a probability distribution from all binaries with all uncertainties and possible ranges of parameters taken into account in a natural way through the MC. See the text for further.  \n<!-- image -->  \nFigure 4. This figure is the same as the upper left panel of Figure 3 but for a virtual Newtonian sample as a counterpart of the Chae (2023a) Gaia sample. Note that the binned medians of ˜ v are nearly flat in the virtual sample.  \n<!-- image -->  \nnearly flat from the lowest to the highest bin of s/r M . Consequently, the median ⟨ ˜ v ⟩ = 0 . 55 of the virtual New-  \nn sample for the entire range is lower than the value 0 . 58 of the counterpart Gaia sample.  \nThe overall medians of ˜ v for the new and the pure binary samples differ somewhat from that for the Chae (2023a) sample. This can be attributed to the fact that the three samples have different f multi and ⟨ M tot ⟩ (as shown below). The sample with the limited range of 2 < s < 30 kau exhibits a minor variation with s/r M . The overall median ⟨ ˜ v ⟩ for this sample is clearly higher than those for the other samples. This is consistent with the fact that the binned median increases with s/r M .  \nBecause each MC draw generates values of ˜ v and s/r M from their respective probability distributions taking fully into account the uncertainties and the values from all binaries are stacked on the plane, there is no need in the present approach, as in the acceleration-plane approach of Chae (2023a), to remove binaries based on a fractional uncertainty of ˜ v as suggested by Banik et al. (2024). In other words, if a binary system has a relatively larger uncertainty (i.e. a broader probability dis-  \nibution) of ˜ v , its MC-drawn values are scattered more broadly than others and thus naturally given a lower weight at fixed s/r M .  \nIn the approach taken by Banik et al. (2024), they applied a hard cut on the fractional uncertainty of ˜ v and then treated all values equally ignoring their individual uncertainties. As Chae (2024a) has already demonstrated, applying such a cut on a binary sample does not change the gravitational anomaly. In other words, the MC-based acceleration-plane analysis is largely immune to any sampling (excluding of course erroneous data) because all uncertainties are naturally taken into account. The same is true for the MC-based analysis of the profile ˜ v ( s/r M ) as will be demonstrated in this work.  \nHowever, I consider a kinematic cut based on an uncertainty of ˜ v for the purpose of illustrating its effects in an analysis directly using ˜ v . For this purpose, a formal uncertainty is estimated as  \nσ ˜ v = ˜ v √ ( σ v p v p ) 2 + ( σ v c v c ) 2 , (18)  \nwhere σ v c follows from the uncertainty of M tot . I consider the kinematic cut suggested by Banik et al. (2024)  \nσ ˜ v < 0 . 1 max(1 , ˜ v/ 2) . (19)  \nI take σ v c /v c = 0 . 05 assuming a uniform uncertainty of 10% in M tot . I note that this formal uncertainty is of limited value because it is not known which systems are hierarchical and thus system-specific realistic uncertainties cannot be determined. Banik et al. (2024) considered an MC error propagation to estimate the uncertainties of ˜ v in their sample. However, they did not take into account the uncertainties of hierarchical systems either. Thus, the uncertainty given by Equation (18) or that estimated by Banik et al. (2024) cannot be considered accurate, and it is necessary to consider an MC method to take into account hierarchical systems as done in this work.  \nIn Figure 5, binaries from the Chae (2023a) sample are shown with different colors split by the kinematic cut of Equation (19). Figure 5 also shows the distributions of v p , v c , and M tot with respect to s/r M . Figure 5 shows that when all binaries are considered (black solid lines), ⟨ ˜ v ⟩ (i.e. the binned median of ˜ v ) increases with s/r M as a consequence of ⟨ v p ⟩ deviating from the Keplerian scaling ∝ s -1 / 2 while ⟨ v c ⟩ follows the Keplerian scaling as ⟨ M tot ⟩ is flat. However, for the binaries satisfying the kinematic cut, ⟨ ˜ v ⟩ (blue dashed line) does not vary with s/r M as was first noted by Banik et al. (2024). For those binaries ⟨ v p ⟩ follows the Keplerian scaling. Thus, at first glance they appear to obey Newtonian dynamics. This  \nwas indeed the claim by Banik et al. (2024). This would be the case if ⟨ M tot ⟩ did not vary and thus ⟨ v c ⟩ followed the Keplerian scaling.  \nHowever, as the bottom panels of Figure 5 show, the binary sample with the kinematic cut is biased so that ⟨ M tot ⟩ increases with s/r M and thus ⟨ v c ⟩ deviates from the Keplerian scaling. This means that the Newtonian prediction of ⟨ ˜ v ⟩ for this sample is different from the general sample for which ⟨ M tot ⟩ does not vary. This indicates that the apparently Newtonian flat behavior of ⟨ ˜ v ⟩ with the kinematic cut is driven by the variation of ⟨ M tot ⟩ with s/r M . This will be addressed quantitatively in this work.  \nFurther insights may be gained from considering the virtual Newtonian sample. Figure 6 shows the distributions in the virtual Newtonian sample. For this sample, when all binaries are considered (black solid lines), ⟨ ˜ v ⟩ does not vary with s/r M as a consequence of both ⟨ v p ⟩ and ⟨ v c ⟩ following the Keplerian scaling. However, for the binaries satisfying the kinematic cut, ⟨ ˜ v ⟩ now decreases as a consequence of ⟨ v c ⟩ deviating from the Keplerian scaling and ⟨ v p ⟩ (nearly) following the Keplerian scaling. What is striking is that the effects of the kinematic cut on the scalings of the four quantities are the same for the Gaia and the virtual Newtonian samples. In particular, due to the variation of the median mass with s/r M when the kinematic cut is imposed, ⟨ ˜ v ⟩ is affected in the same manner by the cut however ⟨ ˜ v ⟩ scales with s/r M without the cut.  \nIn the Chae (2023a) sample, the fraction of binaries not meeting the kinematic cut is ≈ 18%. In the new sample this fraction decreases dramatically to ≈ 3%. Thus, the effects of the kinematic cut are minimal. In the pure binary sample, the fraction is just ≈ 1% and so its effects are negligible.  \nFigure 7 shows the distributions for the Chae (2023a) limited sample of 5635 binaries with the narrow range 2 < s < 30 kau. For this sample, variations of the quantities with s/r M are weak or absent. Since just the twice large dynamics range (in log scale) reveals a clear variation of ⟨ ˜ v ⟩ with s/r M (Figure 5), the weak/null variation in this range (that was used by Banik et al. (2024)) is clearly misleading. It could easily be interpreted as an indication for Newtonian gravity because the flatness of the ⟨ ˜ v ⟩ profile is its generic prediction (Figure 6).  \nIn principle, it would still be possible to test gravity with the flat ⟨ ˜ v ⟩ value from a limited dynamic range sample. However, such a test requires a precise value of ⟨ M tot ⟩ that in turn depends on f multi . This also means that it may well be possible to obtain a wrong gravity with a wrong value of f multi if it is not properly determined for the sample under consideration.  \nFigure 5. This figure shows the scaling of four quantities with s/r M from one MC result for the Chae (2023a) sample. They are ˜ v (that was already shown in Figure 3), v p (sky-projected relative velocity between the pair of stars), M tot (total mass of the binary system including probabilistically assigned mass of any hidden close binary), and v c (the Newtonian circular velocity due to M tot at the sky-projected separation s ). For the bins defined in the upper left panel of Figure 3, binned medians are shown for all (black solid line) and those (red and blue dashed lines) split by a kinematic constraint suggested by Banik et al. (2024). Once about 18% of binaries (small red dots) are removed by the kinematic cut, the binned medians of ˜ v and v p (blue lines) for the remaining binaries follow the Newtonian predictions. However, as the bottom panels show, binned masses for these binaries rise systematically with s/r M and thus v c systematically deviates from the Newtonian prediction showing that the subsamples are biased in the distribution of masses. See the text for further.  \n<!-- image -->', '2.6. Statistical Analysis': 'The distribution of binaries in the plane spanned by two dimensionless quantities ˜ v and s/r M for a sample is used to test gravity. For the range of s/r M within which M tot has an unbiased distribution (Figure 3), the bins of s/r M are defined as shown in Section 2.5. MC sets of (˜ v , s/r M ) are generated taking into account the uncertainties of M tot with a probability of f multi hosting undetected close companions (as already mentioned in Section 2.5 and described in Section 3.2 of Chae (2023a)). For each MC set, binned medians of ˜ v are obtained.  \nAlong with each MC set for the given sample of binaries, a corresponding MC set is derived from a mock Newtonian sample with simulated PMs replacing the given/measured PMs. For each pair of MC sets, medians ⟨ ˜ v ⟩ obs ,i and ⟨ ˜ v ⟩ newt ,i are obtained for the i -th bin, where i = 1 , · · · , N bin . From the medians, values of the velocity ratio ⟨ ˜ v ⟩ obs ,i / ⟨ ˜ v ⟩ newt ,i are obtained. Since the kinematic acceleration g ∝ v 2 and v p ∝ v statistically  \n(averaged over all possible orientations), this velocity ratio probes the ratio ( g obs /g pred ) 1 / 2 where g obs and g pred are the observed and Newtonian-predicted values of the kinematic acceleration.  \nN draw MC sets are generated and from them I obtain a distribution of a parameter defined by  \nΓ ≡ log 10 γ ˜ v ≡ log 10 ( ⟨ ˜ v ⟩ obs ⟨ ˜ v ⟩ newt ) , (20)  \nwhere γ ˜ v represents the velocity boost (or anomaly) factor and Γ is its logarithmic value. I have (Γ i ) j for the i -th bin ( i = 1 , · · · , N bin ) and the j -th draw ( j = 1 , · · · , N draw ). Since the MC distribution of (Γ i ) j for each i follows approximately a Gaussian distribution, I use the conventional χ 2 statistics (e.g. Wall & Jenkins 2012) for gravity model testing. For this I define a reduced χ 2  \nχ 2 ν ≡ 1 ν N bin ∑ i =1 ( µ Γ i -log 10 γ model ˜ v i ) 2 ( σ Γ i ) 2 +( σ model i ) 2 , (21)', 'Virtual Newtonian sample 20249) (Nbinary': 'Figure 6. Same as Figure 5 but for the virtual Newtonian sample as a counterpart of the Chae (2023a) Gaia sample. Note that unlike Figure 5, for the virtual binaries satisfying Newtonian gravity the binned medians of ˜ v (black solid line) for all binaries follow the flat Newtonian line while the blue dashed line deviates now downward from the Newtonian line. The behaviors of M tot and v c are the same as in Figure 5 because virtual binaries are statistically equivalent to true binaries except that Newtonian gravity is assumed for the virtual sample.  \n<!-- image -->  \nwhere ν is the degree of freedom, µ Γ i and σ Γ i are the mean and standard deviation of (Γ i ) j for j = 1 , · · · , N draw . Because each (Γ i ) j is derived in a fully Monte Carlo way from the Gaia observations along with possible ranges of all involved parameters, µ Γ i and σ Γ i are the measurements for the bins. In Equation (21), γ model ˜ v i is the theoretical velocity boost factor for the i -th bin predicted by a given model, and σ model i is the uncertainty of the model prediction due to the uncertainties of the external field effect (EFE; e.g., Chae et al. 2020, 2021) in a model having it.  \nIn Equation (21), the degree of freedom ν is equal to the number of bins N bins if f multi is fixed as in the pure binary sample. Otherwise, I have ν = N bins -1 because f multi is fitted by the high acceleration (or low s/r M ) data where both standard and nonstandard gravity theories converge.  \nBy its definition (Equation (20)), γ model ˜ v i = 1 and σ model i = 0 for the standard gravity. For non-standard gravity I consider two nonrelativistic MOND (Milgrom 1983) models, AQUAL (Bekenstein & Milgrom 1984) and QUMOND (Milgrom 2010) that permit numerical solutions of the gravity boost factor taking into account  \nthe EFE of the Milky Way. The difference between them is not large but a test with galactic rotation curves marginally prefers AQUAL (Chae et al. 2022). Thus, I will test specifically the AQUAL model.  \nFor the AQUAL model, numerical solutions for the elliptical motions of binaries under the external field of the Milky Way are not available at present. However, numerical solutions for circular orbits of a test particle have been studied extensively by Chae & Milgrom (2022) with various mass models of the gravitational source from point-like objects to realistic disks. Chae & Milgrom (2022) report that AQUAL numerical solutions predict the following fitting function for the ratio of Milgromian gravity to Newtonian gravity  \ng g N = ν ( y β ) [ 1 + { tanh ( βg N , ext /a 0 g N /a 0 )} γ ˆ ν ( y β ) 3 ] , (22)  \nwhere the definition y β ≡ √ ( g N /a 0 ) 2 +( βg N , ext /a 0 ) 2 is used with a Newtonian external field g N , ext , ν ( y ) is the MOND interpolating function, and ˆ ν ( y ) ≡ d ln ν ( y ) /d ln y . Chae & Milgrom (2022) use the simple function (Famaey & Binney 2005) ν ( y ) =  \nFigure 7. Similar to Figure 5 but for the Chae (2023a) limited sample satisfying s > 2 kau. Because of the narrow range of s/r M after excluding binaries suffering from the edge effects, there is no obvious scaling of ˜ v , and its overall median is higher than that shown in Figure 5 for the sample including the s < 2 kau binaries.  \n<!-- image -->  \n0 . 5 + √ 0 . 25 + y -1 (for the acceleration range g N ≲ 10 -7 m s -2 ) for which the relation between the empirical external field g ext and the Newtonian external field is given by x 2 / (1 + x ) = g N , ext /a 0 where x ≡ g ext /a 0 .  \nThe AQUAL model described by Equation (22) is implemented as follows. I take β = 1 . 1 and γ = 1 . 2 that were estimated using orbits with an average inclination 60 · of the external field with respect to the orbital axis. 11 For the eternal field I take g ext /a 0 = 1 . 90 ± 0 . 19 (Chae 2023a, 2024a) with a 0 = 1 . 2 × 10 -10 m s -2 . The 10% uncertainty reasonably covers its current uncertainty.  \nThe AQUAL model can be related to the velocity boost factor γ model ˜ v defined in Equation (20) as follows. From g = v 2 /r and g pred = g N = v 2 c /r for circular orbits, the ratio g/g N at g N /a 0 is equal to ( v/v c ) 2 at r/r M . Now averaging over random orientations and phases, I get ˜ v ≈ v/v c at r ≈ 1 . 2 s . Thus, I have the approximate  \nrelation  \nγ model ˜ v ≈ √ g g N with s r M ≈ 1 1 . 2 √ a 0 g N , (23)  \nwhere g/g N is given by Equation (22). Thus, Equation (23) gives an AQUAL prediction as a function of s/r M to be tested through Equation (21).', '3.1. Test with a virtual Newtonian sample': 'Before presenting the results for the Gaia samples I here test the statistical analysis with a virtual Newtonian sample. This sample is similar to the Chae (2023a) sample except that the measured PMs were replaced with simulated PMs as described in Section 2.3. This sample was produced with a fixed f multi = 0 . 43 that is the value required for the Chae (2023a) sample based on the most up-to-date acceleration-plane analysis (Appendix A).  \nFigure 8 shows the results for the virtual Newtonian sample as the input data. The upper panels show the distributions of the binned medians of ˜ v from N draw = 200 MC sampling for the input data ( ⟨ ˜ v ⟩ obs ) and the corresponding Newtonian simulated data ( ⟨ ˜ v ⟩ newt ). The bottom panels show the distribution  \nFigure 8. The upper panels show the distributions of binned medians of ˜ v from 200 MC results for the input sample and the corresponding virtual Newtonian sample. One set of binned medians along with all values of ˜ v from one MC are shown in Figure 3. For this figure, the input sample itself is a virtual Newtonian sample. The bottom panels show the inferred values of Γ( ≡ log 10 ( ⟨ ˜ v ⟩ obs / ⟨ ˜ v ⟩ newt ), Equation (20)) where ⟨ ˜ v ⟩ obs and ⟨ ˜ v ⟩ newt are the values from the input and virtual Newtonian samples. Newtonian and AQUAL theories are tested through the reduced χ 2 statistics (Equation (21)) with the measured values of Γ. The left column is for the full sample while the right column is for the subsample satisfying the kinematic cut. The Newtonian prediction is indicated by the horizontal black dashed line while the AQUAL prediction is indicated by the magenta curve with a band representing the uncertainty of the external field effect. As expected, the Newtonian model is well acceptable while the AQUAL model is ruled out by the virtual Newtonian sample. Note also that the kinematic cut does not change the χ 2 testing significantly.  \n<!-- image -->  \nof Γ ( ≡ log 10 ( ⟨ ˜ v ⟩ obs / ⟨ ˜ v ⟩ newt ), Equation (20)). The left column shows the results for the full sample for which mass is not biased (see Figure 6). The right column with the kinematic cut of Equation (19) is considered to investigate its effects.  \nIn the bottom panels of Figure 8, the Newtonian and AQUAL predictions are compared with the derived distributions of Γ. As expected, the distributions of Γ agree well with the Newtonian prediction of Γ = 0 whether the kinematic constraint is imposed or not. The derived values of χ 2 ν = 0 . 6 and 1 . 0 for the Newtonian model are well acceptable by the reduced χ 2 ν static. In contrast, the AQUAL model has χ 2 ν = 15 . 5 and 11 . 9 for which the complementary (or survival) probability is  \nP c ≡ 1 -P ( < χ 2 ν ) ≲ 2 × 10 -13 . 12 Thus, if the standard gravity holds for wide binaries, the AQUAL model is expected to be completely ruled out with a Gaussian equivalent confidence well above 5 σ .', '3.2. Results for Gaia samples with a broad dynamic range': "In this subsection I present the results for three Gaia samples (Table 1) with the broad dynamic range 0 . 2 < s < 30 kau. For the Chae (2024a,b) pure binary sample, f multi = 0 is fixed, so the analysis is the same as that for the virtual Newtonian sample presented in Section 3.1. For the Chae (2023a) and the new samples, the best-fit value of f multi can be determined with either lows/r M  \nbins or high-acceleration data through the accelerationplane analysis of Chae (2023a). I find that both methods give consistent results except when mass-biased samples are used as in the cases with the Banik et al. (2024) kinematic cut. Even in the latter cases, the differences are moderate and the results are not significantly affected whichever method is used. I will use the f multi values determined through the acceleration-plane method as presented in Appendix A because it is more proper in principle.  \nFigure 9 shows the distributions of ⟨ ˜ v ⟩ obs , ⟨ ˜ v ⟩ newt , and Γ for the Chae (2023a) sample. The best-fit value of f multi for this sample is 0 . 43. The distributions for the mass-unbiased sample (see Figure 5) shown in the left column exhibit almost the exact opposite properties of those for the virtual Newtonian sample shown in the left column of Figure 8. The AQUAL model fits the binned data adequately with χ 2 ν = 2 . 7 while the Newtonian model fails completely with χ 2 ν = 20 . 3 for which the survival probability is P c ≈ 2 . 6 × 10 -20 with a Gaussian equivalent significance of 9 . 2 σ that is similar to 10 σ obtained by Chae (2023a) with a somewhat larger sample.  \nIf some binaries are removed by the kinematic cut of Equation (19), two binned velocities, ⟨ ˜ v ⟩ obs and ⟨ ˜ v ⟩ newt , are lowered progressively with s/r M as already noted in Figure 5. However, the AQUAL model is still not ruled out with χ 2 ν = 3 . 0 ( P c ≈ 0 . 01) by the conventional standards while the Newtonian model is ruled out with χ 2 ν = 7 . 7 ( P c ≈ 3 . 0 × 10 -7 ).  \nMore importantly, I reemphasize that this cut is not necessary in the current MC approach (as in the acceleration-plane MC approach of Chae (2023a)) because all uncertainties and likely ranges are taken into account already. As the right column of Figure 9 shows, the main role of the cut is to systematically twist the normalized velocity profiles in an unnatural way. This is the consequence of the cut-introduced bias in the mass profile. For this reason, the result for the cut-introduced biased subsample should not be regarded as a representative or preferred result. It is presented here merely to clarify the effects of such a cut. For the rest of this work, I will not consider the cut-introduced subsample of the Chae (2023a) sample.  \nIf a sample satisfying the cut such as Equation (19) is desired, it is necessary to apply, in an overall way, stricter selection criteria that naturally select an unbiased sample automatically satisfying the cut while not biasing the mass profile. It turns out that the new sample largely meets the cut while the pure binary sample meets the cut almost completely as already noted in Chae (2024a).  \nFigure 10 shows the results for the new sample. This sample has been selected with different selection criteria and has a significantly lower fraction of hierarchical systems with f multi = 0 . 18. About 97% of this sample already satisfy the kinematic cut. Nevertheless, the observed normalized velocity profile from the full sample (the left column) is statistically equivalent to that from the Chae (2023a) full sample. The subsample satisfying the kinematic cut gives a very similar profile. The AQUAL model is well acceptable by both profiles with χ 2 ν = 1 . 6 and 2 . 4 while the Newtonian model is ruled out with χ 2 ν = 9 . 4 and 9 . 7 with a Gaussian equivalent significance of ≈ 5 . 8 σ .  \nThe above results are based on the individual eccentricity ranges as described in Section 2.4. Now I present results based on less informative eccentricities from statistical distributions to investigate the effects of possible variations of the measured eccentricities. As described in Section 2.4, I consider 'statistical eccentricities' from an observationally-inferred varying power-law distribution and 'thermal eccentricities' from a fixed power-law distribution.  \nFigure 11 and Figure 12 respectively show the results for the Chae (2023a) sample and the new sample based on statistical and thermal eccentricities. With the statistical/thermal eccentricities, the determined values of f multi are somewhat lower. For both samples, the inferred Γ values with statistical eccentricities agree well (or adequately) with the AQUAL model and rule out the Newtonian model. Even with thermal eccentricities, the AQUAL model is clearly preferred in a relative sense.  \nFinally, in Figure 13 I present results for the (Chae 2024a,b) pure binary sample. The results are well consistent with those for the Chae (2023) sample and the new sample. With the individual eccentricities the AQUAL model is well acceptable by the inferred Γ values, but the Newtonian model fails with a high statistical significance ( P c ≈ 8 × 10 -6 for χ 2 ν = 5 . 6) despite the relatively small sample size. Also, with statistical or thermal eccentricities, the AQUAL model is preferred over the Newtonian model.  \nFrom the above results the logarithmic velocity boost factor for s/r M ≳ 0 . 7 is measured to be  \nΓ =      0 . 059 ± 0 . 009 Chae (2023a) 0 . 068 ± 0 . 015 new (individual e ) 0 . 092 ± 0 . 023 pure binary (24)  \nFigure 9. The same as Figure 8 but for a real sample of binaries, i.e. the Chae (2023a) sample. The χ 2 testing outcomes are almost excatly opposite to those for the virtual Newtonian sample shown in Figure 8. In particular, even with the kinematic cut the χ 2 testing conclusion is unaffected: the AQUAL model is adequate while the Newtonian model is ruled out. However, I emphasize that in the present MC approach the kinematic cut is not necessary or motivated as explained in detail in the text.  \n<!-- image -->  \nwith the individual eccentricities and  \nΓ =      0 . 053 ± 0 . 008 Chae (2023a) 0 . 058 ± 0 . 014 new (statistical e ) 0 . 077 ± 0 . 022 pure binary  \n(25)  \nwith statistical eccentricities. Here the quoted uncertainties include the uncertainties of f multi . These results alone rule out the Newtonian model with a significance ranging from 3 . 9 -6 . 6 σ . I do not quote values from the results with thermal eccentricities because they are not based on observational inference but considered merely as a possible systematic variation.  \nTo summarize, three samples with various sample sizes (from ≈ 3500 to 20000), different selection criteria, and different fractions of hierarchical systems (from f multi ≈ 0 to ≈ 0 . 4) are all well consistent with the AQUAL model and rule out the Newtonian model with a varying degree of confidence depending on the sample size. Table 2 gives a summary of the reduced χ 2 values, the survival probabilities, and the logarithmic values of the velocity boost Γ (Equation (20)). From the inferred values of Γ, the gravitational acceleration boost factor or gravitational anomaly is found to be 10 2Γ ≈ 1 . 4 ± 0 . 1  \nwell consistent with recent results (Chae 2023a, 2024a; Hernandez 2023; Hernandez et al. 2024).", '3.3. Test with a Gaia sample with a narrow dynamic range': "Here I carry out a quantitative analysis of the subsample with the limited dynamic range 2 < s < 30 kau from the Chae (2023a) sample. This is the range considered by Banik et al. (2024) while Pittordis & Sutherland (2023) considered an even narrower range of 5 < s < 20 kau. Comparison of Figure 7 with Figure 5 already indicates that the subsample with the limited dynamic range has a nearly flat ˜ v profile with s/r M and cannot easily distinguish between AQUAL and Newton because AQUAL predicts a nearly flat profile for s/r M ≳ 0 . 7. Then, the different medians (not the profiles) of the two models may be used to distinguish the models. For this test to work, f multi must be known accurately for the test sample because of the inherent degeneracy between total mass and gravity in binaries. However, for this narrow range the acceleration-plane analysis cannot be used to determine f multi because the sample lacks the Newtonian regime ( ≳ 10 -8 m s -2 ) data.  \nFigure 10. The same as Figure 9 but for the new sample defined in this work. The χ 2 testing outcomes are qualitatively the same as those for the Chae (2023a) sample. In particular, the AQUAL model is well acceptable while the Newtonian model is ruled out at > 5 σ Gaussian equivalent significance. The kinematic cut has negligible effects because the sample naturally satisfies the cut thanks to the overall higher quality requirements without biasing the mass distribution.  \n<!-- image -->  \nTable 2. The inferred velocity boost factors and the χ 2 test Results  \nNote. Parameter Γ is the logarithmic value of the velocity boot factor defined by Equation (20). Quantity 10 2Γ corresponds to the gravitational acceleration boost factor or gravitational anomaly.  \nFigure 11. Each column is the same as the left column of Figure 9 for the Chae (2023a) sample except that statistical or thermal eccentricity distributions are used instead of the individual eccentricities. The result with statistical eccentricities is well consistent with that with the individual eccentricities. Because thermal eccentricities are not based on pertinent observations for the binary sample, the result with them should not be used to test gravity. It is given here merely as a systematically different one. Not surprisingly, the Γ values in the bottom panel of the right column result in an acceptably large value of χ 2 ν = 4 . 7 for the AQUAL model although the Newtonian model has a much larger value of 10 . 0.  \n<!-- image -->  \nFigure 12. The same as Figure 11 but for the new sample. The results are qualitatively the same as those shown in Figure 11.  \n<!-- image -->  \nFigure 13. This figure shows the results for the pure binary sample with three choices of eccentricities. This highest quality sample selected with very stringent selection criteria (Chae 2024a,b) automatically satisfies the kinematic cut almost completely. The results are well consistent with those with the 'impure' samples shown in Figure 9, 11, 10, and 12 showing that varying degree of data quality requirements does not affect the result of gravity test as long as an MC taking fully into account possible parameter ranges and uncertainties is used.  \n<!-- image -->  \nIf f = 0 . 43 ± 0 . 05 determined from the Chae (2023a) sample with 0 . 2 < s < 30 kau is used, the results for the subsample will of course be consistent with those for the full sample. It is then interesting to ask whether the limited range sample can prefer Newton as claimed by Banik et al. (2024). Thus, I consider a fixed value of f multi = 0 . 65 that is a typical value suggested by Banik et al. (2024). Note, however, that the Banik et al. (2024) sample and the present sample cannot be considered statistically equivalent because the Banik et al. (2024) sample includes kinematic contaminants such as chance alignments (and thus their analysis had to introduce extra free parameter(s).)  \nFigure 14 shows the ˜ v profiles, the binned measurements of Γ, and the χ 2 ν tests of the AQUAL and the Newton models. Remarkably, with f multi = 0 . 65 this narrow range data are perfectly consistent with the Newtonian prediction with χ 2 ν = 1 . 5 or 0 . 8. In contrast, the AQUAL model is ruled out with χ 2 ν = 15 . 0 or 9 . 0 (up to a Gaussian equivalent significance of ≈ 8 σ ). These results are fully consistent with those obtained by Banik et al. (2024) though formal statistical significance is not as high as their 19 σ . Thus, though with a different sample, this test with the same narrow dynamic range confirms the Banik et al. (2024) conclusion. Indeed, this conclusion is the consequence of the narrow dynamic range and the high occurrence rate of hierarchical systems.", '3.4. A Deeper Test with Distributions of ˜ v in Bins of s/r M': "So far all tests have been done with parameter Γ (Equation (20)), i.e. the median of log 10 (˜ v obs / ˜ v newt ), in bins of s/r M . The median is relatively robust against any potential outliers that may arise from some unknown systematic in ˜ v obs and/or ˜ v newt . The Newtonian model has already failed the test with Γ while the MOND model represented by a generic numerical solution of the AQUAL field equation is acceptable. When a model fails a first-order/tier test, there is little point of considering a deeper or fuller test.  \nThus, a deeper test beyond Γ may be warranted only for MOND gravity based on the results of this work so far. However, I also consider testing Newton for the purpose of comparison. A popular test in the recent literature of wide binary gravity tests has been to compare the distribution/histogram of ˜ v obs (the observed value) with that of ˜ v mod (the prediction of the model under consideration) in some ranges of s (e.g. Pittordis & Sutherland 2019, 2023; Clarke 2020; Banik et al. 2024; Hernandez et al. 2024). However, because binary systems have different total masses, it is more accurate to consider bins of s/r M .  \nWhile it is straightforward to obtain ˜ v mod from physically possible orbits for Newtonian gravity as described in Section 2.3, it is not the case for Milgromian gravity, in particular for a fully non-linear field equation such as AQUAL. Here I take advantage of the characteristics of the obtained gravitational anomaly from the studies of this and previous works (Chae 2023a, 2024a; Hernandez et al. 2024), which show consistently that the observed gravity anomaly occurs in a rather narrow range 1 ≲ s ≲ 7 kau (mostly in 2 ≲ s ≲ 5 kau) for binaries with typical total masses M tot ≈ 1 . 0 -1 . 5 M ⊙ , and then becomes pseudo-Newtonian with a boosted effective gravitational constant G eff ≈ 1 . 4 G . In the context of MOND gravity, this pseudo-Newtonian behavior is the consequence of the specific gravitational environment thanks to a super-critical external field of the Milky Way.  \nIn this specific environment common 13 to all nearby binaries, I approximate Milgromian gravity by a model with varying Newton's constant depending on s/r M . The effective Newton's constant depends on s/r M through G eff = γ g ( g N ) G with g N ≈ a 0 (1 . 2 s/r M ) -2 (see Equation (23)), where the function γ g ( g N ) represents the right-hand side of Equation (22). This model becomes Newtonian as s/r M → 0 and pseudo-Newtonian with G eff = 1 . 37 for s/r M ≳ 1. In the transition range 0 . 1 ≲ s/r M ≲ 0 . 7, this model is less accurate and should be taken with a grain of salt.  \nI consider three representative bins of s/r M to test the distribution/histogram of ˜ v : 0 . 040 < s/r M < 0 . 055 (deep Newtonian), 0 . 3 < s/r M < 0 . 5 (transition), and 0 . 8 < s/r M < 1 . 6 (MOND). The first bin can be used to test the validity of the MC algorithm itself because gravity is known a priori for the bin. In other words, all gravity models (whether Newtonian or Milgromian) must agree with the Gaia data. Note that for this purpose I consider a narrow range of s/r M ≪ 1 with the upper limit 0.055 chosen empirically to minimize the possibility of including systems with r ≳ r M (where r is the 3D separation deprojected from s ) as a result of projection effect, and the lower limit 0.040 chosen to avoid the edge effect (see the range shown in Figure 5). Note also that for the last bin the projection effect is not a concern because all systems with s/r M ≳ 1 are guaranteed to be in the MOND regime since r ≥ s . The upper limit of 1 . 6 (corresponding to s ≈ 11 kau for a one solar mass binary) in the last bin is chosen as a compromise between number statistics and the width of  \nFigure 14. This figure shows gravity test results for the 'Chae (2023a) limited sample' the limited dynamic range s > 2 kau. Because it lacks the Newtonian regime binaries and thus f multi could not be determined, a fixed value of f multi = 0 . 65 (a value from Banik et al. (2024)) is used. For this test, statistical eccentricities are used considering that Banik et al. (2024) do not have individual eccentricities for their binaries. Indeed, for this high value of f multi the Newtonian model is acceptable while the AQUAL model is ruled out to agree with the Banik et al. (2024) conclusion. However, the input value of f multi = 0 . 65 is ruled out by the determined value of f multi = 0 . 36 ± 0 . 05 with statistical eccentricities for the sample of the full range 0 . 2 < s < 30 kau.  \n<!-- image -->  \nthe range. Fortunately, this limit also helps preventing any possible contamination of the kinematic data for widest binaries. (However, below I will also consider a full MOND range of 0 . 75 < s/r M < 2 . 50.) The middle bin of 0 . 3 < s/r M < 0 . 5 is taken from the NewtonMOND transition regime 0 . 1 ≲ s/r M ≲ 0 . 7.  \nFigure 15 shows the results of comparing the Newtonian and MOND predictions with the observed distribution of ˜ v in three samples that include different fractions of hierarchical systems in the range 0 ≤ f multi ≤ 0 . 43. For the deep Newtonian regime, the Newton/MOND predicted shape agrees nearly perfectly with the Gaia distribution for all samples regardless of f multi . It can be seen clearly that the distribution gets gradually broadened with a longer tail as f multi increases. As the dotted curves indicate, the broadening and the tail arise from the hierarchical systems. This gradual broadening appears to be accurately predicted by the MC modeling. I consider the K-S test (see, e.g., Wall & Jenkins 2012) to check the overall visual agreement quantitatively. The K-S test is performed with the Python module scipy.stats.ks 2samp between the model (Newton  \nor MOND) MC set and the observed MC set, which includes the random scatter due to hierarchical systems (see below for the results even including the small PM measurement scatters). The K-S p -value provides a survival probability for the null hypothesis that the two samples are drawn from the same underlying distribution. The high values (0 . 2 ≲ p -value ≲ 0 . 6) confirm the agreement between the observed data and the MCproduced model data. Also, the Newton/MOND medians indicated by downward-pointing arrows match well the observed medians. Thus, these results for the deep Newtonian regime robustly verify the validity of the MC algorithm including the kinematic effects of hidden close companions.  \nIn the transition regime (0 . 3 < s/r M < 0 . 5) shown in the middle column of Figure 15, the MOND prediction starts to deviate from the Newtonian prediction as the former is more broadened than the latter. Because the observed distribution is broadened with respect to the deep Newtonian case, the MOND model is preferred over the Newtonian model for all the samples. The MOND predicted medians agree well with the  \nFigure 15. This figure compares the Gaia observed distribution of ˜ v with the Newtonian and the MOND/pseudo-Newtonian predictions in three representative bins of s/r M for the three samples used in this work. The left column shows the results for the deep Newtonian regime. The shape of the distribution is gradually broadened from the statistically pure binary sample with f multi = 0 to the largest sample with f multi = 0 . 43 which is well reproduced by the MC algorithms of both models that are indistinguishable in this regime. As the gravity regime changes from the deep Newton to MOND at fixed f multi , the distribution of ˜ v is also broadened. In the MOND regime shown in the right column, the broadened distribution is adequately reproduced by the pseudo-Newtonian model with a boosted gravitational constant G eff ≈ 1 . 37 for each of the three samples. The inset shows the arrows indicating the medians. In the MOND regime the medians of the MOND distributions (red arrows) agree well with those of the Gaia data while the Newtonian medians do not. Also, in the transition regime the MOND model performs better than the Newtonian model in reproducing the overall shape of the observed distribution. The legend indicates the K-S p -value representing goodness of fit of each model. (see the texts for further). The dotted curves indicate the contributions of the hierarchical systems only. They dominate the probability distributions for ˜ v > 1 . 5.  \n<!-- image -->  \nGaia medians while the Newton predicted medians start to show deviations in agreement with the analyses of Γ in the previous subsections. The overall shapes of the observed distributions are adequately described by the MOND model as the high p -values indicate (see below for the discussion on the bumps in the high ˜ v tail). On the other hand, the low p -values represent problems for Newton at least for the largest Chae (2023a) sample.  \nI consider next the results for the MOND regime (0 . 8 < s/r M < 1 . 6) shown in the right column of Figure 15. It can be clearly seen that the observed distribution is now further broadened compared with the corresponding Newtonian prediction for all the samples regardless of f multi . This inconsistency with the Newtonian prediction was expected from the median test results with Γ in the previous subsections. The more interesting question is how well the MOND/pseudoNewtonian distribution would match the Gaia distribution. Compared with the Newtonian distribution, the MOND distribution is significantly broadened with a more slowly declining tail in all the samples, and consequently the median is shifted to a higher value in good agreement with the observed median.  \nMoreover, the shape of the MOND distribution for 0 . 8 < s/r M < 1 . 6 captures reasonably well the overall shape of the broadened Gaia distribution for each of the samples, as the high p -values suggest. On the other hand, the Newtonian model has some tension with the new and pure binary samples and is clearly discrepant with the largest Chae (2023a) sample with p -value of < 0 . 0005. Considering that MOND gravity is represented by a simplified pseudo-Newtonian model rather than extensive sets of orbits derived by numerically solving the AQUAL field equation, the near perfect reproduction of the median and the overall adequate reproduction of the shape are striking.  \nWhile the medians and p-values in the transition and MOND regions indicate that the MOND predicted distributions agree overall with the Gaia observed distributions, there is some concern in the details of the distributions because the mild bumps in the range 1 . 5 < ˜ v < 2 . 0 of the observed distributions are systematically above the MOND predictions. Considering that the MOND model represented by a varying G in the transition regime is not satisfactory even in principle, any further investigation of the ˜ v distributions within the present study may be well motivated only for the MOND regime where the boosted gravitational constant is expected to be universal.  \nTo understand the mild bumps in 1 . 5 < ˜ v < 2 . 0 and fully test the Newtonian and MOND models with the ˜ v distributions, I further investigate the full range of the  \nMOND regime given by 0 . 75 < s/r M < 2 . 50 where the lower limit is chosen to insure a uniform gravitational boost factor (see, e.g., Figure 9) and the upper limit is chosen to avoid the mass selection effect (see the range shown in Figure 5).  \nFigure 16 shows the results for the MOND regime. I consider the split ranges of 0 . 75 < s/r M < 1 . 20 (range 1: median = 0.93) and 1 . 20 < s/r M < 2 . 50 (range 2: median = 1.6) as well as the combined range 0 . 75 < s/r M < 2 . 50. The two split ranges provide independent data to test MOND and can be used to investigate any variation of the shape within the MOND range. It can be seen that the medians and the shapes are similar between the two ranges. Most significantly, the difference between the Newtonian median and the observed/MOND median is almost the same regardless of the range or the sample. This is exactly what MOND predicts. All the p -values from either range indicate that the MOND model is acceptable while the Newtonian model is strongly disfavored, in particular by the largest Chae (2023a) sample. The combined range shown in the right column reveals that the Newtonian model is ruled out with p -value from < 1 . 8 × 10 -6 (Chae (2023a) sample) to < 0 . 02(pure binary sample) while the MOND model is well acceptable by the new and pure binary samples and acceptable within 95% by the Chae (2023a) sample.  \nHowever, the concern of the mild bumps in the range 1 . 5 < ˜ v < 2 . 0 remains in the results for the combined range. Comparison of the results for the split ranges indicates that the bumps are largely contributed by the data from range 2. At larger s/r M , the sky-projected velocities v p are lower and thus have larger relative uncertainties for the same required precision of PMs. For those low values of v p even the small fractional uncertainty ( < 0 . 005) of PMs could matter. In the main median tests with Γ carried out in this work, the small uncertainties of PMs are largely irrelevant (see Appendix B for an explicit demonstration) compared with the dominant MC scatters contributed by hierarchical systems and thus were ignored. Here I investigate whether the ignored PM scatters can have any relevance with the bumps.  \nFigure 17 shows the results with PM scatters included in the MC procedure for the Chae (2023a) and new samples. (Note that the PM scatters were already included in producing the results for the pure binary sample shown in Figure 16 because the sample has no hierarchical systems by the statistical definition.) Compared with Figure 16 all medians are somewhat shifted to the right, but the difference between the Newtonian median and the observed/MOND median remains unchanged in", 'MOND regime': 'Figure 16. In this figure, a full MOND range of 0 . 75 < s/r M < 2 . 50 is further investigated. The results for the split ranges and the full range are presented in the same format as Figure 15.  \n<!-- image -->  \nall cases. The p -values for the Newtonian model are still very low essentially ruling it out with these tests alone. Now the tension of the bumps for the MOND model is significantly reduced for the new sample. For the Chae (2023a) sample some tension still remains in the combined range. The remaining minor tension may be due to various factors such as some mismatch between the Chae (2023a) sample and the hierarchical modeling, some unremoved kinematic contaminants at large separations, etc. However, it is clear that the tension is not serious enough in a statistical sense.  \nTable 3 presents the median ˜ v values in the six bins of s/r M used in the previous subsections. The Newton and MOND predicted medians are separately given for the  \nbinary and hierarchical possibilities. When these medians are compared with the Gaia observed medians, the given fractions of hierarchical systems should be used to take an average of the binary and hierarchical predictions.', '4.1. Comparison with Previous Results': 'In this work I have considered the normalized velocity profile as a function of the normalized separation, i.e. ˜ v = ˜ v ( s/r M ). This is a new independent statistical analysis in the plane spanned by two dimensionless quantities ˜ v and s/r M . Unlike Chae (2023a, 2024a), the present method is based on the popular parameter ˜ v that', 'MOND regime: with PM scatters': "Figure 17. Same as the first two rows of Figure 16 but with the small measurement errors ( < 0 . 005) of PMs included in the MC procedure. The high-˜ v tails are somewhat affected by the measurement errors. However, the Newtonian model remains excluded while the MOND model is now slightly more acceptable by the Chae (2023a) sample.  \n<!-- image -->  \nhas been widely used in the recent literature of wide binary gravity tests (e.g., Banik & Zhao 2018; Banik et al. 2024; Pittordis & Sutherland 2018, 2019, 2023; Hernandez et al. 2024). However, this work analyzes ˜ v in the bins of s/r M rather than s . This is because the scaling of ˜ v with s/r M is equivalent to the scaling of g/g N with g N /a 0 as described in Section 2.2.  \nThe present results unambiguously reveal the boost of ˜ v for s/r M ≳ 0 . 7 with respect to the Newtonian prediction. The logarithmic value of the velocity boost factor ranges from Γ = 0 . 05 -0 . 09 depending on the sample and the observational input (in particular eccentricity). The implied boost factor of the observed gravity with respect to Newtonian gravity ranges from γ g = 10 2Γ = 1 . 3 -1 . 5. This is in excellent agreement with several recent results (Chae 2023a, 2024a; Hernandez 2023; Hernandez et al. 2024).  \nThe aforementioned studies reporting gravitational anomaly all considered a broad dynamic range from the fully Newtonian regime to the fully MOND regime and identified the anomaly in the MOND regime after checking/requiring no anomaly in the Newtonian regime. In contrast, Banik et al. (2024) considered a narrow dy-  \nnamic range excluding the fully Newtonian regime and claimed no anomaly using a MOND model as the benchmark of modified gravity. To understand the origin of their result I have investigated a subsample with the same dynamic range of 2 < s < 30 kau taken from the Chae (2023a) sample. With the normalized velocity profile method, I find it unfeasible to simultaneously constrain f multi and the boost/anomaly factor (Γ) with such a narrow dynamic range because MOND gravity models predict a nearly flat Γ for such a narrow dynamic range mostly in the MOND regime and thus there is no clear distinction with Newton.  \nNevertheless, Banik et al. (2024) claimed that their methodology allowed a simultaneous constraint on f multi and Γ 14 based on counts-in-cells in the plane spanned by ˜ v and s , not s/r M . However, because their binary masses are in the broad range 0 . 5 ≲ M tot / M ⊙ ≲ 3 and thus the MOND radii are in the corresponding broad range, their distribution of binaries in the s -˜ v plane mixed  \nTable 3. The binned medians of ˜ v : Gaia observed values versus Newton and MOND predictions  \nNote. 'C23a' and 'new' refer to the Chae (2023a) sample of 19716 systems and the new sample of 6389 systems defined in Table 1. The medians of ˜ v are given for the six bins of s/r M uniformly spaced in the log space as defined in the upper left panel of Figure 3. The Gaia observed medians are compared with Newton('N') and MOND('M') predictions. The medians are given separately for the binary ('-bin') and hierarchical ('-multi') possibilities for each sample. Because it is unknown which ones are binary or hierarchical, all systems in each sample are used to derive '-bin' and '-multi' medians. In deriving each MC set, the medians are calculated either including the individual Gaia PM scatters ('/Yes') or not ('/No'). In comparing the model predictions with the Gaia observed medians, f multi = 0 . 43 / 0 . 41 and 0 . 18 should be used for the 'C23a'(No/Yes) and 'new' samples, respectively.  \ngravity data in a complex way. Further methodological problems of Banik et al. (2024) including too small (compared with data errors) cell sizes are discussed in detail in Hernandez et al. (2023).  \nAlthough f multi cannot be determined for the sample with the range 2 < s < 30 kau, I have considered a fixed f multi = 0 . 65 that is a value 'fitted' by Banik et al. (2024). For this high value, the flat ˜ v profile indeed 'agrees well' with the Newtonian prediction and 'rules out' the AQUAL model because AQUAL predicts a higher median value of ˜ v (Figure 14). Thus, though two (Chae (2023a) and Banik et al. (2024)) samples were selected with different criteria, wide binaries selected within a narrow dynamic range appear to confirm standard gravity if f multi ≈ 0 . 65 is really the case for typical Gaia DR3 samples. This is, of course, not the case because f multi ≈ 0 . 65 is ruled out by the Newtonian regime subsample. See the next subsections for further discussions.", '4.2. Possible Sources of Systematic Errors?': "In this work (and Chae (2023a, 2024a)) I have considered various samples and observational inputs and showed that the gravitational anomaly is robust against  \nsample/input variations. Here I carry out some auxiliary analyses to address any remaining possible concerns.  \nFor the analyses of ˜ v I have considered 6 bins uniformly spaced in log 10 ( s/r M ) that are sufficiently wide to have good number statistics and can still capture any possible variation from the flat Newtonian line in the Newton-MOND transition range of 0 . 8 ≲ s ≲ 6 kau. To investigate the effects of binning on the statistical significance of the inferred gravitational anomaly, here I consider variations in binning. For this test, I consider 4 and 9 bins uniformly spaced in log 10 ( s/r M ) using the new sample that is intermediate in number statistics.  \nFigure 18 shows the results with the alternative bins. These results are qualitatively very similar to the result with 6 bins. As for χ 2 statistics, the AQUAL model has 1 . 1 ≤ χ 2 ν ≤ 2 . 1 and survival probabilities are in the range 0 . 10 ≤ P c ≤ 0 . 36. The Newtonian model has χ 2 ν = 13 . 9, 9 . 4, and 6 . 6 for ν = 3, 5, and 8, respectively. The survival probabilities are similar within the range (0 . 5 ≲ P c ≲ 1) × 10 -8 . These results show that the reduced χ 2 test is largely independent of the specific choice of N bins unless unreasonably small or large numbers are used.  \nOne possible issue that Chae (2023a, 2024a) have not dealt with explicitly is the possible effects of the bi-  \nFigure 18. This figure shows the results with different numbers of bins than the standard choice N bins = 6 for the new sample shown in the left column of Figure 10. The χ 2 ν values of the models change but the survival probabilities remain consistent showing that the gravity test results are not affected by the number of bins in the MC approach.  \n<!-- image -->  \nnary system (i.e. center of mass) velocity v sys relative to the observer. If this velocity is sufficiently large, it can create non-negligible spurious relative velocities between the binary components that are comparable to the intrinsic relative velocities due to the orbital motions. They are called the perspective effects in the literature (Shaya & Olling 2011; Pittordis & Sutherland 2018). For the perspective effects to be important, v sys s/d must be comparable to the relative orbital velocity between the binary stars that has the same order-of-magnitude as the Newtonian circular velocity v c ( s ).  \nFor most binaries in my samples v sys s/d ≪ v c ( s ) is already met, so the perspective effects are not a significant concern. Nevertheless, to explicitly test any possible perspective effects, I consider subsamples with the explicit condition v sys s/d < ηv c ( s ) for η = 0 . 08 and 0 . 04. For this purpose, I use the new sample because RVs as well as PMs are needed. Figure 19 shows the results. It is clearly seen that the magnitude of the gravitational anomaly and χ 2 ν statistics are little affected by the constraints. This is not surprising because 92% of the binaries satisfy even the strong requirement with η = 0 . 04.  \nFinally, I revisit the issue of whether variations of f multi and/or eccentricity with s can be responsible for the measured gravitational anomaly. This issue was dis-  \ncussed in Chae (2023a). In light of the pure binary sample of Chae (2024a) and the new sample selected in this work, I discuss further here.  \nObservational inference of eccentricities using the relative displacement and the relative velocity vector on the sky is now routinely done (e.g. Tokovinin & Kiyaeva 2016; Hwang et al. 2022), and there is no reason to suspect the Hwang et al. (2022) results on the El-Badry et al. (2021) binaries. Nevertheless, I have considered not only individually inferred values but also values from statistical distributions. As a further precaution, I have also considered the thermal distribution. Can there still be a wild variation from these? If any, it may depend on projected separation s and/or the hierarchy. The latter concern may be lessened by considering a sample with a low value of f multi as in the new sample or the pure binary sample. The former concern may be largely removed by considering a sufficiently narrow range of s .  \nThe occurrence rate of hierarchical systems f multi has been assumed to be constant across binary subsamples in this work and other wide binary tests of gravity in the literature. When the range of s is broad as in my samples, this could raise some concern. In particular, if f multi dramatically increases with s across the NewtonMOND transition separation ≈ 2 kau, that could re-  \nFigure 19. In this figure I consider subsamples of the new sample (shown in the left column of Figure 10) that explicitly exclude binaries with relatively larger system (i.e. center pf mass) velocities v sys with respect to the Sun. Because systems with v sys s/d comparable to the internal circular velocity v c can have spurious relative velocities between the components (the so-called 'perspective effects'), I consider η = 0 . 08 and 0 . 04 with the requirement v sys s/d < ηv c . With these requirements only the last bin is somewhat affected, but gravity test results are unaffected in any significant way because the Newton-MOND transition range 0 . 15 ≲ s/r M ≲ 0 . 8 are free of the perspective effects. Note also that my samples already large meet a reasonable requirement (e.g. 96/92% of the new sample satisfy the tested requirement with η = 0 . 08 / 0 . 04).  \n<!-- image -->  \nve the gravitational anomaly in principle (see the next subsection further). However, such a possibility is unlikely for two reasons. First, f multi concerns only unresolved hidden components because all resolved triples and higher-order multiples with detectable brightness have already been excluded from the binary samples. But, all binary stars satisfy the same photometric, astrometric, and kinematic criteria by the selection requirement, and thus the probability of harboring a hidden companion in a star should not depend on the separation s because measured stellar masses are not correlated with s . Second, observed statistics even including resolved multiples do not indicate a strong variation of f multi with s (see, e.g., Figure 11 of Tokovinin (2014), Figure 14 of Hartman et al. (2022)).  \nTo test robustly the gravitational anomaly largely free of the above possible remaining concerns, I consider a narrow range of s that precisely covers the NewtonMOND transition range. From a Milgromian gravity model such as AQUAL, the Newton-MOND transition is expected to occur within the acceleration range  \n10 -10 ≲ g N ≲ 10 -8 m s -2 . For a binary system with M tot = 1M ⊙ , the transition range is approximately 0 . 8 < s < 6 kau or 0 . 1 ≲ s/r M ≲ 0 . 8 for r M ≈ 7 kau. This range is sufficiently narrow that the eccentricity distribution is nearly uniform. Moreover, any possible variation of f multi must be tiny and so I take the uniform value already determined for each source/full sample.  \nAs a further precaution, I consider the limited distance range 100 < d < 200 pc for the current test. There are two reasons for this choice. First, this range out of the full range 0 < d < 200 pc allows the maximum number of binaries within a factor of 2 for the distance difference between binaries. Second, for the combination of s ≲ 6 kau and d > 100 pc, s/d ≲ 3 × 10 -4 and thus any possibility of the perspective effect can be safely avoided.  \nFor the above specifications, I consider a robust and simple slope test without reference to any modified gravity model. The motivation is clear. If the standard gravity holds for this narrow range of 0 . 8 ≲ s ≲ 5 kau, the slope of the anomaly parameter Γ (Equation (20)) must  \nbe consistent with zero regardless of the overall amplitude. A nonzero slope is a sure signature of gravity breaking down. As mentioned above, for this test the new and the pure binary samples are preferable because f multi is small or negligible. However, I also include the Chae (2023a) sample.  \nFigure 20 shows the slope test results for the range of 0 . 1 < s/r M < 1. I consider 6 uniform bins in the logarithmic space -1 < log 10 ( s/r M ) < 0. In all three cases, there is a clear systematic rise from Γ = 0 at s/r M ≈ 0 . 1 to Γ > 0 at s/r M ≈ 0 . 8. The bands representing the combined uncertainties of the slope and the y-intercept do not encompass the Newtonian flat line. The estimated slopes are b = 0 . 093 ± 0 . 035 (new sample), 0 . 169 ± 0 . 045 (pure binary), and 0 . 076 ± 0 . 020 (Chae 2023a). These slopes are discrepant with zero at a significance ranging from 2 . 7 -3 . 8 σ . Here the Bayesian uncertainties are estimated with the Markov chain MC method. An example for the new sample is shown in Figure 21.  \nThe slope test results robustly confirm the breakdown of standard gravity. For the narrow range of 0 . 8 ≲ s ≲ 6 kau (within an order of magnitude) there is no observational reason that the systematic rise of Γ should occur. Although I do not show it here (because it is obvious from various results shown in the results section of this work and the results from Chae (2023a, 2024a)), if a similar test with another narrow range such as -1 . 4 < log 10 ( s/r M ) < -0 . 4 (the same width as above), the measured slopes are well consistent with zero.", "4.3. Newtonian 'Prediction' on the Statistical Properties of Wide Binaries": "All the above results including the narrow range slope test results shown in Section 4.2 have already shown that the breakdown of standard gravity is inevitable. As an auxiliary argument and a matter of completeness, here I quantify what standard gravity requires for Gaia wide binaries. These requirements can be considered as the Newtonian 'predictions' on observations.  \nStandard gravity requires a dramatic change in f multi and/or eccentricity distribution in the transition range of 1 ≲ s ≲ 5 kau. The required changes from the regime s < 1 kau to the regime s ≳ 5 kau are summarized in Table 4. Unless the current estimate of eccentricities by Hwang et al. (2022) is grossly in error, f multi must increase from s ≲ 1 kau to s ≳ 5 kau by ∆ f multi ≈ +0 . 4 in all three samples. This seems impossible mainly because f multi concerns only unresolved hidden companions and all stars satisfy the same photometric, astrometric, and kinematic criteria regardless of s , so that stars should  \nhave nearly equal probabilities of hiding unresolved companions as long as masses are independent of s as is the case for the samples. Moreover, the observed statistics even including resolved companions (e.g. Tokovinin 2014; Hartman et al. 2022) do not indicate at all such a dramatic increase.  \nIf f multi is assumed to be not increasing with s for the whole probed range 0 . 2 < s < 30 kau, eccentricity distribution (and consequently the mean value) must become subthermal ( α < 1) for large separation systems, in the opposite direction to the current estimate by Hwang et al. (2022). The mean eccentricity must change from s ≲ 1 kau to s ≳ 5 kau by ∆ ⟨ e ⟩ ≈ -0 . 10. However, as Figure 24 of Chae (2023a) summarizes various survey results, there is a clear observational trend that ⟨ e ⟩ increases with s .  \nBut, could triples and higher-order multiples affect the eccentricity distribution in a complex way? For example, Tokovinin (2022) reports a mean eccentricity of 0 . 54 ± 0 . 02 for wide outer orbits of nearby resolves triples. This value is similar to the required mean values written in the fourth column of Table 4. However, such an argument cannot work because even if one half of the Chae (2023a) sample are triples the average over pure binaries and triples must be > 0 . 60. Also, it is unreasonable to imagine that resolved triples affect the eccentricity distribution preferentially for s ≳ 5 kau. Moreover, the pure binary sample requires a similar or even lower mean eccentricity although the possible effects of triples are negligible or minimal.", '5. SUMMARY AND CONCLUSIONS': 'In this work I have investigated the normalized velocity profile ˜ v as a function of the normalized radius s/r M as a new study complementary to and independent of two recent studies (Chae 2023a, 2024a). I have considered three samples of binaries including a new independent sample selected with requirements on distances and RVs without using the El-Badry et al. (2021) R values. The samples cover a broad range in number from 3500 to 20000 and in the implied fraction of hierarchical systems from f multi ≈ 0 to ≈ 0 . 4. Through the MC method binaries were distributed probabilistically in the plane spanned by two dimensionless quantities s/r M and ˜ v and at the same time corresponding mock binaries from the virtual Newtonian world were also distributed in the same plane. Then, I derived values of the logarithmic velocity boost parameter Γ (Equation (20)) in bins of s/r M using the binned median ˜ v values from the two MC distributions. The uncertainties of Γ were naturally estimated from a number of MC runs.', 'slope test for the narrow transition range 0.7Ks<5.6 kau with 100 <d<200 pc and individual e ranges': "Figure 20. This figure shows the results from a robust slope test for binaries in the narrow range of the Newton-MOND transition regime. For this narrow range of 0 . 1 ≲ s/r M ≲ 0 . 8 (0 . 7 ≲ s ≲ 5 . 6 kau for the median mass of 1M ⊙ ), any possible concerns in variation of f multi and/or eccentricity with s are minimal. As a further precaution, I also consider the narrow distance range of 100 < d < 200 pc to minimize any possible dependence with distance. These results are pure measurements without reference to any modified gravity theory. The results are particularly robust for the new and pure binary samples because f multi is very low. The measured slopes rule out the Newtonian expectation of zero with statistical significance ranging from 2 . 7 -3 . 8 σ estimated with an MCMC analysis as shown in Figure 21. The light blue bands in the bottom panels represent the combined uncertainties of the slope and the y -intercept. The AQUAL model is well within this band while the Newtonian model is not.  \n<!-- image -->  \nTable 4. Newtonian 'prediction' on the wide binary statistics for 5 ≲ s ≲ 30 kauNote. (1) The required change is with respect to the current estimate for 0 . 2 < s ≲ 1 kau. (2) f multi is estimated for the 'statistical eccentricity distribution' with the acceleration-plane method. (3) α refers to the index in the power-law eccentricity distribution p ( e ) = (1 + α ) e α . (4) The last column quote is based on Hwang et al. (2022).  \nThe derived values of Γ are pure measurements without reference to any modified gravity theory because only Newtonian reference calculations were used to derive them. As shown with a sample of realistic virtual Newtonian binaries, measured values of Γ and their slope as a function of s/r M are expected to be consistent with zero in a universe obeying standard gravity. Because the MC distributions of Γ are consistent with the Gaussian probability density function and each bin has a good number of data points, I have used the widely  \nused χ 2 statistics to test gravity models. I specifically tested the Newtonian and AQUAL models. I have also done a simple slope test with a narrow range of s/r M .  \nFrom the measured values of Γ and statistical tests based on them the following have been found.  \n- 1. The measured values of Γ show a systematic variation with s/r M with respect to the Newtonian flat line. The trend can be easily recognizable even before detailed statistical analyses. Γ varies from  \nFigure 21. An example of the MCMC fit of the measured Γ values in the narrow transition range shown in Figure 20 is exhibited for the new sample. The public Python MCMC package emcee (Foreman-Mackey et al. 2013) is used. The Newtonian model (center of the blue cross) is well outside the probable parameter space.  \n<!-- image -->  \n- ≈ 0 at s/r M ≲ 0 . 15 to 0 . 068 ± 0 . 015 (stat) +0 . 024 -0 . 015 (syst) for s/r M ≳ 0 . 7 or for s ≳ 5 kau. Here the representative value is derived from the new sample of 6389 binaries and systematic uncertainties are due to variations in samples and eccentricities.\n- 2. The gravitational anomaly (i.e. acceleration boost) factor given by γ g = 10 2Γ is measured to be γ g = 1 . 37 +0 . 10 -0 . 09 (stat) +0 . 16 -0 . 09 (syst).\n- 3. The Newtonian model has χ 2 ν = 9 . 4 for the new sample with the survival probability of P c = 5 . 7 × 10 -9 that corresponds to a Gaussian equivalent significance of 5 . 8 σ . For the largest sample of ≈ 2 × 10 4 binaries selected with somewhat relaxed criteria, the Newtonian model has χ 2 ν = 20 . 3 with a Gaussian equivalent significance of 9 . 2 σ .\n- 4. The AQUAL model has χ 2 ν = 1 . 6 for the new sample with the survival probability of P c = 0 . 16 that is fully acceptable. For sample and eccentricity variations, it has 0 . 7 ≤ χ 2 ν ≤ 3 . 1 that are acceptable by conventional standards.\n- 5. The variation of Γ occurs in the narrow range of 0 . 8 ≲ s ≲ 6 kau. For best-controlled samples of binaries from this range and the limited distance range 100 < d < 200 pc, a simple slope test robustly rules out the Newtonian null value in line with the χ 2 test results.\n- 6. The shape of the probability density distribution of the normalized velocity ˜ v in a deep Newtonian  \nregime of strong acceleration is well reproduced by the Newtonian model. However, in a MOND regime of low acceleration, the observed distribution mismatches severely the Newtonian prediction but matches the MOND (represented by a pseudo-Newtonian model with a boosted effective gravitational constant G eff = 1 . 37 G ) prediction (Figure 15).  \nFrom this work and Chae (2023a, 2024a) I have now considered three independent methods with numerous samples. All results consistently point to the immovable gravitational anomaly as long as a sufficiently wide dynamic range is used. The null result by Banik et al. (2024) with f multi > 0 . 6 can be largely attributed to their choice of the narrow dynamic range not encompassing the Newtonian regime and a significant portion of the transition regime 15 (in particular considering the edge effect), but there are other methodological issues as well (Hernandez et al. 2023). Indeed, an independent group uses a sufficiently wide dynamic range and reports essentially the same anomaly through independent analyses based on independent samples (Hernandez 2023; Hernandez et al. 2024). The current evidence for the gravitational anomaly is well above 5 σ and the evidence comes from multiple analyses.  \nThe reported gravitational anomaly is a pure measurement without reference to any modified gravity model. It is striking that the gravitational anomaly agrees remarkably well with a Milgormian (MOND) gravity model. There are no observational reasons that gravity suddenly appears to change from s ≲ 1 kau to s ≳ 5 kau by 40%. It appears that Milgromian gravity is verified as a phenomenological model from current observations. In this respect, it is interesting to note two nearly concurrent independent studies that reached the same conclusion through observations of asymmetrical tidal tails of open clusters (Kroupa et al. 2024) and indefinitely flat rotation curves of rotationally supported galaxies (Mistele et al. 2024). Perhaps, what is even more mysterious is that the world does not know what the underlying fundamental theory (presumably encompassing general relativity; Einstein 1916) is for the low-acceleration MOND phenomenology.", 'ACKNOWLEDGMENTS': "Preliminary results of this work were originally presented by an invited talk in an online meeting 'Challenges of modern cosmology' held on January 18th,  \n2024. 16 The author thanks the anonymous referee for suggesting an additional analysis presented in Section 3.4 and valuable comments to improve the presentations. The author also thanks Xavier Hernandez for  \nsuggesting an improved presentation of figures, and Kareem El-Badry for pointing Tokovinin (2022) to the author for the statistics of resolved triples. This work was supported by the National Research Foundation of Korea (grant No. NRF-2022R1A2C1092306).", 'REFERENCES': 'Hwang, H.-C., Ting, Y.-S., Zakamska, N. L. 2022, MNRAS, 512, 3383  \n(Nbinary  \n<!-- image -->  \nsample: distance  \n200 pc, PM rel error < 0.005, distance rel error < 0.01, RV rel error < 0.5 (Nbinary = 6389)  \nFigure 22. Red and blue points represent the median orthogonal deviations ⟨ ∆ ⊥ ⟩ in three orthogonal bins in the acceleration plane from MC results for the Gaia observed binaries and the corresponding virtual Newtonian binaries. The reader is referred to Chae (2023a) for the definition of the acceleration plane and ⟨ ∆ ⊥ ⟩ . The parameter δ obs -newt represents the difference in ⟨ ∆ ⊥ ⟩ between the observed and virtual binaries. The left columns are the results for the best-fit values of f multi while the other columns represent the uncertainties.  \n<!-- image -->', 'A. RESULTS FROM THE ACCELERATION PLANE ANALYSIS': "Here I present the results from the acceleration-plane analysis in light of a revision of the code from Chae (2023a) and samples. The code from Chae (2023a) is slightly updated (Chae 2023b) in drawing individual eccentricities based on the Hwang et al. (2022) measurements as described in Section 2.4. Figure 22 shows the orthogonal deviations from the g = g N line in the plane spanned by logarithmic accelerations for the Gaia observed binaries and the corresponding virtual Newtonian binaries and the gravitational anomaly parameter δ obs -newt defined by the difference of the two. While these results represent measurements of gravitational anomaly, in the present study they mainly serve to provide measurements of f multi to be used in the normalized velocity profile analysis. I find f multi = 0 . 43 ± 0 . 05 and f multi = 0 . 18 ± 0 . 09 for the two samples (Chae 2023a and new) with individual eccentricities. For statistical eccentricities, they are somewhat lower.  \nFigure 23. This figure shows the acceleration-plane results in 7 bins for the 'Chae (2023a)' sample. The best-fit cases are shown with individual or statistical eccentricities. Note that the algebraic MOND model without an EFE represented by the dotted curve is extremely highly excluded. Any modified gravity model phenomenologically mimicking the algebraic MOND is thus ruled out already.  \n<!-- image -->  \nI also present the results for 7 bins in Figure 23 and 24 for the two samples. In these figures, only best-fit f multi cases are shown but for both individual and statistical eccentricities. Finally, Figure 23 shows the result for the pure binary sample updating Figure 11 of Chae (2024a).", 'B. EFFECTS OF PM SCATTERS': 'In this work, all samples have been selected with the requirement that the measurement uncertainties of all PM components are smaller than 0.5%. In the main parts of testing gravity with median quantities, the tiny measurement uncertainties of PMs have been ignored since scatters of hierarchical systems dominate the uncertainties of ˜ v . Here I present essential results of testing gravity with the acceleration-plane analysis and the Γ (Equation (20)) analysis including PM scatters on top of the scatters of hierarchical systems as a matter of completeness and because PM scatters are considered in investigating the tail of the ˜ v distribution in a study beyond the analyses with Γ as discussed in Section 3.4. When PM scatters are included, it is important to endow the same uncertainties to the mock PMs as to the measured PMs. I find that whenever this is done, the MC procedure works well in reproducing the input gravity law when tested with mock samples.  \nFigure 26 shows the acceleration-plane results for the Chae (2023a) and new samples with PM scatters included in the MC procedure. It can be seen that both the observed and the Newton-predicted accelerations are somewhat shifted upward at lower accelerations due to the added scatters, but the gravitational anomaly (i.e. the difference between the two) remains essentially unaffected compared with the results of Figure 23 and Figure 24. In detail, I note that the fitted value of f multi for the Chae (2023a) sample is slightly lower.  \nFigure 27 further shows the results of testing gravity with Γ. For the Chae (2023a) sample the Newton-predicted ⟨ ˜ v ⟩ now rises more clearly with s/r M , but the gravitational anomaly is little affected because the observed ⟨ ˜ v ⟩ also rises with s/r M by an extra amount thanks to the included PM scatters. Also, the statistical significance measured by χ 2 ν is only weakly affected.', 'new 6389) sample (Nbinary': "Figure 24. Same as Figure 23 but for the 'new sample'.  \n<!-- image -->  \nFigure 25. Same as Figure 23 but for the 'pure binary sample'.  \n<!-- image -->", 'with PM scatters': 'Figure 26. Same as the left columns of Figure 23 and Figure 24 but with PM scatters included in the MC procedure.  \n<!-- image -->  \nFigure 27. Same as the left columns of Figure 9 and Figure 10 but with PM scatters included in the MC procedure.  \n<!-- image -->'}