© The Authors 2026

1. Introduction

Wide binaries (WBs) are common stellar systems whose dynamics offer a unique opportunity to test gravity in the low-acceleration regime, in addition to informing stellar and exoplanet studies. Following Gaia data releases (DRs), they have been the subject of a considerable number of studies, particularly after the seminal Gaia catalogue was published by El-Badry et al. (2021). In this Letter we concentrate on the use of WBs as probes of gravity in the regime of low gravitational potential (Hernandez et al. 2012). Wide binaries experience gravitational attraction comparable to that in the outer regions of galaxies. Since in the solar neighbourhood dark matter is negligible (Söding et al. 2025), WBs are suitable probes for directly testing the Newtonian gravitational law in the low-acceleration regime. One limitation is the length of their orbital periods (∼10⁴ − 10⁷ yr), which precludes us from observationally reconstructing the full orbits. A series of studies took advantage of the exceptional Gaia astrometry to create large WB samples with precise Gaia projected separations and transverse velocities (Pittordis & Sutherland 2019; Banik et al. 2024; Chae 2023, 2024) but reached opposite conclusions on whether Newtonian gravity best reproduces the observations. Using HARPS radial velocities, Saglia et al. (2025) have shown that high-quality radial velocities are essential to removing multiple systems from the samples (e.g. where a third body orbits one of the WB stars) and that combining the radial velocities with Gaia DR3 parallaxes and proper motions makes it possible to constrain the orbits of the binaries. These authors analysed 32 WBs, finding that 31 of them could be explained with a viable bound Newtonian solution. Radial velocities have also been used by Chae (2025a), who used Gaia’s best-quality radial velocities and a Bayesian approach to compare Newtonian and MOdified Newtonian Dynamics (MOND) solutions, finding that the latter are preferred for separations beyond 2 kAU. Similar conclusions were reached by Chae (2025b), who applied a Bayesian analysis to the HARPS sample published by Saglia et al. (2025).

In this study we tested whether WBs with large separations can be reproduced by bound Newtonian orbits (‘Keplerian’ orbits hereafter), applying the approach illustrated in detail by Saglia et al. (2025) to a set of archival VLT-ESPRESSO observations of WBs with projected separations larger than 13 kAU. This sample (see Table B.1) complements the WBs studied by Saglia et al. (2025), which have maximum projected separations of 12.3 kAU.

2. Sample and observations

Data were retrieved from the ESO Archive. The stars were observed between December 2020 and March 2022 with the ESPRESSO high-resolution spectrograph at the VLT (Pepe et al. 2021). The spectra are of excellent quality, with a final signal-to-noise ratio (S/N) that varies between 170 and 400 around 570 nm in the combined exposures. ESPRESSO covers the entire 380−790 nm wavelength range at a high resolution (R ∼ 140 000) and is designed to achieve the highest radial-velocity precision (Pepe et al. 2021). In practice, the radial-velocity precision of these observations is limited to a few m s⁻¹ by the S/N of the spectra and by the stellar characteristics (e.g. a high rotation for some targets), fully adequate for the purposes of this study.

Our WBs were taken from the list of solar-type WBs constructed using the Tycho–Gaia DR2 catalogues and were chosen to have projected separations (s) of > 13 kAU (Andrews et al. 2017). The selected binaries all have very similar components, differing by less than 0.3 mag in brightness. Several observations were obtained for each star (a minimum of three and a maximum of ten), for a total of 52 stars in 26 binary systems. All stars are either on the main sequence or just leaving it.

3. Methods

We used the methodology presented in Saglia et al. (2025). Multiple radial velocities, together with metallicity [Fe/H]¹, were used to identify multiple systems and chance alignments (i.e. non-genuine WBs).

We first analysed the radial velocities. Similarly to HARPS, ESPRESSO radial velocities were determined by cross-correlating a digital mask with the observed spectrum. The numbers of observations and masks used are given in Table B.2. The radial-velocity uncertainty was computed by the ESPRESSO data reduction pipeline (Modigliani et al. 2020), and for these observations the error associated with a single measurement is 3−9 ms⁻¹ (see Table B.2). Out of the 26 initial binaries, 12 were discarded because at least one of the components of the system exhibited significant radial-velocity variability. The list of these stars, along with the reasons for discarding them, is provided in Appendix A.

