High-Precision Differential Radial Velocities of C3PO Wide Binaries: A Test of Modified Newtonian Dynamics (MOND)

TL;DR

This study tests Modified Newtonian Dynamics (MOND) in the solar neighborhood by exploiting high-precision differential radial velocities from 100 wide binaries in the C3PO survey. A pixel-integrated forward-modeling approach yields differential RVs of $\sim$8–15 m s$^{-1}$ per binary, enabling fully three-dimensional orbital inferences when combined with Gaia astrometry. A hierarchical Bayesian model jointly infers the orbital elements for all systems and the global MOND acceleration scale $a_0$ under two interpolating functions ($b=1$ and $b=2$). The results show tension with the canonical $a_0$ value, with a stronger rejection for $b=1$ ($3.1\sigma$) and a milder, yet notable, tension for $b=2$ ($1.9\sigma$), indicating that the inferred $a_0$ depends on the interpolating function and that MOND in its standard form may not universally describe wide-binary dynamics.

Abstract

Wide-binary stars, separated by thousands of AU, reside in low-acceleration regimes where Modified Newtonian Dynamics (MOND) predicts deviation from Newtonian gravity. However, Gaia radial velocities (RVs) lack the precision to resolve the small velocity differences expected in these systems, limiting previous MOND analyses to two-dimensional kinematics. In this paper, we introduce a technique to measure differential RVs of wide binary stars using high resolution, high signal-to-noise spectra. We apply this method to measure differential RVs of 100 wide-binaries from the C3PO survey and achieved precisions of $8-15$ m/s per binary pair, a $\sim 10-100 \times$ improvement (median $\sim 24 \times$) over Gaia DR3. Combining these measurements with Gaia astrometry, we construct a hierarchical Bayesian model to infer the orbital elements of all wide-binary pairs and the global MOND acceleration scale ($a_0$). We test two commonly used interpolating functions in MOND formulation: the simple form ($b=1, μ= x/(1+x)$) and the standard form ($b=2, μ= x/\sqrt{1+x^2}$). Our results indicate tension with MOND at the presently accepted $a_0$ value: for $b=1$, the canonical value is excluded at $3.1σ$, while for $b=2$, the exclusion is at $1.9σ$.

Figures (12)

• Figure 1: Per-order differential RV measurements as a function of wavelength for two representative wide-binary pairs. Top panel: Magellan/MIKE observations of the pair Gaia DR3 6473063860275183616 and Gaia DR3 6472230636619614464. The small plot above shows a zoomed in region of this plot. Bottom panel: VLT/UVES observations of the pair Gaia DR3 2919963104221492480 and Gaia DR3 2919963104216198528. Each point represents the $\Delta$RV measured from a single echelle order, with error bars showing the $1\sigma$ uncertainty estimated from the $\chi^2$ curve. The horizontal red line indicates the combined $\Delta$RV obtained from bootstrap resampling of the median, which down-weights outlying orders without requiring sigma clipping.
• Figure 2: Left: Comparison of differential RV uncertainties between our measurements ($\sigma_{\rm meas}$, blue) and Gaia DR3 ($\sigma_{\rm Gaia}$, red) for 84 C3PO binary pairs with Gaia RVs available. Our uncertainties are concentrated at small values, while Gaia uncertainties are broadly distributed at much larger values. Right: Distribution of the ratio $\sigma_{\rm Gaia}/\sigma_{\rm meas}$. The median ratio is $\xi \approx 24$, meaning our measurements have median improvement of $\sim$24 times than Gaia DR3 RVs.
• Figure 3: Transition from Newtonian to MOND regime for two different interpolating functions ($b=1$ is green and $b=2$ is blue). The pink shaded region indicates the range of true separation estimated by our model for the sample we used. The MOND radius ($a_N = a_0$) is marked with red dashed line and the Newtonian limit is marked as gray dotted lines. The transition is sharper for $b=1$ than $b=2$.
• Figure 4: Relative orbits and Keplerian elements for wide binary systems. Left: Three-dimensional view showing the relative orbit of star B around star A (both yellow circles). The coordinate system is defined with the $z$-axis pointing away from the observer (RV direction), and $x$ and $y$ axes aligned with right ascension ($\alpha$) and declination ($\delta$). The orbital orientation is specified by inclination $i$ (angle between orbital plane and sky plane), longitude of ascending node $\Omega$ (where the orbit crosses the sky plane moving away from the observer), argument of periastron $\omega$ (angle from ascending node to periastron), and true anomaly $\nu$ (angle from periastron to current position). The instantaneous separation between stars is $r$. Right: Orbital plane view of the same diagram showing the relationship between eccentric anomaly $E$, true anomaly $\nu$, and the auxiliary circle. The distance $ae$ represents the offset of the focus from the ellipse center, where $a$ is the semi-major axis and $e$ is the eccentricity.
• Figure 5: Graphical representation of the hierarchical Bayesian model for inferring MOND parameter $a_0$ and orbital elements. The global MOND acceleration scale $\log_{10}a_0$ (blue) is shared across all binary systems and inferred jointly from the full dataset. The box denotes plate notation, indicating that the enclosed structure is replicated for each of the $N=100$ binary systems. For each system, we infer six orbital elements: semi-major axis $a$, eccentricity shape parameter $\alpha$ (which determines $e$ via a Beta distribution), mean anomaly $\phi$, inclination $i$, argument of periastron $\omega$, and longitude of ascending node $\Omega$, along with stellar masses $M_1$ and $M_2$. These parameters determine intermediate quantities (gray): true anomaly $\nu$, instantaneous separation $r_{\rm true}$, external field angle $\theta_{\rm ext}$, rotation matrix $R$, and stellar velocities $v_1$, $v_2$. The model predicts four observables (yellow): projected separation $r_{\perp}$, differential radial velocity $\Delta v_r$, and differential proper motions $\Delta\mu_\alpha$ and $\Delta\mu_\delta$. Red nodes indicate probability distributions linking parameters to intermediate quantities or observations. $N_4$ denotes a 4-dimensional Gaussian likelihood over the observables. This model structure applies to both interpolating functions tested: $b=1$ (simple) and $b=2$ (standard), where $b$ controls the sharpness of the transition between Newtonian and MOND regimes in Equation \ref{['eq:mu_general']}.