V1298 Tau with TESS: Updated Ephemerides, Radii, and Period Constraints from a Second Transit of V1298 Tau e

Planetary radii are expected to evolve over time due to a variety of endogenous and exogenous physical processes such as gravitational contraction, atmospheric heating and mass loss, and core-envelope interactions (e.g., Lopez & Fortney 2013; Owen & Wu 2013; Jin et al. 2014; Chen & Rogers 2016; Ginzburg et al. 2018). The most dramatic changes are believed to occur at early stages when planets are still contracting and radiating away the energy from their formation, and when host stars are heating planetary atmospheres with high levels of X-ray and ultraviolet radiation. Since the size evolution of any individual planet is believed to be slow relative to typical observational baselines, the best way to make inferences about the size evolution of exoplanets is by measuring the sizes of large numbers of planets across a range of ages.

NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) has made significant inroads toward this objective. TESS's observations of ∼90% of the sky have allowed for exoplanet transit searches around stars ranging from the pre-main sequence to the giant branch. It is through targeted surveys of young stars such as the THYME (e.g., Newton et al. 2019), PATHOS (e.g., Nardiello et al. 2020), and CDIPS (e.g., Bouma et al. 2020) surveys, along with case studies of individual systems (e.g., Benatti et al. 2019; Plavchan et al. 2020; Hedges et al. 2021; Zhou et al. 2021) that the timeline for planetary radii evolution can be pieced together.

The V1298 Tau planetary system is one particularly valuable benchmark for understanding the size evolution of exoplanets. V1298 Tau is a pre-main sequence, approximately solar-mass star that was observed in 2015 by NASA's K2 mission (Howell et al. 2014). Analysis of the K2 data revealed the presence of four transiting planets, all with sizes between that of Neptune and Jupiter (David et al. 2019a, 2019b). There are no other known examples of exoplanetary systems with so many planets larger than Neptune interior to 0.5 au, despite the high completeness of the Kepler survey to large (>5 R_⊕), close-in planets. This observation raises the possibility of a causal connection between the extreme youth of V1298 Tau and the uncommonly large sizes of its planets.

The youth of V1298 Tau was initially established on the basis of its strong X-ray emission (Wichmann et al. 1996), high photospheric lithium abundance (Wichmann et al. 2000), and proper motion measurements (Frink et al. 1997). An additional recent blind search for co-moving stars using Gaia DR1 astrometry data found V1298 Tau was co-moving with 8 other stars (Group 29 in Oh et al. 2017). Luhman (2018) conducted a kinematic study of the Taurus star-forming region using Gaia DR2 and found new members of this group. With this new sample, they derived an age of ∼40 Myr. However, more recent analyses based on Gaia EDR3 astrometry suggests V1298 Tau may belong to either the D2 or D3 subgroups of Taurus, both of which have estimated ages ≲10 Myr (Gaidos et al. 2022; Krolikowski et al. 2021). Other studies focused specifically on the V1298 Tau system have estimated its age to be 23 ± 4 Myr from comparison with empirical and theoretical isochrones (David et al. 2019a), or 28 ± 4 Myr from isochrone fitting to the Luhman (2018) Group 29 membership list given Gaia EDR3 data (Johnson et al. 2021). While the precise age of V1298 Tau remains uncertain, most estimates fall in the 10–40 Myr range and we adopt t ≈ 20–30 Myr.

Given the system's youth and potential to reveal information about the initial conditions of close-in planetary systems (e.g., Owen 2020; Poppenhaeger et al. 2021), V1298 Tau has been the target for extensive follow-up observations. These include efforts to constrain planet masses with radial velocities (Beichman et al. 2019; Suárez Mascareño et al. 2021), measure the spin–orbit alignments of planet c (Feinstein et al. 2021) and planet b (Gaidos et al. 2022; Johnson et al. 2021), measure or constrain atmospheric mass-loss rates for the innermost planets (Schlawin et al. 2021; Vissapragada et al. 2021) and an approved program to study the planetary atmospheres using the James Webb Space Telescope (Desert et al. 2021).