We also used Gaia DR3 to investigate the presence of faint stars in the fields that share the WB’s parallax and proper motions. BD-20655 has a faint (G = 17.18) nearby companion with the same parallax and proper motions; therefore, we excluded this pair (BD-20655 and BD-20658) from the final sample.

The mean re-normalised unit weight error (RUWE) index for the 26 surviving stars is 1.03 ± 0.14, indicating that they behave as single stars with reliable astrometric solutions. Only TYC 5346-457-1 has a RUWE of 1.44, hinting at a possible unresolved multiplicity. Since the radial velocity is very stable (1 m s⁻¹) over seven observations, we kept this star in the sample.

For the determination of stellar parameters and metallicity, we followed the methodology presented in Magrini et al. (2022) and Tsantaki et al. (2025). The method iteratively combines the results of spectral analysis with photometric and astrometric data, comparing them with stellar evolution models to refine the stellar parameters and thus improve the accuracy and consistency of the results. The method starts from parameters derived solely from spectroscopy: effective temperature (T_eff), surface gravity (log(g)), metallicity ([Fe/H]), and microturbulent velocity (ξ). Using this initial set, a new log(g) value is provided based on luminosity and distance data from Gaia, and the stellar mass is determined via a comparison of the stellar parameters with a grid of Parsec isochrones. The spectroscopic analysis is then repeated with log(g) fixed to its new value, yielding an updated set of stellar parameters and masses. The procedure is repeated for two iterations, after which convergence is typically achieved, resulting in parameters that show improved agreement with the isochrone grid. The outcomes of this process are the stellar effective temperature, log(g), metallicity, ξ, mass, radius, and age, Table B.2 summarize our results. The typical precision (1σ) is 4% for the mass and 3% for the radius.

All components of the binaries differ by less than 1σ from the binary mean metallicity. We note that, even if our metallicity estimates are affected by systematic zero-point shifts, the primary quantity of interest is the similarity between the two stars. Evolutionary models also provide age estimates, and one WB (CD-439363 and CD-439366) is very young, less than 1 Gyr. CD-439363 is also a fast rotator, with a cross-correlation peak width of 32 km s⁻¹, which confirms its young age. Since common proper-motion pairs can survive the first hundred million years without being gravitationally bound binaries, we also excluded this WB from the final sample. All other WBs have age estimates well in excess of 1 Gyr, and the relative ages have a typical uncertainty of 25%. This uncertainty grows for lower-mass stars, which can formally appear older than the best-estimate age of the Universe (Table B.2) because their extended main-sequence lifetimes, together with the high sensitivity and limited variation of the fitting parameters during this evolutionary phase (T_eff, luminosity, and [Fe/H]), reduce the precision of the derived ages.

In conclusion, of the 26 WBs observed, 12 have highly variable radial velocities, 1 has a Gaia companion, and 1 is young. That leaves 12 suitable systems, corresponding to 46% of the initial sample. The median radial velocity standard deviation of the final sample is 4 ± 6 ms⁻¹.

The effective temperature and gravity were used to compute the radial-velocity blueshift produced by convective motions in the stellar atmosphere. These corrections were measured by cross-correlating the synthetic 3D atmosphere spectra with the ESPRESSO digital mask (Ludwig et al. 2009; Leão et al. 2019). Gravitational redshift was computed from the inferred masses and radii. The final radial velocity was obtained by subtracting the convective and gravitational shifts from the measured (ESPRESSO) radial velocity: RV = RV_E − RV_con − RV_GR.

The measured radial velocities (RV_E), gravitational redshifts (RV_GR), and convective velocity shifts (RV_con) are reported in Table B.2. Following Leão et al. (2019), in the analysis we added quadratically to the RV_E variance a systematic uncertainty of 40 m s⁻¹, which comes from the uncertainty of the zero-point calibration. Table B.3 reports mean distances from parallaxes, angular separations, masses, projected separations, and velocity differences for the 12 binaries of this study.

4. Bound Newtonian solutions

The geometry and analysis method we followed have been fully detailed in Saglia et al. (2025). These authors show how, adding the third velocity component to the Gaia projected motions and using the Gaia parallaxes and the derived masses with their uncertainties, it is possible to constrain the Keplerian solutions.

