© The Authors 2024

1 Introduction

Multi-planetary systems account for about 25% of confirmed planets to date¹. The exoplanetary system known to host the largest number of planets so far is Kepler-90, with eight planets (Shallue & Vanderburg 2018). The vast majority of these systems are composed of small planets, super-Earths and mini-Neptunes, in the innermost part of the systems (<100 days periods). These packed systems are the ideal laboratory for testing and calibrating models of planetary formation: originating from the same original disc, the composition of these planets gives us clues about the region of the disc where they are likely to have formed. The orbital architecture of the system gives us information about the orbital evolution they may have undergone (Mishra et al. 2021, 2023; Luque et al. 2023).

In these compact systems, the gravitational interactions between planets can be strong, as predicted by Agol et al. (2005). These interactions result in variations of differing magnitude (from a few seconds to nearly half an hour) in the transit times of the planets. Measuring these transit time variations (TTV, Miralda-Escude 2002; Agol et al. 2005; Holman & Murray 2005) requires high photometric accuracy and a high measurement cadence in order to accurately determine transit times as well as ingress and egress. Studies of such systems can also benefit from a long time coverage to better identify which planet is transiting at a given time and track the transit time evolution. As TTVs are gravitational effects, they can be used to measure the mass of planets from photometric data alone, providing insights into their possible composition.

To date, there are 39 systems known to host five or more planets, but only a handful with a resonant chain (Lissauer et al. 2014). This phenomenon occurs when multiple successive planets are in (or close to) mean motion resonance (MMR, Lissauer et al. 2011; Fabrycky et al. 2014), namely, the planet period ratios are close to a ratio of small integers (r + 1) : r. Among these systems, we observe chains of at least three planets only in: Kepler-90 (Cabrera et al. 2014), HD 219134 (Vogt et al. 2015), TRAPPIST-1 (Luger et al. 2017; Gillon et al. 2017), K2-138 (Christiansen et al. 2018; Acuña et al. 2022), HD 110067 (Luque et al. 2023), Kepler-80 (Shallue & Vanderburg 2018), TOI-1136 (Dai et al. 2023), TOI-178 (Leleu et al. 2021; Delrez et al. 2023), and HD 158259 (Hara et al. 2020). Such observations question formation models that predict a breaking of the chain soon after the dispersion of the protoplanetary disk via dynamical instabilities. The exact nature of these instability processes is however still debated (Goldberg & Batygin 2022). The study of these systems therefore allows us to probe the architecture of these particular systems at a time when the system is in a stable phase.

K2-138 is a prime example of a highly packed resonant system. Originally described in Christiansen et al. (2018) using K2 data, it was found to host five sub-Neptunes close to 3:2 MMR. A better characterization of the planets’ masses was obtained by Lopez et al. (2019) using HARPS observations. Using Spitzer, Hardegree-Ullman et al. (2021, hereafter HU21) confirmed an outer sixth planet, not in the resonant chain. Finally, Acuña et al. (2022, hereafter A22) investigated the composition of the planets and found that the inner planets present an increasing water content with distance from the host star, while the outer planets have an approximately constant water content.

The cadence of the K2 photometry (8–10 min) did not allow for a TTV analysis of the system for which TTV amplitudes between 2.1 and 6.5 min are expected (Lopez et al. 2019). We therefore observed the K2-138 system with CHEOPS (Benz et al. 2021) to obtain high-precision photometry of the target and further explore the presence of TTVs. The higher cadence of CHEOPS, 60 sec in the case of K2-138 should therefore be sufficient to evaluate TTVs of even 2 min.

The paper is organized as follows. First, we describe the observations used in this work in Sect. 2 and the CHEOPS data modeling in Sect. 3. We then present improved ephemerides of the system in Sect. 4 and our TTV analysis in Sect. 5. Finally, we present our conclusions in Sect. 6.

2 Observations and modeling

2.1 Stellar and planetary properties

A22 conducted a full spectroscopic analysis using stellar modeling to derive the atmospheric and fundamental parameters of the host star. The improved host star’s parameters helped determining the possible composition of the planets in the system. To use completely homogeneous priors during our analysis, we used the stellar and planetary parameters derived in that work. For sake of clarity, the stellar parameters are given in Table 1. We note that the values reported by A22 for the planets are consistent with those in Lopez et al. (2019).

The system is thought to be made up of six planets: b, c, d, e, f, and 𝑔. Planets b through f are found to be very close the 3:2 resonance, thus forming one of the longest known chains. The periods of these planets are, respectively, 2.35, 3.56, 4.40, 8.26, 12.76, and 41.97 days.

