Eclipse Timing the Milky Way’s Gravitational Potential

Measurements of the accelerations of stars provide the most direct probe of the mass distributions (the stars and the dark matter) of galaxies. For more than a century now, the fundamental parameters that describe the Milky Way (MW) have been determined using kinematic estimates of the accelerations of stars (Oort 1932; Kuijken & Gilmore 1989; Bovy & Tremaine 2012; McKee et al. 2015; Schutz et al. 2018) that live within its gravitational potential. Recent advances in technology have led to the development of several different techniques to measure Galactic accelerations directly using extremely precise time-series measurements.

Extreme-precision spectrographs that can achieve an instrumental precision of ∼10 cm s⁻¹ (Pepe et al. 2010; Wright & Robertson 2017) have opened up a new avenue to measure the Galactic acceleration directly (Silverwood & Easther 2019; Chakrabarti et al. 2020). The analysis of ongoing pulsar timing observations has recently enabled measurement of Galactic accelerations (Chakrabarti et al. 2021) and constraints on the Galactic potential, including a measurement of the Oort limit (the midplane density), the local dark matter density, and the oblateness of the Galactic potential as traced by the pulsars. Gaia astrometry was analyzed to measure the solar system acceleration (Gaia Collaboration et al. 2021), which is in agreement with earlier VLBI measurements (Charlot et al. 2020). In this Letter, we develop a framework to use eclipsing binaries (EBs) to measure Galactic accelerations directly.

There is now a plethora of observed phenomena, in both the gas and the stellar disk, that indicates that our Galaxy has had a highly dynamic history (Levine et al. 2006; Chakrabarti & Blitz 2009, 2011; Xu et al. 2015; Helmi et al. 2018; Antoja et al. 2018). This dynamic picture of the Galaxy has come into especially sharp focus with the advent of Gaia data (Gaia Collaboration et al. 2016). Analysis of interacting-MW simulations shows that there are differences between the “true” density in these simulations and density estimates from the Jeans analysis (which assumes equilibrium) (Haines et al. 2019); this underlines the need for direct acceleration measurements that are based on time-series observations. The acceleration profiles in interacting-MW simulations are highly asymmetric, in contrast to static potentials or isolated MW simulations (Chakrabarti et al. 2020). For a galaxy with a dynamic history like the MW, kinematic analyses based on snapshots of stars’ positions and velocities do not fully capture the complexity of the Galactic mass distribution.

In our earlier work (Chakrabarti et al. 2021), we used pulsar timing to measure accelerations directly and found differences at the level of a factor of ∼2 for the observed line-of-sight acceleration of the pulsars in our sample, compared to static potentials that are based on the Jeans analysis to estimate accelerations, which may be due to out-of-equilibrium effects. The Oort limit we measured (to 3σ) is about 15% lower than that determined from the Jeans analysis (McKee et al. 2015; Schutz et al. 2018). The oblateness of the potential traced by the pulsars is significantly closer to that of a disk rather than a halo. However, the pulsar sample is small (our earlier analysis included 14 binary millisecond pulsars that were timed sufficiently precisely such that we could extract the Galactic signal) and grows slowly, and so we are prompted to further explore the development of new direct acceleration techniques.

EBs have long been amenable to precise characterization, including ≤3% measurements of their masses and radii (Torres et al. 2010). The advent of continuous, high-precision photometry from space telescopes such as Kepler and TESS has led to comparable-or-better levels of precision being achieved for an increasingly large sample of EBs (Southworth 2015, 2021; Prsa et al. 2022); in particular, Kepler's high photometric precision—≲500 parts per million, or ppm, for its long-cadence observation of EBs at magnitudes K ≲15—and long baseline permit measurements of EBs’ eclipse times to subsecond precision (e.g., Clark Cunningham et al. 2019; Hełminiak et al. 2019; Windemuth et al. 2019). The change in the eclipse midpoint time due to the Galactic acceleration grows quadratically with time. At ∼kiloparsec distances, we expect the eclipse midpoint time to have shifted by ∼0.1 s in the decade between the Kepler observations and today, i.e., to measure the Galactic signal, we require an eclipse-timing precision of ∼0.1 s from both the archival Kepler data and a light curve today.

