Handling secular period change and/or center-of-mass acceleration?

[elnjensen]

I'm currently modeling the data for a binary system for which we have some radial velocities taken in the late 1980s and TESS data (eclipses) taken during the past few years, so the two datasets are about 30 years apart. In order to fit both datasets, we have to introduce a dP/dt term; it appears that the orbital period is changing slightly. This is consistent with a modest acceleration of the center of mass of the system. So the system-rest-frame period of the orbit isn't actually changing, but for each orbit there is progressively more light-time delay of the signal, with the amount per orbit changing with time due to the acceleration.

Currently we are handling this by calculating separate Keplerian orbits for the RV data and for each of the two TESS sectors; all orbital parameter variables are shared between the orbits, except for the period, which we model as a linear function of time, fitting for a slope and an intercept. At the moment each dataset has a single orbital period (for a given value of the slope and intercept) calculated at the mean time of the datapoints in that dataset, since the KeplerianOrbit instance requires a single variable for the period. That yields a pretty good fit, but it isn't strictly correct since each time point really has (effectively) a different orbital period, at least in terms of what is observed from Earth.

Is there a better way to incorporate this into the fit, i.e. in a more continuous way? The effect is tiny over a single 26-day TESS sector, but the RV data span 3-4 years.

Thanks for any thoughts.

[dfm]

Interesting question! For this case, I think that probably the simplest approach that would do what you want would be to "warp" the times before calculating the model. A start would be something like this:

import exoplanet as xo
import numpy as np
import matplotlib.pyplot as plt

tref = 250.0
period_at_tref = 10.0
pdot = 0.005
orbit = xo.orbits.KeplerianOrbit(period=period_at_tref)

warped_t = (t - orbit.t_periastron.eval()) / (1 + pdot * (t - tref) / period_at_tref) + orbit.t_periastron.eval()

t = np.linspace(0, 500, 10000)
y0 = orbit.get_radial_velocity(t, K=1.0).eval()
y1 = orbit.get_radial_velocity(warped_t, K=1.0).eval()

fig, axes = plt.subplots(1, 3, sharey=True, figsize=(15, 5))
for ax in axes:
    ax.plot(t, y0, label="constant period")
    ax.plot(t, y1, label="changing period")
axes[0].set_xlim(0, 20)
axes[1].set_xlim(240, 260)
axes[2].set_xlim(480, 500)

axes[0].set_ylabel("rv")
axes[1].set_xlabel("time")

which generates the following figure:

The warping function comes from the mean anomaly equation:

mean_anom = 2 * pi * (t - tperi) / period
          = 2 * pi * (t - tperi) / (period + pdot * (t - tref))
          = 2 * pi * (t - tperi) / period / (1 + pdot * (t - tref) / period)
          = 2 * pi * (t - tperi) / (1 + pdot * (t - tref) / period) / period
          = 2 * pi * (t_prime - tperi) / period

where

t_prime = (t - tperi) / (1 + pdot * (t - tref) / period) + tperi

[elnjensen]

Ah, that’s an elegant solution! We’ll give it a try. I had thought about trying some kind of time transformation, but transforming to the system rest frame seemed overly complicated - this is much better.

Thanks!