2.2 Previous photometric observations

Available observations of the system prior to this work include the original K2 data from campaign 12. This photometry was acquired over a 79 day period, between December 15, 2016 and March 4, 2017, using the long cadence mode of the Kepler spacecraft (Howell et al. 2014). The Spitzer Space Telescope subsequently observed the system using the InfraRed Array Camera (IRAC) for 11 h between March 15 and 16, 2018 (DDT 13253; PI: J. L. Christiansen). These observations were targeting planet 𝑔, extending coverage to the tenth epoch since discovery. Finally, the Hale Telescope at the Palomar Observatory observed K2-138 with the Wide-field InfraRed Camera (WIRC) on two occasions, August 31, 2019 and November 4, 2019. Both observations were aimed at planet d, and observed epochs 182 and 194, respectively (Boyle et al. 2022, hereafter B22).

2.3 CHEOPS observations

K2-138 was observed by CHEOPS between 2020 July 22 and 2020 October 6 over twelve visits with durations between 8h42m and 34h44m as part of program 17 (PI: T.A. Lopez). The original schedule of observations was planned to acquire transits of planets d, e, f, and 𝑔. Due to the brightness of the star (V = 12.2) an exposure time of 60 s was used. In total, CHEOPS observations account for slightly more than five and a half days, covering fifteen transits of four planets in the system: four transits for planet d, three for planet e, six for planet f, and one for planet 𝑔. Details about each visit can be found in Table 2. For the sake of simplicity, each light curve is given an identifier based on the last three digits of the corresponding file name.

Light curves were obtained after the calibration and correction of the images by the CHEOPS automatic data reduction pipeline (DRPv13.1, Hoyer et al. 2020). This version of the DRP automatically delivers the target’s photometry for four different aperture radii: DEFAULT=25 px, RINF=22.5 px, RSUP=30 px, and OPTIMAL, which is computed for each target based the target’s flux and contamination of the field of view. For K2-138, the OPTIMAL aperture was computed to be of 24 px. To compute the latter, the DRP uses the GAIA DR2 catalog (Gaia Collaboration 2018) to simulate the position and flux of putative contaminants in the field of view of the target. For K2-138, no such contaminant is located within any of the apertures. The closest contaminant² appearing in the DR3 survey is located at ~13″, with a G band magnitude of ~20.8, much fainter than that of K2-138, at 12. We thus select the OPTIMAL aperture, which minimizes the RMS of the light curves.

To correct and detrend each CHEOPS light curves from systematic due to CHEOPS orbit or spacecraft effects, we use the pycheops package (Maxted et al. 2022). We mask points for which the photometry deviates by more than 1.5 × MAD (mean absolute deviation). We also mask points where the flux is strongly affected by cosmic ray hits during passages through the South Atlantic Anomaly or when large deviations of centroids of the Point Spread Function (PSF) of the target are measured. To correct for non-astrophysical noise sources such as background or stray light from the roll angle, we detrend the light curve using the basis vectors provided by the DRP. For each visit, we follow the recommendation of the DRP for the decorelations flags to use. The detrending vectors are then passed as free parameters when fitting the transits with a normal distribution centered on 0 (𝒩(0,1)). To detrend from stray light due nearby objects reflecting in the telescope (or “glint”) we fit a spline against the roll angle of the telescope. This is done in pycheops using the add_glint function, which is fitted using uniform distribution between 0 and 2 (𝒰(0,2)). The cutoff values, vectors, glint orders used, and the remaining number of points are detailed in Table 3.

CHEOPS visits are typically setup to catch one transit of a single planet, and we scheduled multiple visits for planets d, e, and f (Table 2). This way, we observed four transits of planet d with CHEOPS, thus extending the transit monitoring up to 250 orbits after the first epoch, and three transits of planet e. The ephemeris for planet f predicted the transit ~2.5 h too late. Combined with the 3.2 h transit duration, three visits completely miss the first half of the transit. Additionally, the egress falls in the observation gaps of the CHEOPS orbits in four visits. In total, out of the six available visits, only two egresses are perfectly caught. Because of the uncertainty on the previous ephemerides, some transits were predicted to occur in the CHEOPS visits, but were not detected. We identified those in the first column of Table 2. Notably, planets b and c were not detected at all, while only a single transit of planet d was not detected.

Given the number of planets in this system and their short orbital periods, the ephemerides from HU21 and A22 predict light curves with multiple transits. Five visits should contain two or more transits of the different planets in the system: 102 (b, d), 201 (b, c, e), 303 (d, f), 304 (b, f), 501 (b, d, f, 𝑔) (see Figs. 1–3).