The outline of this Letter is as follows. In Section 2, we review the various physical mechanisms that can change a binary’s orbital period, including the general-relativistic (GR) precession of an eccentric orbit, tidal decay, tidally and rotationally induced quadrupole moments, and the acceleration exerted on the binary by planetary companions, as well as the acceleration induced by the Galactic potential. Our goal here is to determine the part of the parameter space where contaminants to the Galactic signal are sufficiently minor that we can reliably measure the small shift in the eclipse midpoint time. We also analyze sources from a recent compilation of EBs with custom-extracted light curves (Windemuth et al. 2019) for which the orbital and stellar parameters were presented. In Section 3, we present the expected timing precision for the set of sources from the Windemuth et al. (2019) paper for which contaminants to the Galactic signal should not be significant, as well as an additional 70 sources from the Kepler EB Catalog (Prša et al. 2011; Kirk et al. 2016) for which we estimate sufficiently precise mideclipse times to enable measurement of the Galactic acceleration today. We also discuss a prototypical EB and calculate its simulated HST light curve and expected timing precision. We conclude in Section 4.

For an EB with binary orbital period P_b , various physical mechanisms can induce a change in the observed binary period over time, thus affecting the observed mideclipse time, t_c . These include the contribution from the Galactic gravitational potential ${\dot{P}}_{b}^{\mathrm{Gal}}$, the Shklovskii effect ${\dot{P}}_{b}^{\mathrm{Shk}}$, and the relativistic precession of an eccentric orbit ${\dot{P}}_{b}^{\mathrm{GR}};$ these first three effects also impact the measured time rate of change of the binary period for pulsars, which we analyzed in our earlier work (Chakrabarti et al. 2021). Additionally, circumbinary planets may affect the drift rate of the binary period, ${\dot{{P}_{b}}}^{\mathrm{pl}}$ ; this last term can also affect pulsar timing, but one can place limits on possible planetary companions from existing pulsar timing data for even distant planetary companions, as in Kaplan et al. (2016). For stars that behave as fluid bodies, there are several additional effects that are also important: tidal decay, ${\dot{{P}_{b}}}^{\mathrm{tidal}}$, and rotationally and tidally induced quadrupoles, ${\dot{{P}_{b}}}^{\mathrm{quad}/\mathrm{rot}}$ and ${\dot{{P}_{b}}}^{\mathrm{quad}/\mathrm{tidal}}$. As there are significant and well-known sources of uncertainties in the tidal decay formulation (Ogilvie 2014; Patra et al. 2020), our approach here is necessarily approximate. The sum of these various mechanisms will then lead to the observed time rate of change of the binary period ${\dot{P}}_{b}^{\mathrm{Obs}}$:

$$\begin{eqnarray}\begin{array}{rcl}{\dot{P}}_{b}^{\mathrm{Obs}} & = & {\dot{P}}_{b}^{\mathrm{Gal}}+{\dot{P}}_{b}^{\mathrm{Shk}}+{\dot{P}}_{b}^{\mathrm{GR}}+{\dot{{P}_{b}}}^{\mathrm{tidal}}+{\dot{{P}_{b}}}^{\mathrm{quad}/\mathrm{rot}}\\ & & +{\dot{{P}_{b}}}^{\mathrm{quad}/\mathrm{tidal}}+{\dot{{P}_{b}}}^{\mathrm{pl}}.\end{array}\end{eqnarray} \tag{ 1 }$$

The Galactic acceleration is a_Gal, which for simplicity we take to be a Gaussian centered at ∼ few - 10 cm/s/decade for stars at ∼kiloparsec distances from the Sun (Chakrabarti et al. 2020, 2021). Thus, we write the Galactic contribution to the time rate of change of the binary period as

$$\begin{eqnarray}&&{\dot{P}}_{b}^{\mathrm{Gal}}=\displaystyle \frac{{a}_{\mathrm{Gal}}}{c}{P}_{b}.\end{eqnarray} \tag{ 2 }$$