Here we report on newly acquired TESS observations of V1298 Tau which help to refine the orbital ephemerides of the transiting planets and enable comparison of the planet sizes inferred from two different telescopes with different bandpasses (TESS and Kepler). The goal of this Letter is to provide a quick analysis of the new TESS data to help improve the transit timings for follow-up observations being performed by the community. We describe the observations and light-curve extraction in Section 2. In Section 3, we present our light-curve modeling and method for computing the marginalized posterior probability of a transiting planet's period from two transits observed with a large time gap. In Section 4, we discuss the differences in measured transit parameters between K2 and TESS data. We conclude in Section 5.

During its Extended Mission Cycle 4, TESS is re-observing many of the previous K2 fields. V1298 Tau (TIC 15756231) was observed by TESS in Sectors 43 (UT 2021 September 16 to UT 2021 October 12) and 44 (UT 2021 October 12 to UT 2021 November 6). For Sectors 43 and 44, we used the 2 minute light curve created by the Science Processing Operations Center pipeline (Jenkins et al. 2016) and binned the data down to a 10 minute cadence.

We compared these new light curves to our original FFI light curves. We created our initial light curves from the tica-processed FFIs by modeling the point-spread function (PSF) of V1298 Tau and the two nearby bright sources (white x's in Figure 1), following the PSF-modeling routine in Feinstein et al. (2019). ⁹ In summary, we calculated and maximized the likelihood value of seven parameters per each Gaussian: the x and y width, 2D position, amplitude, a rotational term, and a background term. The Gaussian fits are allowed to vary at each time step. Aperture photometry (example square aperture shown in Figure 1) provided a light curve with more systematics and scatter. We found that modeling the three brightest stars simultaneously, including V1298 Tau, with a 2D Gaussian created the least contaminated light curve. Our extracted light curve is shown in the top row of Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We simultaneously modeled the transits of V1298 Tau bcde and the stellar variability using the open-source packages exoplanet (Foreman-Mackey et al. 2021a, 2021b) and PyMC3 (Salvatier et al. 2016). Transit timings were originally identified using updated ephemerides from Spitzer (J. H. Livingston et al. 2022, in preparation) and we account for potential transit timing variations (TTVs) in our model using exoplanet's TTVOrbit function. This functionality allows us to fit the individual transit time per each transit, while all other orbital and planetary parameters remain the same.

exoplanet.TTVOrbit worked well when fitting V1298 Tau bcd, due to there being multiple transits. For V1298 Tau e, we found the Gaussian Process (GP) model with TTVs optimized around hyperparameters that accommodated transits where there were none; as such, we used a model without accounting for TTVs to fit the parameters for V1298 Tau e. All other transit parameters (presented in Table 1) were initialized using values from David et al. (2019b). We assumed a quadratic limb-darkening law, following the reparameterization described by Kipping (2013); this method allows for an efficient and uniform sampling of limb-darkened transit models.

Table 1. V1298 Tau Light-curve Fitting Results and Predicted Ground-based Transit Midpoint Events for V1298 Tau e

Star	Value	Prior
R⋆[R⊙]	${1.33}_{-0.03}^{+0.04}$	${ \mathcal G }(1.305,0.07)$
M⋆[M⊙]	${1.095}_{-0.047}^{+0.049}$	${ \mathcal G }(1.10,0.05)$
u1	${0.32}_{-0.19}^{+0.20}$	${ \mathcal U }[0,1]$ in q1
u2	${0.16}_{-0.29}^{+0.31}$	${ \mathcal U }[0,1]$ in q2
Prot [days]	${2.97}_{-0.04}^{+0.03}$	${ \mathcal G }(\mathrm{ln}2.87,2)$
ln(Q0)	${0.72}_{-0.21}^{+0.24}$	${ \mathcal H }(\sigma =2)$
ΔQ0	4.09 ± 1.01	${ \mathcal G }(0,2)$
f [ppt]	${0.85}_{-0.19}^{+0.11}$	${ \mathcal U }[0.1,1]$
Light Curve	Value	Prior
μ	$-{1.66}_{-9.28}^{+9.42}$	${ \mathcal G }(0,10)$
ln(σl )	$-{2.499}_{-4.463}^{+4.444}$	${ \mathcal G }(\mathrm{ln}(0.1{\sigma }_{\mathrm{flux}}),10)$
ln(σj )	$-{1.36}_{-5.06}^{+4.96}$	${ \mathcal G }(\mathrm{ln}(0.1{\sigma }_{\mathrm{flux}}),10)$
Planets	c	d	b	e
T0 [BKJD–2454833]	${4648.16636}_{-0.00339}^{+0.00269}$	${4645.41494}_{-0.00157}^{+0.00172}$	${4648.09023}_{-0.00132}^{+0.00129}$	${4648.79668}_{-0.00114}^{+0.00121}$
P [days]	${8.2438}_{-0.0020}^{+0.0024}$	${12.3960}_{-0.0020}^{+0.0019}$	${24.1315}_{-0.0034}^{+0.0033}$	44.1699
Rp /R⋆	0.0337 ± 0.0009	${0.0409}_{-0.0015}^{+0.0014}$	0.0636 ± 0.0018	${0.0664}_{-0.0021}^{+0.0025}$
Impact parameter, b	${0.14}_{-0.10}^{+0.14}$	${0.19}_{-0.13}^{+0.12}$	${0.45}_{-0.04}^{+0.05}$	${0.48}_{-0.07}^{+0.06}$
T14 [hr]	${4.66}_{-0.43}^{+0.49}$	${5.59}_{-0.53}^{+0.57}$	${6.42}_{-0.61}^{+0.66}$	${7.44}_{-0.71}^{+0.79}$
Rp [R⊕]	5.05 ± 0.14	6.13 ± 0.28	9.53 ± 0.32	9.94 ± 0.39
Rp [RJ]	0.45 ± 0.01	0.55 ± 0.03	0.85 ± 0.03	0.89 ± 0.04
TTVs [minutes]	−0.41 ± 25.38	−0.12 ± 4.08	⋯	⋯
Priors	c	d	b	e
T0 [BKJD–2454833]	${ \mathcal G }(4648.53,0.1)$	${ \mathcal G }(4645.4,0.1)$	${ \mathcal G }(4648.1,0.1)$	${ \mathcal G }(4648.8,0.1)$
log(P [days])	${ \mathcal G }(\mathrm{ln}8.25,1)$	${ \mathcal G }(\mathrm{ln}12.40,1)$	${ \mathcal G }(\mathrm{ln}24.14,1)$	${ \mathcal G }(\mathrm{ln}36.70,1)$
log(depth [ppt])	${ \mathcal G }(\mathrm{ln}1.45,0.1)$	${ \mathcal G }(\mathrm{ln}1.90,0.1)$	${ \mathcal G }(\mathrm{ln}4.90,0.1)$	${ \mathcal G }(\mathrm{ln}3.73,0.1)$
Impact parameter, b	${ \mathcal U }[0,1]$	${ \mathcal U }[0,1]$	${ \mathcal U }[0,1]$	${ \mathcal U }[0,1]$
T14 [days]	${ \mathcal G }(\mathrm{ln}0.19,1)$	${ \mathcal G }(\mathrm{ln}0.23,1)$	${ \mathcal G }(\mathrm{ln}0.27,1)$	${ \mathcal G }(\mathrm{ln}0.31,1)$
TTVs [days]	${ \mathcal G }({T}_{0,c},0.1)$	${ \mathcal G }({T}_{0,d},0.1)$	${ \mathcal G }({T}_{0,b},0.1)$	⋯

Note. u₁ and u₂ are the limb-darkening parameters sampled following the reparameterization described by Kipping (2013); P_rot is the rotation period of V1298 Tau; ln(Q₀) is the quality factor for the secondary oscillation used to fit the stellar variability; ΔQ₀ is the difference between the quality factors of the first and second modes; f is the fractional amplitude of the second mode compared the first; μ is the mean of the light curve; σ_i and σ_j are the light-curve noise and in-transit jitter. Priors are noted for parameters that were directly sampled. The distributions are as follows—${ \mathcal G }$: Gaussian; ${ \mathcal H }$: half-normal; ${ \mathcal U }$: uniform. σ_flux is the standard deviation of the light curve.

Download table as:  ASCIITypeset image

Since the tica FFIs do not provide an error estimate, we fit for flux errors within our GP model. We define the flux error as

$\begin{eqnarray}&&{\sigma }_{y}={e}^{\mathrm{ln}({\sigma }_{l})}+{y}^{2}{e}^{2\mathrm{ln}({\sigma }_{j})}\end{eqnarray} \tag{ 1 }$

where y is the normalized flux array about zero, and σ_l and σ_j are used to define the light-curve noise and in-transit jitter, which is designed to capture the added noise produced by starspot crossing events. σ_l and σ_j are also used as the first and second terms in our rotation model, which we defined as a stochastically driven, damped harmonic oscillator, defined by the SHOTerm in celerite2 (Foreman-Mackey et al. 2017).

We modeled the background within our GP. We defined a quadratic trend with respect to time for varying the background flux, where each polynomial coefficient was drawn from a normal distribution. Then, we generated a Vandermonde matrix (A) of time. This is a way of introducing a polynomial least-squares regression with respect to time. The final background flux was calculated by taking bkg = A · trend.

We performed an MCMC sampling fit to each parameter. ¹⁰ We ran 3 chains with 500 tuning steps and 5000 draws. We discarded the tuning samples from the posterior chains before calculating our best-fit parameters. Our results are presented in Table 1, along with our selected priors for each parameter we fit. These results are consistent with our original tica-processed point-spread function modeled light curves. We verified our chains converged via visual inspection and following the diagnostic provided by Geweke (1992).

We present our final GP model for stellar variability, planet transits, and best-fit transit models in Figure 2. There is a ∼1% flare at TESS BKJD ≈ 4659.18 that we do not fit.

3.1. Constraining V1298 Tau e's Period

V1298 Tau (EPIC 210818897) was observed during Campaign 4 of the K2 mission. There was a single transit of V1298 Tau e in the original K2 data, which occurred roughly in the middle of the campaign. Since no other transits were detected, this provides a lower period limit of 36 days. Additionally, there was only 1 transit of V1298 Tau e between the two TESS sectors, which provides a new lower limit of 42.7 days. Using the original transit timing from K2 and this new transit timing from TESS, we developed a new method for constraining the period of V1298 Tau e. For this analysis, we used the EVEREST 2.0 (Luger et al. 2018) version of the K2 light curve.

Determining the period of a planet from two transits with a significant time gap between surveys has previously been constrained by fitting for orbital periods using MCMC sampling, phase-folding all available transits on the derived transit times and periods, and computing a reduced-χ² fit to a flat line (Becker et al. 2019). Orbital periods providing a match to a flat line with a likelihood exceeding some threshold are then ruled out. In our new method, we fit transit models of many discrete periods at each step of the MCMC sampler rather than postprocessing from our posterior.

First, we de-trended a localized 1 day region around the transit midpoint of V1298 Tau e in the K2 and TESS light curves, assuming a constant transit depth and allowing the other transit parameters θ to vary. Then, we fit a discrete period model, allowing all other transit parameters to vary. We set the GP model to sample over discrete periods ranging from 38 to 56 days. We fit for θ assuming a constant transit depth between the K2 and TESS observations. We assumed there is no correlation between the other transit parameters and the period we are fitting for. For each step in our MCMC fit, we compute all possible periods

$\begin{eqnarray}&&P=\displaystyle \frac{1}{q}\left({T}_{\mathrm{mid},\mathrm{TESS}}-{T}_{\mathrm{mid},{\rm{K}}2}\right)\end{eqnarray} \tag{ 2 }$

where T_mid are the transit midpoints from K2 and TESS and q is an integer representing a specific harmonic. We assume a uniform prior, i.e., we have no prior preference for a specific harmonic. At each step of the sampling process, we compute a new light curve with different orbital periods, given by Equation (2). The log likelihood of the new light curves models are calculated as

$\begin{eqnarray}&&\mathrm{log}{{ \mathcal L }}_{q}={\left[\mathrm{log}\,p\left(X| {\theta }^{k},{q}^{k}=n\right)\right]}_{n}\end{eqnarray} \tag{ 3 }$

where X is the TESS light curve and n is the period being tested. We additionally calculate the sum of all log likelihoods for each period

$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=\mathrm{log}\,{{\rm{\Sigma }}}_{q}\,p(q)p(X| {\theta }^{k},q).\end{eqnarray} \tag{ 4 }$

The summation of all log likelihoods is used to compute the posterior likelihood for each sampled value of q. This analysis assumes a circular orbit for planet e and uses stellar density constraints via priors on the stellar mass and radius. ¹¹

We ran 3 chains with 500 tuning steps and 5000 draws. We discarded the tuning steps before our analysis. Our results are presented in Figure 3, where we plot the median period for each tested harmonic against the posterior probability of each harmonic.

Zoom In Zoom Out Reset image size

We find the most likely period of V1298 Tau e to be 44.17 days. We provide all tested periods and posterior probabilities in Table 2. Our period estimate is at a 4σ disagreement with the period measured in a potential radial-velocity signal for V1298 Tau e presented in Suárez Mascareño et al. (2021). This derived period estimate suggests that V1298 Tau e is in a near 2:1 mean motion resonance with V1298 Tau b. If the period of V1298 Tau e is confirmed to be within the presented range, this could indicate that V1298 Tau bcde are in a nearly four-planet resonant chain.

Table 2. Predicted Transit Midpoint Events for V1298 Tau e

P [days]	Posterior Prob.	Observable Dates UT
44.1699 ± 0.0001	0.071	21/12/2021 15:17	03/02/2022 19:22	19/03/2022 23:27
45.0033 ± 0.0001	0.069	23/12/2021 07:17	06/02/2022 07:22	23/03/2022 07:27
45.8687 ± 0.0001	0.066	25/12/2021 00:50	08/02/2022 21:41	26/03/2022 18:32
46.7681 ± 0.0001	0.064	26/12/2021 20:00	11/02/2022 14:26	30/03/2022 08:52
47.7035 ± 0.0001	0.062	28/12/2021 16:54	14/02/2022 09:47	03/04/2022 02:40
48.6770 ± 0.0001	0.061	30/12/2021 15:38	17/02/2022 07:53	07/04/2022 00:07
49.6911 ± 0.0001	0.059	01/01/2022 16:18	20/02/2022 08:54	11/04/2022 01:29
50.7484 ± 0.0001	0.057	03/01/2022 19:03	23/02/2022 13:01	15/04/2022 06:59
51.8516 ± 0.0001	0.055	06/01/2022 00:01	26/02/2022 20:27	19/04/2022 16:53
53.0039 ± 0.0001	0.053	08/01/2022 07:19	02/03/2022 07:25	24/04/2022 07:30
54.2085 ± 0.0001	0.051	10/01/2022 17:08	05/03/2022 22:09	29/04/2022 03:09
55.4692 ± 0.0001	0.049	13/01/2022 05:39	09/03/2022 16:55	⋯
56.7899 ± 0.0001	0.047	15/01/2022 21:03	13/03/2022 16:00	⋯
58.1750 ± 0.0001	0.045	18/01/2022 15:32	17/03/2022 19:44	⋯
59.6294 ± 0.0001	0.042	21/01/2022 13:21	22/03/2022 04:27	⋯
61.1583 ± 0.0001	0.039	24/01/2022 14:44	26/03/2022 18:32	⋯
62.7678 ± 0.0001	0.037	27/01/2022 19:59	31/03/2022 14:25	⋯

Note. Transit dates were calculated using the Transit Service Query Form on the NASA Exoplanet Archive (Akeson et al. 2013). We queried observable transits between 2021 December 3 through 2022 April 30. Dates presented in DD/MM/YYYY format. A machine-readable version of this table can be found here.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Independent ground-based monitoring of the system may be able to observe another transit of the outermost known planet in this system. Using the Transit Service Query Form on the NASA Exoplanet Archive (Akeson et al. 2013), we provide several potential transit midpoint events for all fitted periods in Table 2.

We compare the differences in transit R_p /R_⋆ between the K2 and TESS data in Figure 4. We masked regions in the light curve where transits overlapped. The residuals of the TESS light curve with our model (color) are plotted as well. For V1298 Tau bcd, the transit radii are smaller in the TESS data, while only the measured radius for V1298 Tau e is larger (Figure 4, bottom panel).

Zoom In Zoom Out Reset image size

The error bars from our MCMC fit on the radii of the planets are smaller than those provided by David et al. (2019a). We initialized our MCMC to fit the transit depths with a Gaussian distribution around the fitted values from David et al. (2019a) with a standard deviation of 0.1 (Table 1). The smaller errors could be due to the higher cadence of the TESS data (10 minutes versus 30 minutes) or due to larger spot-crossing events in K2. Larger spot-crossings would result in a greater uncertainty of the transit depth, and this is potentially evident in comparing the transit depth and shape for V1298 Tau e (Figure 4).

4.1. Radii of V1298 Tau bcd

4.2. Radius of V1298 Tau e

We present updated ephemerides for all four known planets in the V1298 Tau system. Our GP model accounts for TTVs for V1298 Tau bcd. The transit timings for V1298 Tau c. deviates from a linear ephemeris by −30 to 30 minutes, and V1298 Tau d deviates by −5 to 5 minutes. Additionally, we note the transits of V1298 Tau bcd occur 1.92 later, and 5.83 and 4.72 hr earlier than what is expected if we extrapolate forward the ephemerides from David et al. (2019a).

We detected a second transit of V1298 Tau e; this new transit time in combination with the transit observed with K2 allowed us to place tighter constraints on the period of the outermost planet. Our revised radius for planet e makes it now the largest planet in the system and extends an intriguing size–separation correlation in V1298 Tau such that planet size increases monotonically with separation.

We find the transit depths of V1298 Tau bcd as observed by TESS are shallower than those observed by K2 by 1–2σ, with the exception of V1298 Tau e, which is ∼3σ larger. While this could possibly be due to ongoing dusty outflows that make the transit depth appear shallower, it is more likely the differences are due to starspot crossing events, asymmetric starspots off the transit chord, or contamination from nearby faint stars in the same TESS pixels as V1298 Tau. Modeling potential starspot crossing events could be accomplished using the 2 minutes and 20 s cadence light curves, which will be available in the coming months.

The youth of these planets could additionally favor hosting haze-dominated atmospheres. However, without mass estimates for V1298 Tau bcde, it is difficult to determine if this is the cause of the different transit depths measured between the K2 data and our present work. Radial-velocity mass measurements are challenging for young planets due to underlying stellar activity.

Through an intensive radial-velocity campaign, Suárez Mascareño et al. (2021) presented a new mass detection for V1298 Tau b and V1298 Tau e. Our updated radius estimate for V1298 Tau b in combination with the mass estimate provided by Suárez Mascareño et al. (2021) would yield a density of 1.29 g cm⁻³, which is slightly higher than what was originally reported. Our updated radius for V1298 Tau e would yield a lower density of 2.04 g cm⁻³, which is still within a 1σ agreement with Suárez Mascareño et al. (2021). However, the period estimate for V1298 Tau e is at a 4σ disagreement for the period estimate provided in this study and at a 2σ disagreement with the minimum period constraint provided by this new TESS data.

A more promising approach to measuring the masses would be through TTVs (Agol & Fabrycky 2018). A full analysis of system parameters and TTVs from both the K2 and TESS light curves, and additional Spitzer transit photometry will be presented in J. H. Livingston et al. (2022, in preparation). Additional transits at longer wavelengths or simultaneous multi-band photometry or spectroscopy could corroborate the potential of constraining the properties of these young atmospheres.

All data, custom light curves, transit fits, and Python routines created and used in this work are publicly available on Zenodo and GitHub (Feinstein et al. 2022).

We thank Rodrigo Luger for developing showyourwork! (Luger et al. 2021) and helping us debug this Letter. We thank Chas Beichman, Sarah Blunt, Jacob Bean, and Darryl Seligman for helpful comments on our TESS proposal (DDT 036) and thoughtful conversations. We thank our anonymous referee for their thoughtful insights which improved the quality of this manuscript. A.D.F. acknowledges support from the National Science Foundation Graduate Research Fellowship Program under grant No. (DGE-1746045).

This research has made use of NASA's Astrophysics Data System Bibliographic Services. This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration et al. 2018).

Facilities: TESS (Ricker et al. 2015), Kepler (Howell et al. 2014).

Software: exoplanet (Foreman-Mackey et al. 2021b), EVEREST 2.0 (Luger et al. 2018), lightkurve (Lightkurve Collaboration et al. 2018), matplotlib (Hunter 2007), PyMC3 (Salvatier et al. 2016), starry (Luger et al. 2019), theano (The Theano Development Team et al. 2016), tica (Fausnaugh et al. 2020), showyourwork! (Luger et al. 2021), astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), astroquery (Ginsburg et al. 2019), numpy (Van Der Walt et al. 2011), celerite2 (Foreman-Mackey et al. 2017), batman (Kreidberg 2015) doi:10.5281/zenodo.5815452.