2.4 Modeling of the CHEOPS light curves

We perform the CHEOPS light curves analysis with the pycheops package (Maxted et al. 2022), routinely used for CHEOPS observations analyses (e.g. Delrez et al. 2023; Hoyer et al. 2023). This is done in two steps. A first estimate of the model parameters are determined using LMFIT³ (Newville et al. 2014), which fits the transit using a least-square scheme. Then the parameters are refined using the emcee package (Foreman-Mackey et al. 2013), which explores the parameter space using an MCMC scheme. This second step allows for a better estimation of the fit, as well as the errors associated to each parameter. We use a normal distribution on planetary and stellar parameters priors based on their mean and 1σ interval values from A22, given in Tables 1 and A.1.

Unlike most CHEOPS visits that are of short duration, aiming to capture unique transits, visit 501 exhibits four transit-like features (Fig. 1). We first identified the planets most likely to be associated with each transit observed in this visit. Based on the ephemerides from A22, we fitted each of the observed planets in their respective visits to get a good estimate of their transits timing. We have used these updated ephemerides when possible (and retained those of A22 in cases where they were not available) to determine which planets were the most likely observed during this visit. We thus identified the first three transits as being due to planets d, f, and 𝑔, respectively. Planet b appeared as the most likely candidate to explain the fourth transit feature in this visit. However, as is detailed in Sect. 3.1, we finally discarded this possibility.

The pycheops package is designed so that only one transit per light curve can be analyzed. To overcome this limitation for visit 501, to fit each transit event, we mask the other unwanted events using the predicted positions of the planets. Therefore, we were able to proceed with the fit as if a single transit was present in each case. This allowed us to follow the same procedure as described above for each transit. pycheops also allows also for multi-visit modeling (e.g. Hoyer et al. 2023). Transits of planet d, f, and 𝑔 of visit 501, were included in their respective multi visit analysis. The resulting fits on visit 501, including the putative planet b, are shown with their ephemerides in Fig. 1.

For each planet, we detailed below the results obtained from the analysis of our K2-138 CHEOPS observations (Sect. 2.3). The results of the pycheops fits, solely based on CHEOPS data, are presented in Table A.1.

Fig. 1 Detгended light curve of CHEOPS visit 501. The colored regions correspond to the predicted transit time (with ±1σ errors) of planets d, f, 𝑔, and b according to the ephemerides from A22. The detrended CHEOPS data are represented by the light gray points, with their 5-min bins shown with the white dots. The black solid curve shows the fit of the four planet transits and arrows show the planets associated to each transit. Planet b is written in parentheses, as the depth of the event (~550 ppm) doesn't appear compatible with prior data (~250 ppm).

3 CHEOPS Transits analysis

3.1 Planet b

With an expected depth of around 250 ppm (see Table A.1 for details), the transit of planet b is at the edge of detectability with CHEOPS. Since the last observed transit of planet b dates back to K2 and that its current ephemeris predict that it should be visible in four of our visits, we nevertheless attempted to find it in our dataset. No visit was aimed directly at planet b, but it should instead appear on visits targeting other planets (see Table 2). We show the expected transits (based on A22’s ephemeris) in Figs. 1 and 2. In the cases of visits 201 and 501, the predicted transit time and its 1σ error is fully contained within the observations, unlike visits 102 and 304.

Since visits 102, 201 and 304 are quite short due to focusing on a single planet, we attempted to fit planet b’s transit signature in the residual of each visit (Fig. 2). We found no clear evidence of a photometric dip corresponding undoubtedly to planet b, therefore, we do not use them in the following TTV analysis. Visit 501 was aimed at capturing a transit of planet 𝑔 and is a much longer observation than a typical CHEOPS visit. We therefore expected good chances of seeing a transit of planet b in this dataset. The best case scenario obtained is the one shown in Fig. 1, where a transit-like feature is observed almost concomi-tantly to the predicted transit of the planet. However, the fitted transit values are not compatible with those of the literature. The modeled transit has indeed a duration T_14,fit = 1.12 ± 0.08 h compared to the expected $T_{14,\text{ prior }}={2.00}_{-0.11}^{+0.09}\text{\hspace{0.17em}h}$, as well as a depth of D_fit = 557 ± 160 ppm to $D_{\text{prior }}={252}_{-21}^{+23}\text{ppm}$. Because of this mismatch on the duration and on the depth, we do not consider it to be a transit of planet b. This short event could be the residuals of some systematics during the detrending.

3.2 Planet c

