1 Introduction

The Transit Time Variation (TTV) technique is a powerful tool to discover multi-bodies orbiting around a star, and can be used to characterise exoplanetary systems by measuring changes in the orbital period due to gravitational interaction between planets [1, 21, 34]. The amplitude of the TTV effect is greatly boosted if the planets are near a mean-motion resonance (MMR, [1, 21]), allowing us to infer, in particular configurations, even the presence of an exoplanet in the Earth-mass regime. The Kepler mission [6], and its extension Kepler/K2 [23], demonstrated the potentiality of this technique, able to characterise planets also for stars too faint for a spectroscopy analysis and radial velocity measurement. Some examples of the application of the TTV technique on multiple planet systems are: Kepler-9 [4, 5, 22], Kepler-11 [30, 31], Kepler-88 [39, 57], K2-24 [42, 43], and WASP-47 [54, 56].

A TTV can be the result of other different physical sources, e.g. star spots, oblateness of hosting star, orbital precession induced by tidal interactions, and general relativity (for a general review see [41]), and non-physical sources such as binning or sampling of the data (e.g. [25, 33, 45, 52]). In this work we focus on TTV only due to gravitational interaction of exoplanets.

The European Space Agency (ESA) adopted Ariel (Atmospheric Remote sensing Infrared Exoplanet Large survey; [40, 44, 46, 53]) as an M4 mission within the Cosmic Vision programme, and its launch is scheduled for 2028. The Ariel mission aims to study the atmospheres and the chemical compositions of exoplanets by surveying about 1000 transiting exoplanets through fast-cadence high-precision spectroscopy and photometry of transits, eclipses, and phase-curves in a wide spectral range from the visible (VIS, 0.5 μ m) to the infrared (IR, 7.8 μ m; detailed description in [11, 40, 53, 60]). The Ariel satellite mount two spectrometers, the Ariel Infrared Spectrometer (AIRS, two prism-dispersed channels that cover the band width 1.95–3.9μ m with R > 100 and 3.9–7.8μ m with R > 30) and NIRSpec (slit-less prism spectrometer with spectral resolving power R > 15 in spectral range of 1.1–1.95 μ m), and three photometers, namely VISPhot (1.1–1.95 μ m) and two Ariel Fine Guidance Sensors (FGS1 at 0.6–0.8 μ m and FGS2 at 0.8–1.11 μ m). See [53], [40], and [37] for further details on Ariel instruments. The fast-cadence photometric data will allow us to measure the planetary transit time at precision and accuracy level of a few seconds. This is crucial to extend the temporal baseline of TTV signals of known planets in multiple-planet systems, breaking the degeneracy on doubtful dynamical solutions and improving the precision on their masses and orbital parameters [8, 43].

In this work we present a study on the Ariel timing performances based on simulations of real targets (Section 2) and some possible science cases (Section 3). This analysis describes one of the possible extended uses of the Ariel’s Core data and it has been developed within the Ariel Working Group High-Precision Photometry (described in details in Szabo et al.; 2020, submitted and in Haswell 2020, in preparation).

2 Transit time analysis

We simulated two science cases, modelled after the real planet hosts 55 Cnc e and K2-24 b, by analysing synthetic light curves assuming the nominal performances of FGS1 and 2. In our analysis we simultaneously fitted the light curves gathered from the two channels to get a single value of the parameters that are not dependent on wavelength, such as the T0. We simulated photometric data at three different cadences: 1 s, that is the nominal cadence of Ariel observations, 30 and 60 s as a comparison, in order to understand when the sampling cadence becomes a limiting factor for our analysis.