We summarize below the additional contributions to ${\dot{P}}_{b}^{\mathrm{Obs}}$, and we follow closely the notation in Rafikov (2009) in which n = 2π/P is the orbital frequency, μ the proper motion in the plane of the sky, d the distance, e the eccentricity, and a the semimajor axis. The so-called Shklovskii effect (Shklovskii 1970) arises due to the transverse motion of the binary and can be expressed as

$$\begin{eqnarray}&&{\dot{P}}_{b}^{\mathrm{Shk}}={\mu }^{2}d\displaystyle \frac{{P}_{b}}{c}.\end{eqnarray} \tag{ 3 }$$

Tidal dissipation inside the star (Ogilvie 2014) gives rise to

$$\begin{eqnarray}&&{\dot{{P}_{b}}}^{\mathrm{tidal}}=-\displaystyle \frac{27\pi }{2}\sum _{i=1,2}\displaystyle \frac{1}{{Q}_{i}^{{\prime} }}\displaystyle \frac{{M}_{j}}{{M}_{i}}{\left(\displaystyle \frac{{R}_{i}}{a}\right)}^{5},\end{eqnarray} \tag{ 4 }$$

where j ≠ i and ${Q}_{i}^{{\prime} }$ are the reduced tidal quality factors for both stars. Here, we assume that there are equal contributions from both stars and that the stars are not tidally synchronized. This is the maximum possible contribution from tidal decay because tidal decay would be suppressed if the system is synchronized. Typical values of the reduced tidal quality factor are ${Q}^{{\prime} }\sim {10}^{6}$ (Ogilvie 2014; Patra et al. 2020), although lower values (Q′ ∼ 4 × 10⁵) have also been inferred for short-period planets like WASP-43 b (Davoudi et al. 2021).

When the binary is eccentric, there are also contributions to ${\dot{P}}_{b}$ due to the apsidal precession of the orbit. The contribution due to the GR precession (Rafikov 2009) is

$$\begin{eqnarray}&&{\dot{P}}_{b}^{\mathrm{GR}}=\displaystyle \frac{36\pi e\cos \omega }{{\left(1-{e}^{2}\right)}^{1/2}{\left(1+e\sin \omega \right)}^{3}}{\left(\displaystyle \frac{{na}}{c}\right)}^{4}.\end{eqnarray} \tag{ 5 }$$

The period variation caused by apsidal precession due to tidally and rotationally induced quadrupoles is

$$\begin{eqnarray}&&\begin{array}{l}{\dot{{P}_{b}}}^{\mathrm{quad}}=\displaystyle \frac{4\pi {\left({\dot{\omega }}^{\mathrm{quad}}\right)}^{2}}{{n}^{2}}e\cos \omega \displaystyle \frac{{\left(1-{{\rm{e}}}^{2}\right)}^{3/2}}{{\left(1+\mathrm{esin}\omega \right)}^{3}},\end{array}\end{eqnarray} \tag{ 6 }$$

where for a tidally induced quadrupole (Fabrycky & Tremaine 2007; Philippov & Rafikov 2013),

$$\begin{eqnarray}&&\begin{array}{l}{\dot{\omega }}_{\mathrm{tide}}^{\mathrm{quad}}=\displaystyle \frac{15}{8}{{nf}}_{1}(e)\displaystyle \sum _{i=1,2}{k}_{2,i}\displaystyle \frac{{M}_{j}}{{M}_{i}}{\left(\displaystyle \frac{{R}_{i}}{a}\right)}^{5},\end{array}\end{eqnarray} \tag{ 7 }$$

where j ≠ i and k_2,i are the apsidal motion constants for both stars (k_2,⋆ ≈ 0.016 for the Sun, Claret 2019), and ${f}_{1}{(e)=(8+12{e}^{2}+{e}^{4})(1-{e}^{2})}^{-5}\approx 8$ for e → 0.

For the rotationally induced quadrupole, assuming that both stellar spin axes are aligned with the orbital angular momentum axis, one has