We expected to catch the planet c transit in visit 201, as shown in Fig. 3, together with those of planets e and b. With an estimated depth of ~500 ppm, planet c should be detectable with CHEOPS. Unfortunately, this visit was strongly affected by gaps and the only transit-like feature appears outside the expected ranges for the three planets (Fig. 2). Given the uncertainties in planets c and e ephemeris, however, it remains close to these ranges. We therefore explore the possibility that the transits of the two planets overlap, the event been masked by the numerous interruptions in the light curve that would hide clear ingress and egress of transits.

To that end, we performed a one-planet (planet e), and two-planet (planets e and c) transit-fitting modeling successively, using the juliet code (Espinoza et al. 2019). The advantage

of juliet, is that it allows the simultaneous fit of two planetary transits and also a straightforward comparison of the results based on Bayesian criteria. For these fits (apart from the transit mid-time), we fixed all the parameters to the central values of planets c and e priors reported in Table A.1. For the transit mid-time of planet e we used a normal distribution defined by BJD−2450000 = 9078.98 ± 0.3. For the mid-time of planet c we defined an uniform distribution between BJD−2450000 = 9078.905−9079.065, which was the wider distribution possible that secured the convergence of the fit. In Fig. 4, we present both solutions and the corresponding mid-times are given in Table A.1. As expected, the mid-times of planet e obtained from both models are consistent. The fitted mid-time of planet c ranges along the first half of the transit of planet e, namely, between the ephemeris’ predicted time (as marked by the red vertical region in right panel of Fig. 4) up to the middle of planet e transit. To determine which is the best model, juliet provides the Bayesian evidences (ln(Z)) of each fitted scenarii (see Table 4). In this comparison, the planet-e only model yields a log-evidence of ln(Z) = 1208.59 ± 0.46, while that of the planets-e and c model is of ln(Z) = 1206.37 ± 0.44, leading to a Δ ln(Z) = 2.2 ± 0.7. The limit between weak and moderate evidence of preference being ln(Z) = 2, we can therefore not statistically prefer one model over the other (Espinoza et al. 2019). This situation is discussed further when we compute the new ephemerides in in Sect. 4.

Fig. 2 Putative transits of planet b. From top to bottom: Light curves of CHEOPS visits 102, 201 and 304, in which planet b transits were expected to occur. Colored regions show the expected transits mid-time with their ±1σ errors. Transits due to other planets than b have already been removed. The solid black curve represent the best fits for planet b, while fixing its depth and duration to D = 252 ppm and T14 = 2.00 h.

Fig. 3 Detrended light curve of CHEOPS visit 201. The solid black curve represents the fitted transit of planet e. The predicted mid-transit time of planet c is shown within ±1σ error as colored regions from two different sources: A22 in green and HU21 in yellow.

3.3 Planet d

Transits of planet d are present in visits 101, 102, 103, and 501. The ephemeris for this planet from A22 predicts the transit roughly 3.5 h too early in each case, but the CHEOPS observations still manage to capture four complete transits (see Fig. B.1). Three of the visits catch the ingress of the transit, while one shows a partial egress. While we were unable to catch both features on a single visit, we are able to derive accurate transit mid-times to use in our subsequent TTV analysis. Notably, planet d is also predicted in visit 303, but not detected. This is likely due to the short duration of the visit, which prevented us from catching the transit.

Fig. 4 Modeling of CHEOPS visit 201. The detrended photometric points and its 10-min bins are represented by the light and dark blue symbols. The predicted transit times from A22 (±lσ) of each planet are represented by the vertical regions. The left panel shows the fit of planet e transit only. The model and its ±lσ uncertainties shown with the black curve and red region, respectively. The right panel shows the fit including transit of both planets c and e. The fitted transit of planet c and e are represented by the red and yellow curves, respectively, while the full model (±lσ) is represented by the black curve (and gray area). The two horizontal lines at the top of the panel represents the fitted transit mid-time for each planet and its ±lσ errors.

3.4 Planet e

This planet had the most accurate predicted transit times among all the planets in the system, with all its transits occurring within the estimated error range. Although visit 201 may have been affected by the potential transits of planets b and c, the other two visits are not anticipated to show any additional transits. However, due primarily to insufficient coverage during the ingress and egress phases of the transit, we faced limitations in the precision of the obtained transit mid-times (see Fig. B.2).

3.5 Planet f