For both targets we estimated the expected Ariel photometric noise per time unit, and for each band-pass, by using the ArielRad code [37]. This program takes into account different sources of noise, such as detector noise (dark current, gain and read-out noise), and other stationary random processes modelled as Poisson noise, i.e. photon noise, Zodiacal noise, and instrument emission). ArielRad cannot include non-stationary noise (time-correlated noise), so it includes a jitter noise as the variance at different wavelength on timescale of 1 hour, and it includes a margin noise to take into account noise uncertainties and possible time-correlated effects (about 20% on top of the photon noise, see [37]). We re-scaled the noise to the actual sampling cadence and we added a noise floor of 20 ppm in quadrature (as shown in eq. 17 of [37]). This noise floor is needed to avoid to indefinitely integrate down the noise [2, 17, 35], and it is 20 ppm at every integration time [37]. Due to the lack of time-correlated model in the noise analysis of ArielRad, we also added a linear trend (as 0.2% on the transit time scale) to the light curve to mimic the presence of systematic effects from astrophysical (e.g. stellar activity) or from instrumental sources (e.g. due to stray light or pointing effects of the satellite). We scaled the FGS2 trend as a random Gaussian with σ = 0.005). We also created light curves with a quadratic term,Footnote 1 and scaled the FGS2 as before. We analysed the statistics of Table 4 of [20] and we found that for 92.2% of the planets (on a sample of 2339) the trend could be described by a polynomial of first order; the rest can be detrendend with a second and third order polynomial (4.7% and 3.2%, respectively). So, we can assume that in most cases a linear trend is a good approximation of the out-of-transit on timescales of 2–3 transit duration. However, the measurement of the transit time and of its uncertainty is strongly affected by the sampling of the ingress and egress phases and by the symmetry of the transit model, so, the trend (and the red noise) contributes only at a second order. We decided to only fit a linear trend even if the light curve has been created with a quadratic one and a red noise term.

We modelled the transit light curve with the batman package [28] and we ran a statistical analysis with the emceeFootnote 2 package [13, 14], a Markov-Chain Monte-Carlo (MCMC) code with an affine-invariant sampler [15].

Table 1 Stellar and planetary parameters adopted in the synthetic light curves. 55 Cnc e parameters from [9, 55], K2-24 b from [42, 43]
Full size table

Firstly, we generated the synthetic transit light curve with the physical parameters reported in Table 1. Then, for each cadence, we fitted common parameters (different parameterization with respect to Table 1) between the two filters, i.e. the stellar density in solar units (ρ⋆), the base-2 logarithm of the period (\(\log _{2} P\), used in combination with ρ⋆ to constraint the a/R⋆), the impact parameterFootnote 3 (b), a combination of eccentricity e and argument of pericenter ω (\(\sqrt {e}\cos \limits \omega \), \(\sqrt {e}\sin \limits \omega \)), and the transit time (T0), while, for each filter, we fitted the planet to stellar radius ratio (k = Rp/R⋆), the quadratic limb-darkening (LD) coefficients, q1 and q2, an uninformative sampling of quadratic LD u1 and u2 in the form q1 = (u1 + u2)2 and q2 = 0.5u1(u1 + u2)− 1 as suggested in [27], a jitter factor in base-2 logarithm (\(\log _{2} \sigma _{\text {jit}}\)), and linear trend. We selected the best (and common) parameterization of the fitting model to reduce the correlation among parameters and allowing the emcee code to properly, and as evenly as possible, sample the parameter space. We wanted to simulate an Ariel transit observation, i.e. different transit depth for each different band-pass depending on different planetary radius that corresponds to different layers of atmosphere. Given that the two FGS channels have different, almost non-overlapping, band-pass (see Figure 2 in [37]), we decided to fit different k, but we generated the synthetic light curves with a single value. We used as priors the stellar density (computed from stellar mass and radius), period P , and eccentricity reported in Table 1, and we set very tight prior if the parameter has a null error (i.e. eccentricity of 55 Cnc e). We did not apply any priors on the LD coefficients.