$$\begin{eqnarray}&&\begin{array}{l}{\dot{\omega }}_{\mathrm{rot}}^{\mathrm{quad}}={{nf}}_{2}(e)\displaystyle \sum _{i=1,2}{k}_{2,i}\displaystyle \frac{{M}_{1}+{M}_{2}}{{M}_{i}}{\left(\displaystyle \frac{{R}_{i}}{a}\right)}^{5}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{i}}{n}\right)}^{2},\end{array}\end{eqnarray} \tag{ 8 }$$

where ${f}_{2}{(e)=(1-{e}^{2})}^{-2}\approx 1$ as e → 0, and Ω_i are the stellar spin rates.

We define the observed line-of-sight acceleration, ${a}_{\mathrm{LOS}}^{\mathrm{Obs}}$, as

$$\begin{eqnarray}&&{a}_{\mathrm{LOS}}^{\mathrm{Obs}}=\displaystyle \frac{c{\dot{P}}_{b}^{\mathrm{Obs}}}{{P}_{b}}.\end{eqnarray} \tag{ 9 }$$

A specific physical mechanism that induces a time rate of the binary period is denoted ${\dot{P}}_{b,i}$, which then leads to a shift of Δt_c,i :

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{c,i}=\displaystyle \frac{1}{2}\displaystyle \frac{{\dot{P}}_{b,i}}{{P}_{b}}{T}^{2},\end{eqnarray} \tag{ 10 }$$

where T is the baseline of the observations covering multiple eclipses. We take this time baseline to be a decade (roughly the elapsed time between the Kepler mission and the present). The sum of these various mechanisms discussed above leads to a total observed ${\dot{P}}_{b}^{\mathrm{Obs}}$ and a total shift in the midpoint time Δt_c . To clarify which physical mechanisms are dominant in various parts of the parameter space, we plot in Figure 1 the Δt_c,i induced by an individual physical acceleration mechanism. For systems that are circularized but not tidally synchronized, only Equation (4) would apply, while for systems that are tidally synchronized but not circularized, Equations (5) and (6) would apply.

Zoom In Zoom Out Reset image size

We have focused here on eclipse-timing measurements of stellar EBs rather than transiting exoplanets because EBs’ timing precision is generally better due to their deeper eclipses and shorter ingress–egress durations. The timing precision depends linearly on the transit depth, which is proportional to the square of the radius ratio of the two stars. Additionally, the timing precision is proportional to the square root of the ingress–egress durations; these durations are also proportional to the sizes of the stars (Carter et al. 2008; Winn 2010).

2.1. The Accessible Parameter Space for Measuring Galactic Accelerations with Eclipse Timing

To measure Galactic accelerations over a decade timescale, it is necessary to be able to measure a shift in the eclipse midpoint time to about ∼0.1 s for sources that are at ∼kiloparsec distances from the Sun. We can expect the vertical Galactic acceleration to scale approximately linearly with vertical height; we gave a fitting formula for the vertical dependence in earlier work by analyzing pulsar timing observations (Chakrabarti et al. 2021). This linear dependence is also expected from earlier kinematic analysis (Holmberg & Flynn 2000). We can expect the radial component of the acceleration to scale as ${V}_{c}^{2}/R$, where V_c is the circular speed and R is the Galactocentric radius.

Using the Kepler EB Catalog (Prša et al. 2011; Kirk et al. 2016), we identified detached EBs with no significant evidence of eccentricity by requiring that the catalog’s morphology parameter morph < 0.5 and primary/secondary eclipse separation parameter 0.49 < sep < 0.51. For each EB, we pulled its long-cadence light curve and observable quantities (e.g., periods and eclipse durations) from the Kepler EB Catalog, and we used the Price & Rogers (2014) relations to estimate the uncertainty on the mideclipse times. These criteria select a sample of 70 detached EBs; for these, we estimate Δt_c ≤ 1 s, and we have inspected their light curves by-eye. The Kepler light curves of these sources typically exhibit few-hundred-ppm photometric precision. We used the batman Python package (Kreidberg 2015) to fit transit models to these data; note that we are interested solely in the mideclipse timing for our purposes and not on the accuracy of the recovered physical parameters.