Planet f appears most frequently in our dataset, being observed during visits 301, 302, 303, 304, 305, and 501. The first five visits were dedicated to observing this planet, while visit 501 targeted planet ɡ. Despite the number of visits dedicated to observing planet f, the ephemeris used to prepare these observations deviated significantly as illustrated in Fig. B.3. This results in predicted transit times later than observed. Consequently, we captured only one egress of the transit during visit 501. This discrepancy led to mid-transit times with deviations from a constant orbital period of up to ~25 min, with an average error of ~ 11 min, which allows for a TTV analysis of the system to be carried out.

3.6 Planet ɡ

When planning observations with CHEOPS, only two transits of the planet with a 42-day orbit had been observed with K2 (Christiansen et al. 2018; Lopez et al. 2019), and an additional transit was captured with Spitzer (HU21). However, the projected uncertainties of the transit times were significantly larger than 1.6 h based on these observations. To address this, the visit for planet ɡ was extended beyond the usual duration of CHEOPS visits, allowing us to also capture transits of planets d and f (see Sects. 3.3 and 3.5). After the transit modeling, we were able to determine the mid-transit time of planet ɡ with a precision of approximately ~ 12.25 min at epoch 32. This extends the observation coverage of the planet beyond its last observation at epoch 10 by Spitzer. Such comprehensive coverage is essential for estimating TTVs in the system, therefore, we have included this transit in our analysis.

4 Updated ephemerides

The new ephemerides of the planetary transits of this system make use of the K2, Spitzer, the Hale telescope, and the newly acquired CHEOPS data, including the transits from planets d, e, f, and ɡ. The K2 data and their respective mid-transit times were first used in Christiansen et al. (2018) and Lopez et al. (2019), but unfortunately, the individual transit mid-times were not reported by the authors. Thus, we reproduced the procedure used in both studies: the individual K2 transit mid-times were measured by fixing the model transit parameters to the best-fit values given in previous studies (Christiansen et al. 2018; Lopez et al. 2019; Hardegree-Ullman et al. 2021; Boyle et al. 2022) and allowing only the mid-transit time to vary. The uncertainties were calculated by computing the residuals from the best-fit model and performing a bootstrap analysis using the closest 100 timestamps, re-fitting the mid-transit time at each timestamp permutation. All the mid-times are listed in Tables C.1 and C.2 for each planet. For the rest of the data (Spitzer and Hale), the individual transit mid-times were provided directly by the respective authors, therefore, there was no need for a reanalysis.

With all the available transit times of each planet, we built the respective observed minus calculated (O–C) diagrams of the transits (Figs. C.1–C.5), where we compared the observed transit times with our updated ephemeris as described below. For each planet, the updated values of the average period and reference epoch (T₀) were obtained by a weighted linear fit of the transit times. The updated values of the average periods and T₀ are given in Table 5. The time residuals are also listed in Tables C.1 and C.2. The reported errors of the time deviations consider the uncertainty of the updated ephemeris and the respective error of the observed mid-time.

To improve the ephemeris of planet c, we utilized the transit times from K2 and three reference epochs (T₀) provided in Lopez et al. (2019), HU21 and A22 designated as epoch 0. All transit times are listed in Table C.1. If the mid-time of the potential CHEOPS transit in visit 201 is excluded, fitting the transit times with a linear function yields an updated ephemeris that diverges notably from those of HU21 and A22. This discrepancy is illustrated in Fig. C.1, which presents the O–C diagram of planet c’s transit times using our updated ephemeris. Our refined ephemeris aligns with the scenario of a non-detection during CHEOPS visit 201, as the planet was likely to have transited shortly before or at the very beginning of the light curve. Conversely, as previously discussed and as shown in Fig. C.1, assuming the detection of the transit in visit 201 results in a fitted time consistent with the predictions of HU21 ephemeris. Given the inability to definitively rule out any of these scenarios with our current data, any future transit monitoring of this planet ought to consider the large uncertainties in transit time predictions.

The ephemerides of planets d and e were refined using all available previous data in addition to observations from CHEOPS, as detailed in Tables C.1 and C.2 respectively. For comparison, we also show in their updated O–C diagrams (Figs. C.2 and C.3) the time projections from previous transit ephemerides. For planet d, our CHEOPS times do not agree with the linear ephemeris defined by considering only the K2 and B22 transits. When incorporating the CHEOPS data, the B22 transit times present a significant deviation of −158 ± 14 min and −48 ± 12 min for planets d and e, respectively, with respect to our updated linear ephemeris. This result is discussed in the next Sect. 5.

The six new CHEOPS visits for planet f effectively double its number of observations. Despite the incomplete coverage during these visits (as discussed in Sect. 3.5), we were able to improve the ephemeris of the planet. As illustrated in Fig. B.3, the ephemeris from A22 exhibited drift, resulting in predictions that were consistently too early by 121 ± 76 to 164 ± 77 min.