The emcee code outputs the posterior distribution of the fitted parameters from which we extract a best-fit solution and an uncertainty as the high density interval (HDI at 68.27%, equivalent of 1-σ) for the transit model parameters, including the most relevant here, i.e. the transit time (T0) and its uncertainty (\(\sigma _{T_{0}}\)). As best-fit solution we compared the median, the mode,Footnote 4 and the maximum likelihood estimation (MLE, as the parameter sample that maximises the likelihood within the HDI). In our cases the differences among these best-fit parameter sets are negligible (in terms of values and uncertainties) and we decided to use as representative of the best-fit solution the median of the posterior distribution, because of its symmetry with respect to the HDI. We want also to stress that finding the best-fit solution was out of the purpose of this work and we wanted to focus on the determination of the uncertainty of the parameters, mainly the \(\sigma _{T_{0}}\).

2.1 55 Cnc e

The first case, 55 Cnc e, has been selected to assess the performances of Ariel on transit timing at the bright end. The magnitude and spectral type of the host star 55 Cnc (V = 5.95 mag, M⋆ = 0.91M⊙, R⋆ = 0.94R⊙; [55]) and the size of planet e, \(R_{\mathrm {e}} \sim 2 R_{\oplus }\) [9] represent a perfect case to assess the feasibility of a solar-like star with a small transiting planet. Its transit has very short, and almost vertical, ingress and egress phases, reducing the impact of the LD effect to the timing of the transit model.

We obtained stellar and planetary parameters from [55] and [9] and available through the NASA Exoplanet ArchiveFootnote 5. We computed the coefficients of a quadratic LD-law (u1 and u2) fitting the ATLAS [29] and PHOENIX [24] models through the toolFootnote 6 developed by [12]. We assumed a box-shaped filter for both FGS1 and FGS2 within the spectral range of 0.8 − 1.0μ m and 1.05 − 1.2μ m, respectively. We adopted the values of u1 and u2 for the fitting as the average of the values computed for each combination of stellar parameters and their extremal errors for both models, obtaining: u1 = 0.33,u2 = 0.23 for the FGS1 and u1 = 0.26,u2 = 0.25 for the FGS2. We also tested LD coefficients for band I and J computed with exofastFootnote 7 [10]: u1 = 0.39,u2 = 0.22 for the I band and u1 = 0.22,u2 = 0.29 for the J band. We used the LD coefficients of I for FGS1 and those of J for FGS2 and we found no difference with the results obtained from the LD coefficients computed from the stellar models. See in Table 1 the summary of the stellar and planetary parameters used to create the synthetic light curves.

The analysis (see best-fit model for each cadence in Fig. 1) yields an estimated uncertainty in the transit time (\(\sigma _{T_{0}}\)) of 12 s, 12 s, and 14 s for cadence at 1 s, 30 s, and 60 s respectively, that is the sampling cadence for bright targets is not a limiting factor.

Fig. 1
figure 1

Transit light curves of 55 Cnc e with simulated Ariel noise (blue dots with gray error bars) and quadratic trend (orange line) along with the best-fit transit model with fitted linear trend (green line) for FGS1 (left) and FGS2 (right). First row for 1 s cadence, second and third row for 30 s and 60 s, respectively

Full size image

2.2 K2-24 b

The second case, K2-24 b, is a Neptune-size transiting planet in a multiple-planet system showing anti-correlated TTVs with another transiting planet (K2-24 c). The star is much fainter than 55 Cnc (V = 11.28, K = 9.18) with a mass and radius slightly bigger than the Sun (M⋆ = 1.07M⊙, R⋆ = 1.16R⊙).

We studied the transit timing performances of Ariel with the same two FGS channels as before and we also analysed the impact of a few Ariel observations on the recovering of the planetary ephemeris and on the improvement of the other orbital parameters (see Section 3). We used the stellar parameters from [43] and we computed the quadratic LD coefficients, as we did for 55 Cnc e. We found u1 = 0.29,u2 = 0.25 for the FGS1 and u1 = 0.23,u2 = 0.26 for the FGS2. We also tested LD coefficients from exofast: u1 = 0.32,u2 = 0.27 for the I band and u1 = 0.17,u2 = 0.30 for the J band. A summary of the stellar and planetary parameters are reported in Table 1.