Figure 2 depicts the Galactocentric coordinates for these 70 EBs, which we determine from the Gaia eDR3 data set (Lindegren et al. 2021). We also show here the set of 425 EBs from Windemuth et al. (2019) that have a low enough tidally induced quadrupole moment that the Galactic signal should be measurable. For this set of about 500 sources, about half have timing precision better than 0.5 s. Relative to our pulsar sample (Chakrabarti et al. 2021), there are a significantly larger number of EBs above the Galactic disk and more EBs above the Galactic disk at larger radial distances from the Sun. This suggests that EBs, in addition to pulsars, may allow for a tighter constraint on the Oort limit. In our earlier analysis using pulsar timing, we did not see any clear patterns in the residuals of the measured line-of-sight accelerations at the pulsar locations relative to commonly used static models. The larger number of EBs may manifest clearer residuals (that may arise due to, e.g., out-of-equilibrium effects, dark matter substructure, or a warp or a lopsided mass distribution) in the line-of-sight acceleration.

Zoom In Zoom Out Reset image size

We note that the uncertainties cited in Windemuth et al. (2019) for the specific set of sources that we consider here (which correspond to sources for which the Galactic potential produces the dominant dynamical effect) typically lead to timing precisions that are about a factor of 2–3 worse than our estimates here. If ∼0.1 s timing precision is not possible on a per-star basis, we should still be able to produce a statistical measurement by increasing our sample size correspondingly (as N²). Follow-up observational work will also require the use of TESS light curves (or other observations after Kepler) to forecast the eclipse times to better than 1 minute eclipse-timing accuracy. Data of eclipsing binaries (Petrosky et al. 2021) from all-sky surveys, such as the time-series component of the Wide-field Infrared Survey Explorer (WISE), may also help to improve the timing precision relative to Kepler observations alone. Indeed, it may be necessary in our follow-up observational work to obtain multiple transits rather than a single transit to obtain the timing precision that we are aiming for.

To demonstrate that 0.1 s eclipse timing is possible in principle with HST, we use PandExo (Batalha et al. 2017) to generate a synthetic HST light curve covering the 25% deep primary eclipse of this bright (K = 11.9; J_mag = 10.93), P_b = 23 day EB KIC 4144236. With 2 minute exposures, 28 exposures per HST orbit, and 6 orbits, our synthetic light curve has a typical per-point flux precision of 104 ppm; by fitting a batman transit model to it and estimating the timing uncertainty using a Markov Chain Monte Carlo (MCMC) analysis, we are able to determine the mideclipse time to 0.1 s. Figure 3 shows our synthetic light curve, best-fit transit model, and MCMC posterior for ${\sigma }_{{t}_{c}}$.

Zoom In Zoom Out Reset image size

Additional variability in EB light curves can in principle reduce the accuracy of eclipse-timing measurements. One can account for ellipsoidal variations and rotation signals from the Kepler light curves via, e.g., basis spline fitting (Vanderburg & Johnson 2014). These signals’ timescales are generally long relative to the eclipse durations, so they will manifest as low-order polynomial flux changes in light-curve data. Follow-up spectroscopy and radial velocities can constrain EBs’ eccentricities to <1% and further rule out third bodies that might contribute to eclipse transit viariations (ETVs). Our analysis here, sampling from the observed demographics of circumbinary planets, suggests that dynamic effects from planets would not contaminate the Galactic signal for sufficiently large samples of acceleration measurements. One can estimate (as we have done above for circumbinary planets) the contribution of a third body to ETVs. For the sample of 392 hierarchical triple stellar systems that are resolved by Gaia and have wide separations, the contribution to Δt_c,i is much smaller than ∼0.1 s. However, there is also a population of compact triple stellar systems where the magnitude of ETVs is large (Borkovits et al. 2015). Most binaries with P > 10 days do not have tertiary companions (Tokovinin et al. 2006; Tokovinin 2020). To select for sources where other dynamical effects are subdominant, our sources typically have periods larger than 10 days, and we may therefore expect this effect to be small, although any observational follow-up should vet for third light.

