Dark Energy as a Critical Period in Binary Motion: Bounds from Multi-scale Binaries

The two-body problem under the influence of both dark energy and post-Newtonian modifications is studied. In this unified framework, we demonstrate that dark energy plays the role of a critical period with $T_{\Lambda} = 2\pi/c \sqrt{\Lambda} \approx 60~\text{Gyr}$. We also show that the ratio between orbital and critical period naturally emerges from the Kretschmann scalar, which is a quadratic curvature invariant characterizing all binary systems effectively represented by a de Sitter-Schwarzschild spacetime. The suitability of a binary system to constrain dark energy is determined by the ratio between its Keplerian orbital period $T_\text{K}$ and the critical period $T_\Lambda$. Systems with $T_\text{K} \approx T_\Lambda$ are optimal for constraining the cosmological constant $\Lambda$, such as the Local Group and the Virgo Cluster. Systems with $T_{\text{K}} \ll T_\Lambda$ are dominated by attractive gravity (which are best suited for studying modified gravity corrections). Systems with $T_{\text{K}} \gg T_\Lambda$ are dominated by repulsive dark energy and can thus be used to constrain $\Lambda$ from below. We use our unified framework of post-Newtonian and dark-energy modifications to calculate the precession of bounded and unbounded astrophysical systems and infer constraints on $\Lambda$ from them. Pulsars, the solar system, S stars around Sgr A*, the Local Group, and the Virgo Cluster, having orbital periods of days to gigayears, are analyzed. The results reveal that the upper bound on the cosmological constant decreases when the orbital period of the system increases, emphasizing that $\Lambda$ is a critical period in binary motion.