For 3 of the 12 WBs studied, we could not find a viable Keplerian orbit. The total velocity difference for these WBs (pairs 7, 11, and 12) is much higher than the maximal Newtonian velocity for bound orbits. This can be seen in the last column of Table B.3, which gives the ratio between the total velocity difference between the components and $\nu =\sqrt{2GM/r}$, the highest velocity difference allowed by Keplerian orbits. For the three WBs without a viable Keplerian solution, we used $\nu =\sqrt{2GM/s}$.

Table B.4 summarises the best-fit distances and the orbital parameters for the nine WBs for which a Keplerian solution could be found. Out of these nine pairs, three (pairs 4, 8, and 9) require an eccentricity above 0.9. For pair 1 we find a peculiar geometry: despite the very small radial velocity difference, the binary is edge-on – the two stars are crossing in front of the observer. The median ratio between the 3D separation (r) and the projected separation (s) for the nine WBs with Keplerian orbits is 1.7, slightly higher than expected from a uniform angle distribution (1.15). The highest r/s ratio is ∼5.4 for the ‘edge-on’ pair (pair 1), confirming the exceptional geometry of this binary.

Saglia et al. (2025) find that the orbital solutions in their sample show an excess of systems with high inclination angles (i), a larger than expected number of orbits close to periastron (ϕ − ϕ₀ ∼ 0), and a thermal eccentricity distribution. Adding the nine WBs of this work does not change these findings.

5. Discussion

About one-quarter of the analysed systems cannot be reproduced with Keplerian solutions, which makes it necessary to examine the possible causes of this discrepancy. For these WBs, the uncertainties in V_tot are less than 10%, and a few percent in ν; we therefore rule out the possibility that the observed high velocity differences are caused by observational errors. Since we cannot demonstrate that these WBs are in bound orbits, before invoking new physics, we explored alternative explanations and the possible causes of this disagreement, which we outline as follows.

5.1. Contamination by unseen companions or giant planets:

Although this hypothesis is very difficult to rule out, the precise radial velocities and the low RUWE parameter would rule out the presence of unseen stellar companions. One should, for instance, invoke special geometry, and this might be the case for pair 7, which contains TYC 5346-457-1. This star has the highest RUWE, and the radial-velocity difference between the components is ∼1/8 of the total velocity difference (see Table B.3). This might indicate a low inclination angle, i.e. an almost face-on orbit. A potential companion of TYC 5346-457-1, co-planar with the face-on WB orbit, would explain the high RUWE and would not produce radial velocity variations. This is the only ‘pair’ in our sample that may have a third companion.

Similar considerations apply regarding the possible presence of massive short-period planets. Of the binaries that cannot be fit with Keplerian solutions, only one star has three ESPRESSO observations, while the other stars have between five and ten. It is unlikely that a massive planet in a close orbit would not be detected, unless in a very peculiar orbit, while a massive exoplanet in a wide orbit would not increase the radial velocity difference by a substantial amount (Jupiter produces a radial-velocity variation of only ∼13 m s⁻¹ on the Sun). Finally, we compared our catalogue with the updated (September 2025) database of confirmed planet-hosting stars (www.exoplanets.eu; Schneider et al. 2011) but found no matches. Considering also the relatively low occurrence of giant planets in short-period orbits, we can reasonably rule this out as a possible explanation.

5.2. Common proper motion pairs from clusters and associations:

We investigated whether any star in our sample is associated with known open clusters or moving groups. We cross-matched our catalogue with the list of open clusters and moving group members compiled by Hunt & Reffert (2023), which is based on Gaia DR3 parallaxes and proper motions, adopting a cross-match radius of up to 5 arcsec. No matches were found. These WBs are not members of known open clusters or moving groups.

5.3. Perturbed, unbound orbits:

Another hypothesis is that WBs are isolated systems. However, since our WBs are several billion years old, their orbits have likely been perturbed by encounters with other stars or molecular clouds and by Galactic disk instabilities (Jiang & Tremaine 2010; Hamilton & Modak 2024). Hamilton & Modak (2024) provide typical timescales for those interactions as a function of the binary’s orbital parameters, finding that an orbit semi-major axis of a ∼ 10⁴ AU marks the transition between the timescale dominated by encounters and that dominated by galactic tides. The disruption timescale for a = 10⁴ AU and a 2 solar mass system is ∼5.7 Gyr (Eq. (18) of Hamilton & Modak 2024). These authors also conclude that orbital eccentricity does not increase with time, while orbital semi-major axes increase with age as an effect of these perturbations. The first conclusion implies that the high eccentricities we now observe were not smaller in the past; the second implies that we may expect an evolution of semi-major axes with stellar age. Finally, Peñarrubia (2021) and Hamilton & Modak (2024) agree that binaries with semi-major axes > 10⁵ AU are disrupted in a short time (≪1 Gyr).