With only three previous observations, the ephemeris of planet ɡ were not very well constrained with projected uncertainties of hours, as illustrated in Fig. C.5. The new CHEOPS observation allows us to correct the ephemeris: it accounts for a deviation of around 45 min with respect to HU21 ephemeris for the CHEOPS’s epoch and also constrains the uncertainties of the projected transit times down to only few minutes (see Fig. C.5 and Table 5).

5 TTV analysis

Despite the fact that the number of transits observed are not enough to perform a full TTV analysis of the system and thereby enable a probe of the masses and/or eccentricities distribution of planets using photometry alone (e.g., Leleu et al. 2021; Hoyer et al. 2022), we conducted a first-order TTV analysis. To that end, we used the TTVfaster⁴ code (Agol & Deck 2016) for the TTV simulations combined with the nested sampling method (Buchner 2023). This code implements analytical formulas instead of numerical integrations to model the TTVs of multi-planetary systems, allowing for more efficient transit timing modeling of the planets in packed systems. This approach helps us constrain the expected TTVs and validate our updated orbital periods as reported in Table 5. Such modeling is particularly crucial, for example, in verifying the significant TTVs observed in planet d.

Our approach is well suited when exploring non-Gaussian parameters spaces. This is expected in our case, as we have a large number of fitted parameter compared to the number of transit times we have. Indeed TTVfaster requires seven input parameters per planet. For this we use the implementation of the UltraNest⁵ package (Buchner 2021), in particular, its stepsampler.SliceSampler method, which is more efficient for high-dimensional spaces. We analyzed the TTVs of the planets based on the timing information detailed in Sects. 3 and 4 and reported in Tables C.1 and C.2, combined with the masses values from A22. In addition, we used our updated planetary parameters as priors. With respect to the parameters we lacked updates for, such as those of planets b and c, we utilized the values reported in A22. Specifically, for planet b, we relied on the three time stamps provided by A22, as we have no significant detection of any transit of this planet; thus, we were not able to update its ephemeris. Therefore, little information from this planet can be used to constrain the planet-system parameters. The defined priors distribution for each parameter, as well as the posteriors, are detailed in Table 6. We report in this table the distribution of the relevant fitted planetary parameters: mean orbital period (P), mass (M), and eccentricity (e). In addition, the RMS and the larger amplitude of the best-fit TTV model are also reported.

The time residuals of the planetary system simulation are presented in the O–C diagrams of Fig. 5, along with the best-fit TTVs model. The posterior values of the fitted parameters are reported in Table 6. As also reported in the table, the models predict mean TTV amplitudes that are lower than 10 min for planets b, c, e, f, and ɡ, as depicted in Fig. 5. For planet d, however, the temporal deviations are as large as 60 min, with a modeled TTVs RMS of 25 min. The fitted planetary parameters of the modeled TTVs are fully consistent with the defined priors, particularly the mean orbital period, eccentricities, and masses. Planet d is the only notorious exception, where its fitted mass deviates to a small M_d ~ 0.8 M_⊕, and a relatively high eccentricity of e_d = 0.16 ± 0.04. These values were necessary to fit the large TTVs observed with CHEOPS.

It is clear that the transit times cannot be described purely with a strict constant period. Notably for planet d, as reported in Table C.1, B22’s transits deviate −158 ± 14 min and −48 ± 12 min from the updated ephemeris, while three out of the four CHEOPS transits deviate by more than 10 min (but with uncertainties that are on the same order). Assuming these deviations are a consequence of gravitational interactions between the planetary bodies, we can predict that the actual mean orbital period of planet d should therefore be between the B22 value and our updated estimate. This detection of large TTV indicates the significant dynamical interaction of planet d with the other planets in the system, as demonstrated by our simulations. New observations would be required to constrain the presence of these TTVs, which are expected for planets located in MMR (e.g., Agol et al. 2005; Xie 2013), as is the case of planet d. These observations would also be required to properly sample the super-period of planet d, and further constrain its eccentricity.

For planet f, the time residuals shown in Fig. C.4 seem to suggest a hint of TTVs. In fact, the CHEOPS transits show a more scattered distribution around the new linear ephemeris than the K2 points, though the time uncertainties prevent us to claim this is a significant TTV detection. Moreover, the TTV simulations only predict small periodic oscillations with amplitudes of 7.4 ± 3.2 min around a constant period. Thus, more precise measurements are required to confirm the presence of such deviations. No deviations were predicted from the simulations for planet ɡ as it is located further out in the system; therefore, it is subject to only weak gravitational influence from the innermost planets.