Introduction
The surprising discovery that cosmic expansion is accelerating has led to the proposal that a mysterious form of dark energy exists, making up 68% of the energy of the observable Universe (Perlmutter et al. 1999;Peebles & Ratra 2003).The simplest model for dark energy, the cosmological constant Λ, assumes a constant energy density that causes the galaxies further away from us to accelerate.Signatures of dark energy emerge from different probes, namely: from supernovae (Riess et al. 1998;Perlmutter et al. 1999) and baryon acoustic oscillations in the late Universe (Cole et al. 2005;Eisenstein et al. 2005;Adil et al. 2021) and from the cosmic microwave background (CMB) radiation in the early Universe (Planck Collaboration VI 2020).All of these probes rely on measurements of the cosmic expansion rate detected at large distances.However, the effect of dark energy could also occur in observations at much lower distances, as recently found by Benisty et al. (2023a) and discussed in numerous other studies in relation to the impact of dark energy in the Local Universe (Chernin et al. 2000(Chernin et al. , 2002(Chernin et al. , 2004(Chernin et al. , 2006(Chernin et al. , 2009(Chernin et al. , 2010(Chernin et al. , 2015;;Baryshev et al. 2001;Kim et al. 2020;Karachentsev et al. 2003;Teerikorpi et al. 2005;Peirani & Pacheco 2008;Teerikorpi & Chernin 2010;Chernin 2015;Silbergleit & Chernin 2019).The cosmological constant changes the predicted mass of the Local Group (LG), assuming that it mainly consists of the Milky Way and Andromeda forming a two-body system (Carlesi et al. 2017;Gonzalez et al. 2014;Li & White 2008;van der Marel et al. 2012;Chernin et al. 2009;Hartl & Strigari 2022;Phelps et al. 2013;McLeod et al. 2017;Lemos et al. 2021).Therefore, it is important for the contribution of dark energy to be taken into account even at such short distances (Peñarrubia et al. 2014).
In an isolated two-body system (as shown in Fig. 1), Λ changes the radial force by a modification with linear dependence of the separation between the two bodies: where r is the separation, M is the total mass of the binary system, G is the Newtonian gravitational constant, and c is the speed of light.It is then possible to quantify the impact of Λ versus the Newtonian force via the ratio between the two terms on the righthand side of Eq. (1): where the Keplerian period, depending on the semi-major axis of the orbit a, is given by: and the critical period related to the cosmological constant is: Fig. 1.Interaction between two bodies in an effective central potential, as considered in this paper: gravitational force (red arrow) as an attracting antagonist against the repulsive force from the cosmological constant (blue arrow).The eccentric anomaly is denoted as η, the true anomaly as θ, and the semi-major axis of the orbit is a.
If κ < 1 (or T K < T Λ ), the Newtonian force prevails.If κ > 1 (or T K > T Λ ), the repulsion caused by the cosmological constant dominates.This condition was included in Benisty et al. (2023a), where it was shown that for the Milky Way and Andromeda to be considered as a binary system, κ ≈ 10%.This implies a dominating Newtonian force for this system on the one hand.However, on the other hand, dark energy still continues to pull Andromeda apart.Since the relative velocities are in the Newtonian regime, Benisty et al. (2023a) did not consider the post-Newtonian (PN) correction to the solution for this binary motion.
In this paper, we extend the approach of Benisty et al. (2023a) to solve the equations of motion for general binary systems, including both the effect of a cosmological-constant term ("post-cosmological", PC) and the PN terms.Thus, contrary to existing approaches, which mostly focus either on the limit of vanishing cosmological contribution to test alternative theories of gravity or vice versa, our unified approach can produce the orbital characteristics for both limits, as well as the transition regime in between.Applying this unified formalism to different bounded systems, including those considered in Jetzer & Sereno (2006), Baker et al. (2015), and Brax & Davis (2018), we find that the upper bound on Λ is closer to the value inferred from the CMB, when the period of the system is close to T Λ .
This paper also discusses the fundamental reason, based on differential geometry, why assuming T Λ as the critical period provides the relevant scale for the impact of dark energy on two-body motions.In the literature, we often come across the vague argument that Λ is only relevant on Gpc scales, implying that a characteristic length scale is most significant, when, in fact, it is the time scale of the orbital period that defines the significance of the dark energy term in binary motion.This is what we aim to show in this work.
The structure of the paper is as follows.Section 2 introduces different astrophysical systems and compares the estimated impact of the cosmological constant on their motion to highlight suitable probes that can put the strongest constraints on Λ on small scales.Section 3 derives the analytical solution of the equations of motion in the presence of both dark energy and the first PN correction.Section 4 calculates the constraints on dark energy from several binary systems introduced in Sect.2, for which the required observational data is available.Subsequently, Sect. 5 summarizes the results.

Astrophysical binary systems
Many astrophysical systems across various size and mass scales can be modelled as binary systems and used as probes to constrain theories of gravity and the impact of dark energy.To evaluate the suitability of a binary system for such studies, two characteristic quantities have to be calculated.The first is the Newtonian gravitational potential Φ, which we scale as: typically inserting r = a, which is the semi-major axis of the orbit.The second quantity is the Keplerian period, T K , as defined in Eq. ( 3).Then, the distinction between different regimes can be made based on the ratios between T K with respect to T Λ and β to 1.

Binary systems considered in this work
In Fig. 2, we plot the following astrophysical systems that can be approximated as binary systems with different β and κ: -Solar System: Since these systems have periods of the order of years or tens of years, which is much smaller than T Λ (T K ≈ 1 yr ≈ 10 −11 T Λ ), the effect of dark energy in these systems is negligible.Kagramanova et al. (2006) constrained dark energy from the solar system and obtained an upper bound of Λ < 10 −37 m −2 .
-S stars: The centre of the Milky Way hosts the closest supermassive black hole, Sgr A*.The stars orbiting Sgr A* are called S stars, with their locations and velocities monitored over the past decades 1 .A large fraction of these stars have orbits with high eccentricities.Thus, they reach high velocities at the pericentre and can be used for constraining scalar interactions.The period of the stars varies from about 10 to about 1000 years, with circa 10 −10 T Λ to 10 −8 T Λ .Therefore, the effect of dark energy in these systems is negligible.The precise bound of dark energy from S stars will be explored in Sect. 4.
-Binary pulsars: The first evidence for gravitational waves was provided by the binary Hulse-Taylor pulsar PSR B1913+16 (Weisberg et al. 1981;Weisberg & Huang 2016).Pulsar systems consist of a neutron star with a white dwarf or another neutron star companion.The monitoring of the times of arrival of the pulsar's radio pulses allows us to infer the properties of the orbit, such as the precession and the emission rate of gravitational waves (Damour & Deruelle 1986).Therefore, they are useful tools for testing gravity due to the extreme precision of the radio pulses they emit.For instance, post-Keplerian parameters which take different forms in different theories of gravity have been constrained from pulsar measurements (Liu et al. 2018;De Laurentis & Capozziello 2011;De Laurentis & De Martino 2014, 2015, 2014;Dyadina et al. 2016;Narang et al. 2023;  4)) sets the scale above which dark energy dominates over gravity.For systems with β ≤ 1, the post-Newtonian correction to general relativistic effects can be observed.Brax et al. 2019Brax et al. , 2021a,b;,b;Davis & Melville 2020;Benisty & Davis 2022;Benisty et al. 2023b).The period of pulsars has been measured to be about 10 −3 yr ≈ 10 −14 T Λ .Their constraint on dark energy is explored in Jetzer & Sereno (2006).In Sect.4, we give the latest update on this constraint with respect to the pulsars' period.
-Local Group: Modelling the LG as a binary system of its two largest member galaxies, Milky Way and Andromeda, the binary motion of these two galaxies forms a bounded system with T K ≈ 17 Gyr ≈ 0.27 T Λ .Thus, as detailed in Benisty et al. (2023a), the system can be used to infer an upper bound on Λ after biases due to small-mass perturbers and the embedding environment have been accounted for.Here, PN corrections are negligible, so that the system cannot serve as a probe for modified gravity tests.
-Bullet Cluster (1E 0657-558): As noted in Bisnovatyi-Kogan & Merafina (2019), binary systems on galaxy-cluster scale are promising probes of Λ, which is why we consider the two receding mass clumps, "bullets", in the galaxy cluster 1E 0657-558.Using the distance between the clumps from Clowe et al. (2006) and the total mass of the system from Lage & Farrar (2014), we estimate T K ≈ 1.23 Gyr ≈ 0.02 T Λ , such that its suitability to constrain Λ is worse than the LG but still worth investigating.
-Galaxy cluster MACS J1206.2:To compare the suitability of an unrelaxed cluster collision, as in the case of the Bullet Cluster, with the rotation curves of a relaxed cluster, we estimate the outermost bound orbit around the centre of MACS J1206.2 to be at its virial radius.Together with the virial mass from Umetsu et al. (2012), we arrive at T K ≈ 8.2 Gyr ≈ 0.14 T Λ .Thus, galaxy-cluster rotation curves could be as promising as LG-like systems in terms of constraining Λ.
-Virgo Cluster and LG: Considering the LG and the Virgo Cluster as a binary system and neglecting the impact of the embedding environment, as done in Bisnovatyi-Kogan & Merafina (2019), we find that the system is unbounded because T K ≈ 172 Gyr T Λ .Hence, it is not a suitable probe for theories of gravity, but a promising probe for Λ instead.Yet, we first have to investigate that the formalism as introduced in Sect. 3 for bounded orbits with perturbative PN and PC corrections can be applied to this system.

Advantages and disadvantages
For binary systems with periods very short compared to T Λ and observing times (e.g.pulsars), their period can be measured many times, resulting in the very high level of precision that is necessary to constrain PN corrections to the unperturbed Keplerian orbits in Newtonian gravity.Besides this, short periods allow us to track the orbits and infer their characteristic parameters from observables, as is the case for the planets of the solar system or the S stars around Sgr A*.Since these probes are also found at small distances away from Earth, they do not require to assume any cosmological background model for their interpretation.As a result, they are less dependent on model assumptions, for instance, about a distance-redshift relation or the presence A83, page 3 of 11 of additional dark matter and its distribution.Therefore, these probes are ideal for constraining PN corrections.
In contrast, objects at higher distances, such as the LG or the bullets in the Bullet Cluster are more suitable for constraining dark energy since cosmic expansion plays a non-negligible role in their motion.These binary systems have longer Keplerian periods, which come at a cost, namely, we only observe a tiny fraction of their period and their orbits in space.Therefore, additional models have to be adopted to solve the equations of motion.Another interesting aspect of these probes is the fact that the LG as a binary system and galaxy rotation curves (e.g. the one for UGC 12632 shown in Fig. 2) are almost on par in terms of their constraining quality.In contrast, rotation curves on galaxy-cluster scale seem to yield better constraints on Λ compared to merging cluster clumps (seen on Fig. 2), comparing the estimates of the Keplerian period for the galaxy cluster MACS J1206.2 and the Bullet Cluster.The reason for this behaviour can be traced to the sizes of the semi-major axes and masses to be inserted into Eq.( 3).
For the same reason, galaxy clusters are not superior to the LG in constraining Λ, as we might have thought.Even if their Keplerian periods are almost of the same order of magnitude, we have to replace the highly precise data from Milky Way and Andromeda by cosmological and other model assumptions to constrain Λ from clusters at cosmic distances.Hence, the confidence bounds on Λ do not need not be smaller than those obtained from the LG measurements.

Relation to the Hubble tension
In determining T Λ , there are two choices for Λ to be inserted into Eq.( 4) to obtain the values shown in Fig. 2. The first comes from Planck Collaboration VI (2020), namely, a constraint on Ω m in a flat Friedmann cosmology that directly yields Ω Λ , which determines Λ when inserting a value of the Hubble constant, H 0 , into Using H 0 = (67.36± 0.54) km s −1 Mpc −1 of (Planck Collaboration VI 2020), the critical period is: Using H 0 = (73.6 ± 1.1) km s −1 Mpc −1 measured from the Pantheon+ and SH0ES data set instead (Brout et al. 2022), we obtain Thus inserting two different values of H 0 which are discrepant at the 5-σ level, we note that their corresponding Λ values have a 3.6-σ difference.As we do not focus on the H 0 -tension in this work, we only use H 0 from Planck Collaboration VI (2020).The critical period in Fig. 2 is also calculated from that value.

Fundamental justification
All binary systems introduced in Sect. 2 can be treated by the same formalism as an effective central potential and a test body with reduced mass orbiting in this (PN) potential, with or without including dark energy terms.The joint, general line element ds 2 for the effective central potential and the effect of dark energy is expressed as: in which θ and φ are the usual spherical angles and captures the gravitational-potential and dark-energy contributions.It is thus a de Sitter-Schwarzschild metric that is reduced to the Schwarzschild metric in the limit of Λ = 0. Equation ( 10) directly contains the dimensionless parameter, β, (Eq.( 5)) that will later be used to track PN corrections.Expressing β in terms of the Schwarzschild radius, r s = 2GM/c 2 , gives r s /(2r).
The role of κ becomes clear when we determine the quadratic scalar that quantifies the gauge-invariant curvature of spacetime, namely, the Kretschmann scalar, defined as: in which R λ αβγ is the Riemann tensor.For the de Sitter-Schwarzschild metric, we obtain using r = a in the last step and introducing the Keplerian angular frequency ω K = 2π/T K .The curvature scalar can thus be written as a combination of the two periods T K and T Λ .Alternatively, we can rewrite the periods in terms of angular frequencies, such that κ = (ω Λ /ω K ) 2 , to avoid introducing factors of 2π or express K as: in which K s = 48 G 2 M 2 /(c 4 r 6 ) is the Kretschmann scalar for the Schwarzschild metric.Equation ( 13) shows that κ is the characteristic term to describe the impact of Λ when we embed the twobody motion, effectively modelled as orbits in a Schwarzschild metric, into an expanding de Sitter space.So, the ratio between the Keplerian orbital period to the critical period T Λ , or the corresponding ratio of angular frequencies, arising in the solution of Benisty et al. (2023a) emerges naturally from the more fundamental origin of a curvature invariant in general relativity.This implies that κ as a ratio of time-related quantities determines the impact of dark energy on a binary system and not a ratio of length scales, as is often given.After deriving Eq. ( 12), it is clear that our overview of probes of PC and PN corrections in Fig. 2 can be directly linked to Baker et al. (2015) because this work also used Φ/c 2 and the Kretschmann scalar to characterize tests of gravity on various scales.

Critical frequency versus critical distance
There are different criteria to test the dominance of dark energy.One of them is based on the zero-velocity surface (Kim et al. 2020), while another is based on the zeroacceleration surface (Pavlidou & Tomaras 2014;Pavlidou et al. 2014;Tanoglidis et al. 2015).In principle, however, all of these criteria use the same equilibrium in which the outward-directed influence of Λ is balanced by the inwardly directed gravitational A83, page 4 of 11 attraction as reference.For instance, the zero velocity surface is expressed as: which corresponds to κ = 3 in our framework.Inserting Eq. ( 14) into Eq.( 12), we see that the ratios of critical frequencies are equivalent to a ratio of length scales.Benisty et al. (2023a) shows that the period of the LG is about 30% of the critical period, hence, dark energy is important in this binary motion.This is equivalent to the condition that the distance between the Milky Way and Andromeda galaxies is ∼0.77Mpc, which is of the same order of magnitude as the radius of the zero-acceleration surface of about 2 Mpc.We prefer to use the critical frequency, since the critical distance, r V , depends on the mass of system and Λ.In contrast to this, the frequency condition separates cosmology from the system considered.It is a direct comparison between the critical one, only depending on Λc 2 , and the Keplerian one, only depending on the characteristics of the binary system.Thus, similar to defining a Planck time, the critical period only depends on a cosmological parameter independent of the system under consideration.

Lagrangian and action
As also discussed in Damour & Schaefer (1988), Bisnovatyi-Kogan & Merafina (2023), and Jetzer & Sereno (2006), after transforming into the centre-of-mass system, the Lagrangian per reduced mass is expressed as: with the following three individual components The action L PN 1 is the first PN correction to the Newtonian term L 0 (Damour & Deruelle 1985).The cosmological-constant action in binary motion is L Λ (Carrera & Giulini 2006).We denote the reduced mass as µ and the total mass of the system as M. In addition, r is the separation between M and µ and v their relative velocity.To simplify the notations, ν ≡ µ/M is the ratio between the reduced mass and the total mass with 0 ≤ ν ≤ 1/4 (Damour & Deruelle 1985).For a binary system whose motion occurs in the polar coordinates (r, θ), we relate v and r with the derivatives of r and θ to obtain: Here, n denotes the unit vector in direction of the separation and r and u are the separation and velocity vectors, respectively.Consequently, the energy per units of reduced mass is: with the first PN correction as: The angular momentum per units of reduced mass is:

Effective potential
In order to solve the equations of motion, we isolate ṙ and θ from the conservation of the energy and the angular momentum (see Sect. 3.5 below).To do so, we write: with coefficients The cosmological constant adds a quadratic correction to the potential and a linear term from the interaction with the PN correction.For the second equation, we introduce having γ −1 = 0.
The effective potential can be rescaled with respect to the Schwarzschild radius using the dimensionless quantities: yielding the dimensionless effective potential: where we take ν = 0 without loss of generality.Since the maximum ν attainable is 1/4, all signs in Eq. ( 20) remain unchanged for the entire interval of admissible values for ν.Hence, while ν = 0 only yields the effective potential for a massless test particle, the general case of ν 0 is omitted for the sake of brevity and simplicity of equations shown.
Figure 3 shows Eq. ( 23) for different values of Λ.We take Λ = 0, 10 −3 , and 2 × 10 −3 for the blue, yellow, and green curves, respectively.For all of them, we chose /c 2 = −0.15 and j = 2.65 as arbitrary, representative example values.For small radii, the term ∼r −3 is dominant.The first maximum in the effective potential which is caused by the interplay between the r −3 and the r −2 terms does not change much when changing Λ.At larger radii, the quadratic shape of the potential depends on the value of Λ.
A83, page 5 of 11 Fig. 3. Effective potential including post-Newtonian (PN) and cosmological-constant contributions scaled by the Newtonian potential at the Schwarzschild radius r s = 2GM/c 2 versus the radius scaled by r s .For Λ = 0 (blue curve), stable orbits are reduced to the PN case.For small, perturbative Λ (yellow curve), stable orbits can only become unbounded beyond a certain radius, where the impact of Λ dominates.For large, dominating values of Λ (green curve), there are only unbounded orbits.
For Λ = 0, stable orbits are reduced to the PN case.For small, perturbative Λ = 10 −3 , stable orbits can become unbounded beyond a certain radius, where the impact of Λ dominates.For the value of Λ = 2 × 10 −3 , there are only unbounded orbits.

Equations of motion
Similarly to Damour & Schaefer (1988), we use the transformation on r: that simplifies the potential for the ṙ eliminating the r −3 term.
The equation of motion is reduced to with the following abbreviations: Correspondingly, the transformation simplifies the equation of motion for the true anomaly

Closed orbits
In order to find the relation between the semi-major axis a and the eccentricity e to the variables in Sect.3.5, we use the condition where the ṙ is zero: Imposing this condition on Eq. ( 25) yields the relations for the energy and the angular momentum: −βκ e 2 (6ν + 1) − 14ν + 15 6 . (29) For the limit κ → 0, we retrieve the same relations for the modified energy and angular momentum as for the pure PN case.Analogously, for β → 0, we obtain the modified relations for the pure PC case as detailed in Benisty et al. (2023a).For closed orbits, we parameterise the solution of r as: in which η is the mean anomaly angle and ā and ē are the transformed semi-major axis and eccentricity, respectively.Integrating Eq. ( 25) under the assumption that β and κ are sufficiently small, we obtain: 2π T (t − t 0 ) = η − e t sin η + κ(g 2 e 2 s 2 + g 3 e 3 s 3 ), where s 2 = sin 2η, s 3 = sin 3η, and analogously for higher orders not shown here.We also abbreviate: The modified temporal eccentricity and the period are expressed as: The g 2 and g 3 couplings are due to the κ order and inside their definitions is a modification of order of βκ.Using Eq. ( 27), parameterising r via the relation we obtain the solution for the true anomaly: A83, page 6 of 11 with e φ = e 1 + ν 2 β and k being the pericentre advance parameter introduced in Eq. ( 36).For a complete period running from ν = 0 to ν = 2π, the excess angle is Here ω is the argument of the pericentre.For κ → 0, we get the same known precession term from the PN correction.

Unbound orbits
For very large separations, r, the Λ terms, and thus cosmic repulsion, dominate in Eq. (1).To obtain an analytic solution for this regime, we consequently set up the radial equation of motion for an unbound orbit to be: Taking the derivative of this equation with respect to t and assuming that ṙ 0, we obtain: The ansatz r(t with coefficients A and B determined from initial conditions is easily verified as a solution to Eq. (37).To fix the boundary conditions that determine A and B, let us assume that r(t = 0) ≡ r0 and ṙ(t = 0) ≡ v0 .The solution to Eq. ( 37) is then expressed as: The plus or minus sign can be chosen depending on the initial direction of the motion.Next, we use Eq. ( 24) to transform Eq. ( 39) back into the original coordinate system to find a physically reasonable initial position.Requiring that r(t = 0) ≥ 0, we find that r0 ≥ (2 − 3/4ν) r s .Depending on the energy of the reduced mass, k , the initial velocity is fixed as well, given the initial position.At the maximum ν = 1/4, we find that the initial condition for the position requires to start at a radius of at least 1.625 r s , implying that the reduced mass is on the right hand side of the maximum effective potential belonging to the last stable orbit around the central mass as plotted for dominating Λ in Fig. 3. Next, we solve the angular equation of motion for an unbound system.Similarly to neglecting the terms inversely proportional to a power-law of r, we neglect the first term on the right-hand side of Eq. ( 27) for an effective potential with dominating Λ.Then, θ ∝ jΛ.Treating j as a conserved and known quantity, the solution of the angular equation of motion is simply: assuming that the motion started at t = 0 for consistency with the solution for r(t).For the more general case of intermediate separations, the equations of motion have to be solved numerically.This is, for instance, done in Kim et al. (2020) to study the interplay between gravitation and dark energy for the Virgo Cluster (see Sect. 4 for details).
Noting that an unbound orbit can also be understood as a scattering process with the reduced mass being scattered off the effective central potential, our approach can be connected to the one sketched in Hertzberg & Loeb (2023).The authors also use a de Sitter metric as an embedding for their scattering scenario.Then, they employ a timing argument that the scattering black hole needs to exist for at least twice the Hubble time to constrain Λ of the de Sitter embedding based on the fact that the evaporation time of a black hole is linked to its mass and this, in turn, can be inferred from the scattering process.Most notably, their calculations arrive at a "lower" bound on Λ.Hence, transferring this argument to our unbound system of the LG and the Virgo Cluster, it is obvious that follow-up studies for this system should also obtain a minimum Λ, necessary to turn this binary system into an unbound one (see also Sect.4).

Osculating formalism
Another way to approach the orbits is the osculating formalism (Poisson & Will 2014).It can be used to describe the orbits of a body moving in a gravitational field, taking into account all the perturbing forces acting on it at the moment under analysis.Keplerian orbits are assumed to be tangential to the true orbit at any moment.The equations that govern the orbital parameters in our case can be summarised as follows: , where the perturbing force is described by To shorten the equations, we set c 1 ≡ cos f, s 1 ≡ sin f , in which f is the true anomaly θ with an initial value of θ 0 and ω is the angle between the ascending node and the pericentre.The integration of the osculating equations gives two types of changes: an oscillation with a period equal to the orbital period and a steady drift of the orbital elements that do not average out after few periods.For the secular drift, the integration over [0, 2π] yields We use this formalism to confirm the results obtained for the closed orbits in Sect.3.6, see Eqs. ( 33) and (36).As expected, the terms match the derivation in our unified framework.Here again, we see the explicit dependence on κ relating the orbital elements to T Λ .

Bounds on the cosmological constant
Figure 2 shows the relation between the orbital periods of different systems and T Λ .The closer the orbital periods are to T Λ , the tighter the constraints on the cosmological constant will be.
A83, page 7 of 11 However, as briefly noted in Sect.2, for binary systems fulfiling this criterion, such as galaxies and galaxy clusters, we only have a snapshot of the orbital motion because possible observation times are much shorter than the orbital periods.For systems with much shorter orbital periods, such as pulsars or stars orbiting around Sgr A*, we are able to track the entire orbit for several periods and constrain the corresponding precession.Yet, despite the increased precision, the lower orbital period will weaken the constraint that can be put on Λ from such systems.
In this section, we show the correlation between the precision and the bound on Λ based on the observed periods of binary systems introduced in Sect.2, taking their measurement precision into account.To solve the equations of motions of Sect.3.5 for the individual systems given the measurement precision of the velocities and distances to be inserted into the formalism, we use a χ2 -approach with a flat prior of κ ∈ [0, 100], or Λ ∈ [0, 100] × ω 2 K /c 2 , where ω K is the orbital angular frequency of the system.We use an affine-invariant nested sampler for the minimisation of the likelihoods as implemented in the open-source package Polychord (Handley et al. 2015) and run it in the standard configuration, satisfying the needs of our optimisation problem.

Solar System
Upper bounds from the Solar System were discussed in Wickramasinghe (1999), Kagramanova et al. (2006), Sereno & Jetzer (2006), Iorio (2012) and Liang & Xie (2014).These different methods give different ranges of Λ.In our analysis, we use those planets whose orbital precession ωm has been measured: Mercury, Venus, Earth, Mars, Jupiter, and Saturn.We also take their planetary parameters from a NASA fact sheet 2 .All observations are inserted into: to minimise the difference between the measured and modelled precession per the orbital period ω = k/n, which is a function of the planetary parameters and Λ.The posterior on Λ for Mercury reads: and the one for Saturn is: In Fig. 4, we plot the upper bound for different planets over their orbital periods.Already for the Solar System planets alone, the figure and the numbers show that the upper bound on Λ is decreasing for increasing orbital period.

Galactic Centre
Using the known precession of the S2 star orbiting Sgr A*, we forecast a bound on Λ with a χ 2 similar to the one used for the solar system.The posterior on Λ for S 2 yields: Comparing with the values from the solar system, this bound is of the order of the bounds from the outer planets, since the orbital period of S2 is about the same as Saturn.

Pulsars
PSR J0737−3039A/B is the only known double pulsar with associated very high-precision measurements.The system has been studied continuously using a number of radio telescopes, with improved data acquisition systems and better sensitivity, resulting in much improved timing precision over time.
The latest measurement of PSR J0737−3039A/B is published in Kramer et al. (2021) and includes higher orders in the PN expansion to guarantee a high precision.The masses of the pulsar, m p , and its companion, m c , (which is another pulsar for the PSR J0737−3039A/B event) have to be determined from the observables together with Λ.To do so, we calculate the power radiated in a de Sitter universe (Bonga & Hazboun 2017): where M c = (m p m c ) 3/5 /(m p + m c ) 1/5 is the chirp mass and n is the orbital frequency.The χ 2 for the double pulsar is then: in which ξ = ( ω, Ṗ, q) is the vector of the post-Keplerian parameters containing the measured values and δξ represents the measurement uncertainties.The ratio q ≡ m p /m c contains the mass of the pulsar and its companion.As priors for the post-Keplerian parameters, we use the Gaussian priors reported in the original papers.For the masses, we set a uniform prior of [0, 4] × M .The posterior for Λ then yields: Despite the high precision of the pulsar observations, the upper bound is less tight than the one from the solar system or the S2 star due the shorter period of the system.

Local Group
The LG includes the Milky Way, the Andromeda galaxy, along with at least 78 other known members, most of which are dwarf galaxies.It has a total diameter of about 3 Mpc with the total mass of few 10 12 M .The two largest members, the Andromeda Galaxy and the Milky Way, are both spiral galaxies of about 10 12 M , with each hosting its own system of satellite galaxies.
A simple model of the LG dynamics reduces the LG to these two galaxies as a two-body system.Since the Keplerian period of the Milky Way and Andromeda (17 Gyr) is about 0.27 T Λ , Benisty et al. (2023a) derives an upper bound on Λ to be 5.44 times the Λ-value obtained by Planck.The χ 2 uses the mass of the LG and is taken from Benisty et al. (2023a): where M(r, v r , v t , Λ) is the predicted mass as a function of the separation r towards Andromeda, the tangential velocity v t , the radial velocity v r , and Λ. The, M m ± ∆M m is the mass and its uncertainty inferred from the data as detailed in Benisty et al. (2022).It amounts to: The posterior distribution includes the Planck value within less than one σ.

Virgo Cluster
The Keplerian orbital period for the binary system consisting of the LG and the Virgo Cluster is calculated from the total mass of the system being of the order of M ≈ 10 15 M and the semi-major axis of the orbit a ≈ 16 Mpc.The local group with approximately 10 12 M only makes a minor contribution to the total mass.Besides this, we use the distance between the LG and Virgo as an approximation for the semi-major axis.
For the obtained T K , κ > 1.Consequently, we conclude that this system is unbounded and derive a lower bound on Λ from Eq. ( 2), which yields While this estimate only considers the LG as a test particle in Virgo's extended, effective potential, the complex structure of the Virgo Cluster also allows us to probe the potential for bounded orbits.Using 33 galaxies along the line of sight from our LG to the centre of the Virgo Cluster, Kim et al. (2020) probe the effective potential to determine the radius at which gravity is balanced out by dark energy and H 0 from their infall model towards Virgo.The equations of motions are analogous to the ones in Sect.3.5, so that we can convert their constraint on H 0 to a constraint on Λ using Eq. ( 6) to obtain This value is larger than the lower bound we calculated from Eq. ( 2) and thus in agreement with our theoretical derivations.

Planck versus SH0ES values
Determining Λ from the expansion rate of the universe (Eq.( 6)) yields different values due the different values of the Hubble constant.The value from the Planck observations is: while the value from the SH0ES data yields Even though these two different values are in the range of 10 −52 m −2 , they have a 3-σ difference.

Bounds versus orbital frequency
Figure 4 summarizes all upper bounds on Λ for the different systems in terms of 1-σ upper bounds obtained by minimising the χ 2 with respect to their measured orbital periods.Constraints on Λ improve with increasing precision of the measurements.Therefore, we observe that the upper bounds on Λ for the highly precise pulsar observations are far below the κ = 1-line (dashed black line in Fig. 4) which separates the bound orbits from the unbound ones.For the latter, T K = T Λ by Eq. ( 2), such that the attractive gravitation is balanced by repulsive dark energy (or A83, page 9 of 11 zero-velocity surface as we discussed).Yet, even though the pulsar observations have a higher precision than the planets in the solar system or the S stars, the latter yield a stronger upper bound on Λ.The reason for this can be found in the fact that the periods for the planets in the solar systems are larger and the upper bound on Λ comes closer to the κ = 1-line, which is decreasing for increasing periods.Given that no bound on Λ can cross this line without turning the bound system into an unbound one, it becomes clear that binary systems with increasing periods will yield tighter constraints on Λ.
Considering unbound systems such as the LG and the Virgo Cluster, the tightest lower bounds on Λ can analogously be achieved close to the κ = 1-line for systems with orbital periods close to T Λ .Increasing the periods to much larger values, we arrive at regimes in which the cosmic expansion dominates and the orbiting bodies cannot be considered as a (perturbed) binary system anymore.In this regime, Λ is rather constrained by supernovae or other large-scale observations.
As an outlook, Iorio (2018) discusses the possibility for an upper bound on Λ from future measurements of a pulsar in the vicinity of Sgr A*.According to our criterion of κ, the upper bound from this system will be optimal if the period of the pulsar is as long as possible and, at the same time, achieving the highest possible precision.

Conclusion
In this paper, we set up a general and unifying framework to solve the equations of motions for binary systems including both, the impact of dark energy, and the first post-Newtonian modification.We show that T Λ = 2π/c √ Λ ≈ 60 Gyr is a critical period for these systems.If the orbital period is much lower than T Λ , dark energy only perturbs the binary motion and the gravitational attraction remains dominant.If the orbital period is much larger than T Λ , the system is not bounded and the expansion of the universe, meaning the repulsive dark energy, is dominant.
As detailed in Benisty et al. (2023a), the ratio between the Keplerian orbital period T K and T Λ , turning T Λ into a critical time scale, arose from the Newtonian equations of motion when adding Λ.However, all binary systems we consider can also be effectively described as Schwarzschild-like potentials embedded in a de Sitter space-time.In doing so (see Sect. 3.1), we find that the same ratio between T K and T Λ naturally arises as the change in the Kretschmann scalar when embedding the Schwarzschild metric into the de Sitter space-time to account for the impact of Λ on the binary system, as seen in Eq. ( 13).Thus, the ratio between the orbital period and the critical period is based on differential geometrical foundations beyond the purely Newtonian framework.
Hence, after the systematic analysis of binary systems in this unifying framework, we conclude that the solutions are characterised by two parameters: the dimensionless gravitational potential β = Φ/c 2 , selecting between the Newtonian and post-Newtonian regime, and the ratio κ = (T K /T Λ ) 2 , determining the influence of Λ on the binary motion.
Considering measurements from different binary systems having orbital periods from days to gigayears, we show that the upper bound on the cosmological constant decreases when the orbital period of the system increases.We derived upper bounds3 on Λ from the precession of Solar System planets, the precession of the S2 star around the galactic center, and the precession of the double pulsar PSR J0737−3039A/B.Despite its more precise measurements, the upper bound on Λ from the double pulsar is less tight than the one obtained from S2 or Saturn, which yield the tightest constraints among these systems.This result is explained naturally in our generalised framework: As κ = 1 separates bounded from unbounded systems, the increasing orbital periods of the S2 star and Saturn (compared to the double pulsar) require the upper bounds on Λ to be below the κ = 1-line in Fig. 4 to avoid turning bounded systems into unbounded ones.Consequently, the constraints on Λ from bounded systems become tighter, the closer T K approaches T Λ .
The Local Group (LG) has the longest orbital period in bounded binary motion that we discuss in this paper.Since the orbital period of the LG is about 0.27 T Λ (κ ∼ 10%), the upper bound on Λ is expected to be tighter than the ones from the solar system and S2.Benisty et al. (2023a) found the upper bound on Λ to be 5.44 times larger than the Planck value reported in Planck Collaboration VI (2020), which is an impressive result for an individual system in our local cosmic neighbourhood compared to the all-sky observations from the early universe.As Fig. 2 shows, further applications of our framework introduced here to galaxy rotation curves or even galaxy clusters should yield constraints on Λ within the same precision range, which will be subject of further studies.
As T K < T Λ for bounded orbits, they yield upper bounds on Λ, while unbounded orbits with their T K > T Λ yield lower bounds on Λ. Strictly speaking, the latter case does not have "orbits" anymore, so we can also interpret it as a scattering of the test particle in the effective central potential; this is, for instance, similar to the case considered in Hertzberg & Loeb (2023).One example for an unbounded motion considered here is the system of the LG and the Virgo Cluster.Kim et al. (2020) probe the effective potential of the Virgo Cluster at increasing distances out to the LG to find the distance at which the gravitational attraction by the Virgo Cluster is balanced by the repulsion from Λ. Hence, this is the first binary system from which an upper and a lower bound on Λ can be inferred and probably the one with the highest precision.While Λ = 0 could not have been excluded for the LG in Benisty et al. (2023a), current estimates for the Virgo Cluster, Eq. ( 54), based on the results of Kim et al. (2020) already challenging a vanishing Λ.
Attempting to constrain dark energy models beyond a constant Λ is more complicated, as binary systems with the right order of magnitude in terms of orbital periods, T K ≈ T Λ , are usually subject to a lower observational precision (to be performed at larger distances from our position).The evaluation also involves more model dependencies, as the increased orbital periods only allow us to observe a tiny fraction of the orbit.Additionally increasing the number of degrees of freedom that are to be constrained will lead to weaker confidence bounds due to degeneracies.Thus, we should first aim to constrain Λ as tightly as possible in follow-up studies before turning to more complicated models.

Fig. 2 .
Fig. 2. Scaled Newtonian gravitational potential, β = GM/ac 2 , versus the measured and estimated orbital periods for various astrophysical systems.The critical period T Λ = 2π/c √ Λ ≈ 60 Gyr (Eq.(4)) sets the scale above which dark energy dominates over gravity.For systems with β ≤ 1, the post-Newtonian correction to general relativistic effects can be observed.

Fig. 4 .
Fig. 4. Comparison between bounds on the cosmological constant for different systems versus the orbital period of these systems: Planets in the solar system (dark blue), S2 star around Sgr A* (green), double pulsar PSR J0737−3039A/B (azure), Local Group (red), and Virgo (pink).Each upper bound is the 1-σ upper bound on Λ.For comparison, we put the estimated values of the cosmological constant from Planck Collaboration VI (2020; grey) and SH0ES,Brout et al. (2022; yellow).For longer periods, the upper bound on Λ decreases and proves that systems with orbital periods in the range of T Λ ≈ 60 Gyr have to be used to constrain Λ from binary systems tightly.The dashed black line is κ = 1, which represents T K = T Λ and thus separates bound from unbound systems.