In the K2-24 b case (Fig. 2), when using LD for I and J we obtained for each cadence the same \(\sigma _{T_{0}}\) of 35 s. This result suggests us that in this situation we are photon noise dominated. In the analysis with the LD from box-shaped FGS1 and 2 filters we obtained \(\sigma _{T_{0}}\) of 34 s, 29 s, and 32 s for cadence at 1 s, 30 s, and 60 s respectively. In both analysis the \(\sigma _{T_{0}}\), obtained with only one transit, is about half the best timing uncertainty (\(\sim 78\) s) obtained with Kepler/K2 observations [43].

Fig. 2
figure 2

Same as in Fig. 1, but for K2-24 b

Full size image

We also analysed the case of K2-24 b observed by Ariel with only one channel, i.e. the FGS1. This could mimic the possible failure of one of the instruments, reducing the number of light curves per transit. We generated the synthetic light curve with the parameters in Table 1 with LD coefficients of FGS1 in I band. We repeated the analysis as before and we obtained a \(\sigma _{T_{0}}\) of about 53 seconds, that is still better than Kepler/K2 timing obtained from the simultaneous analysis of four transits for this target. This could be due to the high cadence photometry that is able to homogeneously sample the ingress and egress, the transit phases that have the greatest impact on determining the transit time.

3 Dynamical analysis

3.1 Impact of Ariel observations

K2-24 is a multiple-planet system showing anti-correlated TTVs by the transiting planet b and c, but with a poor sampling of the TTV period on K2 data [42, 43]. This can lead to an imprecise or ambiguous determination of masses and orbital parameters of the planets (see the Kepler-9 case by [4, 5, 22]). [43] found a possible candidate planet d, from the radial velocity (RV) analysis, with a mass of 54 ± 14M⊕ and a period of about 427 days.

We simulates the K2-24 system with the dynamical integrator within the TRADES codeFootnote 8 [4, 5, 32]. We integrated the orbits for the same time range spanning the observations by Kepler/K2 and follow-up (radial velocities, RVs, and additional transits) and computed all the transit times (T0s) and RVs within the same range. We then selected the corresponding T0s and RVs analysed by [43] and we added white noise, based on the actual error bars from the same work. We used TRADES with the emcee module to get masses and orbital parameters from this synthetic data set and to compute the so-called “Observed - Calculated” (O − C) plot, where the residuals of the observed T0 with respect to the reference linear ephemeris are shown as a function of time (see Fig. 3).

Fig. 3
figure 3

Top panel: Observed - Calculated (O − C) plot of K2-24 b (left) and c (right) synthetic transit times for the baseline of the observations. The Observed is given by the synthetic T0s (open-black dots) or simulated T0s (filled-blue dots), while Calculated is given by the T0s computed from linear ephemeris (from synthetic T0s). The gray lines are 100 drawn from the posterior distribution. Bottom panel: residuals between observed and simulated O − C. The x-axis is in unit of days, starting from a reference time BJDTDB − 2454833 = 2071.0 to match the transit times from [43]

Full size image

If we integrate for 10 years further the orbital solution obtained, we have two main effects on the prediction of the transit times. The accumulating errors of the linear ephemeris produce a shift of the transit time with respect to the predicted one. Another source of error on the transit time is due to the uncertainty on the orbital parameters of the planetary system (see. Fig. 4). In 2028, when the Ariel launch is supposed to take place, the shift from the linear ephemeris could lead to almost one day of uncertainty on the transit times. Even if we could recover the ephemeris, we still have more than three hours of uncertainty on transit times due to the errors on the orbital parameters.

Fig. 4
figure 4

Same as Fig. 3, but extending the orbital integration to almost 10 years in the future. The linear trend of the departure from O − C = 0 is due to the short time coverage of the linear ephemeris computed from the synthetic T0s

Full size image