We could therefore argue that the WBs with large semi-major axes or large separations are in non-closed orbits, and that they have been disrupted by encounters or galactic tidal forces. This hypothesis seems reasonable when looking, for instance, at the Jiang & Tremaine (2010) simulations (their Fig. 4) that show a secondary peak in the probability density of escaped binaries for systems with a separation of ∼0.5 Jacobi radii, or ∼170 kAU for 2 solar mass stars, independent of the initial stars’ separation. This simulation shows that even when a WB is disrupted, its stars can remain in quasi-periodic orbits within the Jacobi radius for many gigayears. This separation is larger than that of most of our systems but is comparable at the order-of-magnitude level for several of our WBs (see Tables B.3 and B.4). Also, Hamilton & Modak (2024) show that binaries with initial separations above a few thousand astronomical units are disrupted in a few gigayears (see their Fig. 6). One point that seems at odds with this hypothesis is that the WBs without a Keplerian solution do not have the largest projected separations (s) in our sample (see Table B.3: cf. 10, 22, and 25 kAU respectively for WBs 7, 11, and 12). WB24 in Saglia et al. (2025), the only one that could not be fitted by a Keplerian solution in that sample, has a projected separation of only 5.77 kAU. However, we cannot exclude the possibility that the 3D separation for these pairs is much larger than the projected one. Indeed, this seems to be the case when looking at the probability distribution of the 3D separation for these binaries. Figure 1 shows the r/s versus distance χ²: the most likely 3D separation, r, is several times larger than the projected separation, s, for all WBs without viable Keplerian solutions, and larger than ∼100 kAU for the three WBs of this study. Note also that if these stars were bound, such high r/s values would be extremely unlikely. This supports the hypothesis that they are not bound. We stress that we would miss this information if we used only average distances to the binaries and projected separations.

Fig. 1. r/s (3D/projected separation) ratio vs distance χ2 distribution for the four binaries with no Keplerian solutions: three from this study and one from the HARPS sample (WB24, shown in the top-left panel; Saglia et al. 2025). For all stellar pairs, the most probable r/s are large, from 5 to 50. When multiplying these ratios by the projected separations (s), the 3D minimum separation of the three pairs analysed in this Letter is ∼100 KAU.

Finally, we mention that the orbital eccentricity in WBs has a super-thermal distribution (Hwang et al. 2022), a feature that cannot be explained by orbital evolution (Hamilton & Modak 2024). While a full discussion on the origin of the distribution is beyond the scope of this work, we note that WBs with high eccentricities can be produced as the result of capture (Peñarrubia 2023) or as chance entrapment in tidal tails (Peñarrubia 2021). Our analysis shows that the WBs we studied likely have a common origin, and we verified that they do not belong to known clusters. However, we cannot exclude that their clusters or associations of origin are now dissolved. We believe that our findings add a variable to the eccentricity discussion: since several WBs do not have a Keplerian orbital solution, the super-thermal distribution of eccentricity could have an additional component that has not yet been considered – the observations are interpreted under the wrong assumption because some of these systems are not bound. This could pollute the very wide systems, but we note that the super-thermal eccentricity is already present at projected separations (s) of 1000 AU; at this separation, all the WBs analysed so far can be closely reproduced with Keplerian orbits (Saglia et al. 2025; Chae 2025b).

6. Conclusions

By combining precise radial velocities with Gaia astrometry, even a sample of 12 WBs provides a unique opportunity to study their dynamics. Remarkably, one-third of these systems cannot be reconciled with a Keplerian solution.

Before interpreting these systems as a challenge to Newtonian gravity in the low-acceleration regime, it is essential to explore alternative explanations by improving our understanding of WB evolution and disruption. Current models suggest that the three very wide binaries formed in tighter orbits, which subsequently expanded due to Galactic tides and stellar encounters (Hamilton & Modak 2024). We find that the binaries lacking a Keplerian solution have very large 3D separations, suggesting that they may not be gravitationally bound and could represent disrupted systems. We emphasise that using individual stellar parallaxes, rather than just the projected separation (s) and the system’s average distance, is crucial to revealing the large 3D separations of these binaries, which can be hidden by moderate projected distances.