6 Conclusion

The observations obtained during this CHEOPS campaign represent the first large series of measurement on the system since the initial Kepler data reported three-and-a-half years prior. The newly computed ephemerides for the system represent a leap forward in accuracy. Particularly for planets d, f, and ɡ, the refined values are crucial for targeting future transit events. This is even more valuable as the accumulated deviations for these planets had already reached up to 200 min. For example, in the case of planet f, we examined the deviation between predictions by the end of 2024. Using the ephemerides from A22, we predict T_C = 2 460 673.05 ± 0.12. However, with the new average period derived from our work, we obtain T_C = 2460672.727 ± 0.015. The new ephemeris yield a T_C 8 h prior to that of A22.

We also want to emphasize the absence of transits of planet b in CHEOPS light curves. Indeed, despite the numerous transit events predicted within our observing windows, no transits were confirmed in the CHEOPS dataset. The non-detection of the transits of planet b in CHEOPS could result from a large accumulated error in the ephemeris or, for example, from a dynamical mechanism such as the precession of the planet. Therefore, a dedicated photometric follow-up is required to resolve this debate.

Planets in the K2-138 system typically display TTVs on the order of 10 min, except for planet d, which shows TTV amplitudes of up to 60 min compared to the mean orbital period. However, this is only possible if planet d is of a mass of ~0.8 M_⊕, with a fairly high eccentricity of ~0.16. As planet d is in a MMR with its neigbors, large TTVs can be expected. Therefore, a future analysis based on better sampled TTVs will provide a more detailed view of the architecture of this compact system and its stability.

Fig. 5 Results of TTV simulations of the K2-138 planetary system. The orange points represent the observed transit time residuals for each planet, reported in Tables C.1 and C.2. The best-fit TTV models obtained with TTVfaster and its ±1σ errors are represented respectively by the blue curve and region, and are detailed in Table 6.

Acknowledgements

H.V. and M.D. acknowledge funding from the Institut Universitaire de France (IUF) that made this work possible. SS acknowledges support from the Programme National de Planétologie (PNP) of CNRS-INSU. S.H. acknowledges CNES funding through the grant 837319. We thank the referee for their comments that improved the quality of the paper.

Appendix A System parameters

This appendix compiles the stellar and planetary parameters used in this study. Table 1 shows the stellar parameters of host star and Table A.1 lists the parameters of the planets in the system.

Appendix B Multi-transit figures

Here, we present all the light curves used during multi-transit analyses. All the figures show the best-fit models obtained with pycheops. For each figure the detrented data are shown as gray dots, with their 15 minute bins as black circles. The transit model is shown as a black-gray line. The transit are plotted in chronological order from bottom to top.

Fig. B.1 Transits and the best-fit model for the multi-visit of planet d.

Fig. B.2 Transits and the best-fit model for the multi-visit of planet e.

Appendix C O-C diagram

Fig. B.3 Transits and the best-fit model for the multi-visit of planet f.

Fig. C.1 Observed minus calculated (O-C) diagram of the transit times of planet c. The estimates from K2 (blue symbols) are compared to our updated ephemeris, which uncertainties are represented by the light blue region. The predicted times of the ephemeris from HU21 (green curve) and A22 (yellow curve and region) are also shown. The full duration of the CHEOPS visit 201 is marked with the vertical dashed curve and the corresponding fitted mid-time assuming a transit detection (see Sect. 3.2) is marked with the red symbol.

Fig. C.2 Observed minus calculated (O-C) diagram of the transit times of K2-328 d, compared to our updated ephemeris (light blue region), as described in Sect. 3.3. The times of the transits observed by K2, Hale observatory, and CHEOPS are shown with the blue, red, and orange symbols, respectively. The predicted ephemeris from A22 and B22 are represented with yellow and red curves and regions. The mean ephemeris obtained by the TTV simulations (Sect. 5) is represented by the gray dashed curve.

Fig. C.3 Observed minus calculated (O-C) diagram of the transit times of planet e. The K2 and CHEOPS transit times are shown with the blue and orange symbols and compared with our updated ephemeris (light blue region). The yellow curve and region represent the ephemeris of A22 and their uncertainties.

Fig. C.4 Observed minus calculated (O-C) diagram of the transit times of planet f. The K2 and CHEOPS transit times (blue and orange symbols) were used to refine the ephemeris reported by A22 (yellow curve and region). The uncertainties of the updated ephemeris are represented by the light blue region.