We simulated the K2-24 system and ran a dynamical retrieval with TRADES+emcee in the case of our synthetic data set with two observations with Ariel, one for planet b and one for c. We assumed the same error bar for the two transit times \(\sigma _{T_{0}} = 35\) s, that is the highest (conservative) value obtained from the previous analysis (see Section 2.2). We built the O − C plots (see Fig. 5) from the best-fit solution extracted from the posterior distribution and we found that just one transit per planet is enough to recover the ephemeris and to greatly improve the knowledge of orbital parameters of the system. This can be easily achieved taking advantage of existing ground-based transit survey, such as the The Asiago Search for Transit timing variations of Exoplanets project (TASTE, [16, 38]), the Next-Generation Transit Survey (NGTS, [58]), and the ExoClock ProjectFootnote 9 (within the Ariel Ephemerides Working Group), or space missions, i.e. Transiting Exoplanet Survey Satellite (TESS, [48]), CHaracterizing ExOPlanets Satellite (CHEOPS, [3, 7]), PLAnetary Transits and Oscillations of stars (PLATO, [47]). Thanks to TESS and to the extended operations TESS will also be of great advantage for the ephemeris recovery of almost all the Kepler/K2 targets on the ecliptic and most of the targets on the sky. With only two transit times from Ariel we get an improvement on the uncertainty of the mass of about 22% and 31% for planet b and c, respectively.

Fig. 5
figure 5

Same as Fig. 4, but for the analysis with two T0s by Ariel, one for planet b (left) and one for planet c (right). The thin gray area is the result of 100 random draws from the posterior distribution; this area is thinner than in Fig. 4 as the result of better determination of orbital parameters

Full size image

3.2 Possible TTV signal from Ariel observations

In the most favourable case Ariel will be able to observe 10 transits of the same planet [46, 53]. We wanted to determine the amplitude of the TTV signal (ATTV) produced by an external perturber. We decided to consider as transiting planet the K2-24 b (orbital parameters from [42, 43]) and we used TRADES, by choosing the mass and period of an hypothetical outer perturber from a grid with 30 log-spaced values of mass, ranging from 1M⊕ to 4MNep, and 30 log-spaced values of orbital period, from 25 to 100 days. For each simulation the code integrates the orbits for 4 years (as the nominal duration of the Ariel mission), selects the transit times, with associated error of 35 s determined in Section 2.2, and it computes the linear ephemeris. Then it selects 10 random transits (without replacement), re-computes the linear ephemeris, and calculates the ATTV as the semi-amplitude of the O − C (selected transit times minus the newly computed linear ephemeris). It repeats this for 100 times and computes the median of ATTV. We obtained a map of the ATTV as a function of the mass (Mperturber) and of the period (Pperturber). It is well known that the eccentricity of the perturber (eperturber) boosts the ATTV, so we did the same analysis three times, with three different initial values of eperturber: 0.0, 0.05, and 0.1 (see Figs. 6, 7, and 8).

Fig. 6
figure 6

TTV amplitude (ATTV) for K2-24 b as a function of mass and period of the perturber. Each gray dot is the combination of mass-period in a grid of 900 simulations, where Mperturber = [1M⊕,4MNep] and Pperturber = [25,100] days, both in 30 logarithmic steps. For each simulation we selected 10 transit times within 4 years of orbital integration, we repeated the selection 100 times and computed the ATTV as the the median value of the semi-amplitude of the O − C with respect to a linear ephemeris. In these simulations the initial eperturber was set to 0.0. The white contour lines are the RV semi-amplitude due to the perturber

Full size image
Fig. 7
figure 7

Same as Fig. 6, but with initial eperturber = 0.05. It is evident that a perturber with non-zero eccentricity extend the possible TTV signal regions and increase their amplitudes

Full size image
Fig. 8
figure 8

Same as Fig. 6, but with initial eperturber = 0.1. In this case, the effect of the eccentricity of the perturber is even stronger than that in Fig. 7

Full size image

