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

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 85 wide-binaries from the C3PO survey and achieved precisions of $\sim8 - 15$ m,s${\mathstrut}^{-1}$ per binary pair, a $\sim10-100$ $\times$ improvement (median $\sim24\times$) over \textit{Gaia} DR3. After applying a selection criterion based on the scaled velocity, we retain 57 gravitationally bound systems for further analysis. 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 determine that all 57 systems reside in the MOND transition regime ($0.1 < a_N/a_0 < 10$), with none in the deep MOND regime ($a_N/a_0 < 0.1$). Our test is therefore performed entirely in the transition regime, where the choice of interpolating function plays a significant role. We test two commonly used interpolating functions in MOND formulation: the simple form ($b=1$, $\mu = x/(1+x)$) and the standard form ($b=2$, $\mu = x/\sqrt{1+x^2}$). Given our adopted treatment of the external field effect and interpolating functions, our results suggest tension with MOND at the presently accepted $a_0$ value: for $b=1$, the canonical value is excluded at $2.5\sigma$, while for $b=2$, the exclusion is at $1.5\sigma$.