Fig. C.5 Observed minus calculated (O-C) diagram of the transit times of planet g. The K2, Spitzer and CHEOPS transit times (blue, red and orange symbols) allowed to refine the ephemeris reported by HU21, yellow curve and region). The uncertainties of our updated ephemeris are represented by the light blue region.

References

All Tables

All Figures

Fig. 1 Detгended light curve of CHEOPS visit 501. The colored regions correspond to the predicted transit time (with ±1σ errors) of planets d, f, 𝑔, and b according to the ephemerides from A22. The detrended CHEOPS data are represented by the light gray points, with their 5-min bins shown with the white dots. The black solid curve shows the fit of the four planet transits and arrows show the planets associated to each transit. Planet b is written in parentheses, as the depth of the event (~550 ppm) doesn't appear compatible with prior data (~250 ppm).

In the text

Fig. 2 Putative transits of planet b. From top to bottom: Light curves of CHEOPS visits 102, 201 and 304, in which planet b transits were expected to occur. Colored regions show the expected transits mid-time with their ±1σ errors. Transits due to other planets than b have already been removed. The solid black curve represent the best fits for planet b, while fixing its depth and duration to D = 252 ppm and T14 = 2.00 h.

In the text

	Fig. 3 Detrended light curve of CHEOPS visit 201. The solid black curve represents the fitted transit of planet e. The predicted mid-transit time of planet c is shown within ±1σ error as colored regions from two different sources: A22 in green and HU21 in yellow.
In the text

Fig. 4 Modeling of CHEOPS visit 201. The detrended photometric points and its 10-min bins are represented by the light and dark blue symbols. The predicted transit times from A22 (±lσ) of each planet are represented by the vertical regions. The left panel shows the fit of planet e transit only. The model and its ±lσ uncertainties shown with the black curve and red region, respectively. The right panel shows the fit including transit of both planets c and e. The fitted transit of planet c and e are represented by the red and yellow curves, respectively, while the full model (±lσ) is represented by the black curve (and gray area). The two horizontal lines at the top of the panel represents the fitted transit mid-time for each planet and its ±lσ errors.

In the text

	Fig. 5 Results of TTV simulations of the K2-138 planetary system. The orange points represent the observed transit time residuals for each planet, reported in Tables C.1 and C.2. The best-fit TTV models obtained with TTVfaster and its ±1σ errors are represented respectively by the blue curve and region, and are detailed in Table 6.
In the text

	Fig. B.1 Transits and the best-fit model for the multi-visit of planet d.
In the text

	Fig. B.2 Transits and the best-fit model for the multi-visit of planet e.
In the text

	Fig. B.3 Transits and the best-fit model for the multi-visit of planet f.
In the text

Fig. C.1 Observed minus calculated (O-C) diagram of the transit times of planet c. The estimates from K2 (blue symbols) are compared to our updated ephemeris, which uncertainties are represented by the light blue region. The predicted times of the ephemeris from HU21 (green curve) and A22 (yellow curve and region) are also shown. The full duration of the CHEOPS visit 201 is marked with the vertical dashed curve and the corresponding fitted mid-time assuming a transit detection (see Sect. 3.2) is marked with the red symbol.

In the text

Fig. C.2 Observed minus calculated (O-C) diagram of the transit times of K2-328 d, compared to our updated ephemeris (light blue region), as described in Sect. 3.3. The times of the transits observed by K2, Hale observatory, and CHEOPS are shown with the blue, red, and orange symbols, respectively. The predicted ephemeris from A22 and B22 are represented with yellow and red curves and regions. The mean ephemeris obtained by the TTV simulations (Sect. 5) is represented by the gray dashed curve.

In the text

	Fig. C.3 Observed minus calculated (O-C) diagram of the transit times of planet e. The K2 and CHEOPS transit times are shown with the blue and orange symbols and compared with our updated ephemeris (light blue region). The yellow curve and region represent the ephemeris of A22 and their uncertainties.
In the text

	Fig. C.4 Observed minus calculated (O-C) diagram of the transit times of planet f. The K2 and CHEOPS transit times (blue and orange symbols) were used to refine the ephemeris reported by A22 (yellow curve and region). The uncertainties of the updated ephemeris are represented by the light blue region.
In the text

	Fig. C.5 Observed minus calculated (O-C) diagram of the transit times of planet g. The K2, Spitzer and CHEOPS transit times (blue, red and orange symbols) allowed to refine the ephemeris reported by HU21, yellow curve and region). The uncertainties of our updated ephemeris are represented by the light blue region.
In the text