These new techniques—pulsar timing, time-domain optical spectroscopy using extreme-precision radial velocity observations or “optical timing,” and eclipse timing—now enable precision measurements for “real-time” Galactic dynamics. Establishing a baseline for distant EBs will prove fruitful for upcoming space-based missions such as Roman, as the signal grows quadratically with time.

We summarize our main findings below.

1. 1.  
   We show that, in principle, it is now possible to use HST (or other space missions that can achieve comparable photometric precision) of Kepler eclipsing binaries to measure the Galactic acceleration. By measuring an eclipse today, we can detect the ∼0.1 s shift in the midpoint of the eclipse time due to the Galactic gravitational potential. If ∼0.1 s timing precision is not possible on a per-star basis, we should still be able to produce a statistical measurement by increasing our sample size correspondingly (as N²).
2. 2.  
   We have analyzed contributions to the period drift rate from sources other than the Galactic potential, including the relativistic precession of an eccentric orbit, tidal decay, and the rotationally and tidally induced quadrupole. Our approach here in modeling these contributions is necessarily approximate due to the uncertainties in modeling tidal decay, and our calculation of the tidal decay contribution is an upper limit assuming the stars are not tidally synchronized. For low periods, there can be significant contributions from the GR precession and the tidally and rotationally induced quadrupole even for low eccentricities.
3. 3.  
   We calculate the contributions to the change in the eclipse midpoint time due to these dynamical effects for the sample of EBs analyzed earlier by Windemuth et al. (2019) using their orbital and stellar parameters. For this data set, assuming that the stars are tidally synchronized, there are 425 sources with Δt_c,i < the Galactic signal, such that the Galactic acceleration is indeed measurable.
4. 4.  
   We create a synthetic population of circumbinary planets by sampling from their observed demographics and calculating their contribution to Δt_c,i. We find that there is a sufficiently small overlap with the Galactic signal such that circumbinary planets are not a significant contaminant to the Galactic signal.
5. 5.  
   Using analytic relations for the timing precision, we find that there are about 200 EBs in the Kepler field that have timing precision better than 0.5 s and for which all other dynamical effects are subdominant to the Galactic signal. These EBs populate a different region of physical parameter space than the pulsars earlier analyzed by Chakrabarti et al. (2021). We may be able to leverage the size of the larger EB sample (compared to the pulsar sample) to place constraints on nonequilibrium effects and/or dark matter substructure by analyzing residuals in the line-of-sight acceleration relative to static models and to improve the precision of the Oort limit and local dark matter density relative to the earlier pulsar timing analysis.
6. 6.  
   We used the EB KIC 4144236 (K = 11.9; J_mag = 10.93) as a worked example, analyzing a simulated HST WFC3 light curve (∼100 ppm flux precision in 2 minutes of exposures) of its 25% deep primary eclipse. The light curve has a photometric precision of 104 ppm, and we are able to constrain the mideclipse time to 0.1 s.
7. 7.  
   Unlike kinematic methods, the signal from direct acceleration methods like this one grows with time and, in this case, quadratically with time. Thus, the Galactic signal from EBs for upcoming Roman observations will be substantially larger than what it is today (when the measurement has just become possible), given the baseline established a decade ago by Kepler.

The Center for Exoplanets and Habitable Worlds and the Penn State Extraterrestrial Intelligence Center are supported by the Pennsylvania State University and the Eberly College of Science.

S.C. gratefully acknowledges support from NSF AAG 2009574 and hospitality provided by the CCA at the Flatiron Institute. R.R.R. is supported by STFC grant ST/T00049X/1. We thank N. Sehgal for coining the word “optical timing.” We are grateful to S. Tremaine for helpful comments on the paper and for catching an error on the value of the apsidal motion constant. We thank J. Winn, R. Dawson, G. Ogilvie, A. Prsa, and J. Pepper for helpful discussions and the anonymous referee for helpful comments.

This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.