Although the hypothesis of isolated systems can break at any separation, WBs with 3D separations below 10 000 AU may represent the best limit for testing gravity, because at this separation, the binary destruction time is expected to be comparable to or longer than the stellar age (Hamilton & Modak 2024). Ongoing VLT–ESPRESSO observations will help address this problem. Low-mass stars, which experience similar accelerations at smaller separations, provide another promising target sample despite being fainter and harder to observe. The upcoming Gaia DR4, with improved parallaxes and proper motions, will further refine Keplerian orbital constraints.

Acknowledgments

The authors thank an anonymous referee for very useful comments and suggestions. Based on data from the European Space Agency (ESA) Gaia mission, processed by the Gaia Data Processing and Analysis Consortium (DPAC), https://www.cosmos.esa.int/web/gaia/dpac/consortium. MTM acknowledges the support of the Australian Research Council through Future Fellowship grant FT180100194 and through the Australian Research Council Centre of Excellence in Optical Microcombs for Breakthrough Science (project number CE230100006) funded by the Australian Government. L.B. and L.M. acknowledge support from INAF through the Large Grants EPOCH and WST, funding for the WEAVE project, the Mini-Grants Checs, the PRIN PNRR Project Cosmic POT (2022X4TM3H, MUR) and from the HORIZON-INFRA-2024-DEV-01-01 (Grant No. 101183153). J.C. acknowledges support from ANID via Proyecto Fondecyt Regular 1231345, and BASAL project FB210003.

References

Appendix A: Observed systems not suitable for our test

• CD-28476 and CPD-28132: CD-28476 has a RV variability of 160 ms⁻¹ over 5 observations.

• TYC8862-728-1 and TYC8862-1029-1: TYC8862-1029-1 has only 16 ms⁻¹ RV rms, but the observations in the the second year cluster around a value about 30 ms⁻¹ lower than the observations of the first year. We consider this star a likely long-period binary.

• TYC-6462-304-1 and TYC6462-326-1: The RV TYC-6462-304-1 varies between 23 and 38 km⁻¹.

• TYC7080-660-1 and TYC7080-472-1: TYC7080-472-1 shows an RMS of 15 kms⁻¹.

• HD57376 and HD57257: HD57257 shows substantial variability (100 ms⁻¹) for observations taken 10 days apart.

• CD-68467 and CD-68468: CD-68468 shows 287 ms⁻¹ radial velocity rms.

• TYC7745-788-1 and TYC7745-1033-1: TYC7745-788-1 shows 127 ms⁻¹ spread, with a large difference between one year and next year.

• CD-456209 and HD296796: both stars show large radial velocity variability, probably a quadruple system.

• TYC7644-884-1 and TYC7644-399-1: TYC7644-884-1 shows a variability of 140 ms⁻¹: observations taken 1 year later differ by 300 ms⁻¹

• TYC565-731-1 the RV between the two observations varies by 85 ms⁻¹, with an associated uncertainty of 2 ms⁻¹. The star is likely double.

• TYC8445-1018-1 in four observations within 16 days the radial velocity varies between 993 and 1083 ms⁻¹, but the observations have large associated uncertainties (8-14 ms⁻¹). Our chemical analysis, in addition, shows a difference of 0.5 dex in [Fe/H] with respect to its companion TYC 8445-1073-1.

• TYC9059-134-1 has an rms of 36 ms⁻¹ on 3 observations within a few weeks, each with an associated uncertainty of 9 ms⁻¹. Its companion, TYC9059-376-1 has only one observation.

• BD-20655 and BD-20658: the Gaia DR3 catalogue shows that BD-20655 has a close, faint companion with same parallax and proper motion.

Appendix B: Tables

All Tables

All Figures

Fig. 1. r/s (3D/projected separation) ratio vs distance χ2 distribution for the four binaries with no Keplerian solutions: three from this study and one from the HARPS sample (WB24, shown in the top-left panel; Saglia et al. 2025). For all stellar pairs, the most probable r/s are large, from 5 to 50. When multiplying these ratios by the projected separations (s), the 3D minimum separation of the three pairs analysed in this Letter is ∼100 KAU.

In the text