The darkest regions of the P-M parameter space (shown in Figs. 5, 6, and 7) where the amplitude ATTV is much larger than the timing error \(\sigma _{T_{0}}\) for the individual transit observations (\(\sim 35\) s, Section 2.2) are those where the TTV signal could be detected at high confidence with Ariel. At least in some orbital configurations close to a MMR, and especially for eccentric orbits, it will be possible for Ariel to robustly detect a TTV signal induced by an outer pertuber with one Earth mass (or greater) with just 10 transit observations.

4 Conclusions

According to our simulations, the Ariel photometric and timing precision is suitable to study known multiple-planet systems showing TTV signals in order to extend the TTV phase coverage, breaking the degeneracy on the retrieved orbital parameters, and in particular to improve the mass estimation. The latter is fundamental to put strong constraints on planetary formation, migration, and evolution processes (e.g. [36]). However, to avoid losing transits within the Ariel observing window due to lack of precision on the ephemeris and/or on the orbital parameters, it will be necessary to constrain the ephemeris by taking advantage of ground- (e.g. ExoClock Project, TASTE, and NGTS) and space-based telescopes (e.g. TESS, CHEOPS, and PLATO).

We warn the reader that even if ArielRad takes into account different sources of stationary noise and includes a noise margin, it cannot model time-correlated noise. So, we can consider the values of the uncertainty on the T0 obtained in our analysis as optimistic. As soon as the on-flight performances of the satellite and of instruments will be precisely evaluated we will be able to perform a deeper and more robust analysis of the impact of correlated noise on the determination of the transit time and of its uncertainty. A possible way to model correlated noise (or coloured noise, such as pink and red noise) would be to use non-parametric models such as ARMA or ARIMA processes, creating, for example, a common coloured noise for both FGS1 and FGS2 and an additional coloured noise different for each filter. We want to stress that currently the instrumental correlated noise is unknown and unpredictable and also the stellar noise can vary star-by-star, so it is not so trivial to properly model these noise components.

In the future it would be interesting to repeat the analysis taking into account multiple transits of the same planet observed with Ariel (e.g. Tier 3 targets, about \(\sim 10\%\) of Tier 1) and assess the improvement on \(\sigma _{T_{0}}\), even if the main contribution derives from the sampling rate of the transit ingress and egress. Furthermore, an MCMC analysis with a different algorithm, such as Nested and Dynamic Nested Sampling [19, 49, 50] implemented in dynesty [51] would help to sample in a more efficient and robust way a high dimensional space as in the modelling of many transit light curves in different bands and modelling the dynamics of multi-planet systems.

Furthermore, Ariel will be able for some orbital configurations to independently detect TTV signals with about 10 transit observations, allowing us to identify perturbing planets with masses spanning the Earth-Neptune range. Combining RV data with transit times helps to put strong constraints on bulk density and orbital parameters of the transiting planet (e.g. Kepler-19 and Kepler-9, , respectively [5, 32]). However, if the RV data have not enough precision to detect the external perturber, and taking into account the precision in the mass measurement of 17% for K2-24 c in [43], we can say that “in any case” the detection of a TTV signal in Ariel data will allow us to determine planetary masses of perturber with a precision better than 20% in the Earth-Neptune regime.

The Ariel target list is evolving with time and with new discoveries, so at the moment is very difficult to give a quantitative estimates of the number of known multi-planet system with TTV that Ariel will observe and characterise. We can expect about several tens of planets with significant TTV improvements, and a similar number of targets with TTV, whose contribution from Ariel will be important. [18] predicted that about less than 4000 planet systems will be discovered during the nominal mission of TESS and about half of these will be in multi-planet systems. 30 of these multi-planet systems will have a measurable TTV signal and only 10 will be characterised with precise mass measurement from TTV only. Differently from missions like TESS or PLATO, Ariel will point a single target, as CHEOPS, and it will observe about 1000 targets, roughly a quarter of the TESS predicted value. So, we can assume for Ariel a similar number of the predicted TTV systems for TESS, but only before launch we will know which planets will be observed and how many of these will be systems with TTV, and it will mostly depend on the selection process and on the schedulability of the targets.