A new dynamical modeling of the WASP-47 system with CHEOPS observations

Among the hundreds of known hot Jupiters (HJs), only five have been found to have companions on short-period orbits. Within this rare class of multiple planetary systems, the architecture of WASP-47 is unique, hosting an HJ (planet -b) with both an inner and an outer sub-Neptunian mass companion (-e and -d, respectively) as well as an additional non-transiting, long-period giant (-c). The small period ratio between planets -b and -d boosts the transit time variation (TTV) signal, making it possible to reliably measure the masses of these planets in synergy with the radial velocity (RV) technique. In this paper, we present new space- and ground-based photometric data of WASP-47b and WASP-47-d, including 11 unpublished light curves from the ESA mission CHEOPS. We analyzed the light curves in a homogeneous way together with all the publicly available data to carry out a global $N$-body dynamical modeling of the TTV and RV signals. We retrieved, among other parameters, a mass and density for planet -d of $M_\mathrm{d}=15.5\pm 0.8$ $M_\oplus$ and $\rho_\mathrm{d}=1.69\pm 0.22$ g\,cm$^{-3}$, which is in good agreement with the literature and consistent with a Neptune-like composition. For the inner planet (-e), we found a mass and density of $M_\mathrm{e}=9.0\pm 0.5$ $M_\oplus$ and $\rho_\mathrm{e}=8.1\pm 0.5$ g\,cm$^{-3}$, suggesting an Earth-like composition close to other ultra-hot planets at similar irradiation levels. Though this result is in agreement with previous RV+TTV studies, it is not in agreement with the most recent RV analysis (at 2.8$\sigma$), which yielded a lower density compatible with a pure silicate composition. This discrepancy highlights the still unresolved issue of suspected systematic offsets between RV and TTV measurements. In this paper, we also significantly improve the orbital ephemerides of all transiting planets, which will be crucial for any future follow-up.

ions: WASP-47, the subject of this paper (see below), Kepler-730 (Zhu et al. 2018;Cañas et al. 2019), TOI-1130 (Huang et al. 2020), WASP-132 (Hord et al. 2022), and the most recently reported TOI-2000 (Sha et al. 2021).Together, these five planets represent less than 1% of the known HJs2 .Notably, WASP-47 is unique among the mentioned systems for showing clear transit time variations (TTVs), allowing researchers to jointly exploit TTVs and RVs to measure the planetary masses and other orbital parameters in a much more reliable way than using the two techniques separately (Malavolta et al. 2017).WASP-47 has been dubbed "the gift that keeps on giving" (Kane et al. 2020) due to the many layers of scientific investigation it has stimulated.The first member of this planetary system, HJ WASP-47b3 (1.14 ± 0.02 M jup , 1.127 ± 0.013 R jup , P 4.159 d) was discovered by Hellier et al. (2012) from ground-based photometric data (SuperWASP).A few years later, WASP-47 was observed by K2 during its Campaign 3, enabling the detection of two additional companions (Becker et al. 2015): WASP-47e (6.8 ± 0.7 M ⊕ , 1.81 ± 0.03 R ⊕ , P 0.789 d), an inner super-Earth, and WASP-47d (13.1 ± 1.5 M ⊕ , 3.58 ± 0.05 R ⊕ , P 9.031 d), an outer Neptunian.Around the same time, the results from the first extensive RV campaign (Neveu-VanMalle et al. 2016) revealed the presence of a fourth non-transiting Jupitersized planet on a much longer orbit, WASP-47c (1.25±0.03M jup , P 588 d).Subsequent RV analysis with fresh data have further refined and/or constrained the planetary masses (Dai et al. 2015;Sinukoff et al. 2017;Vanderburg et al. 2017).The opportunity of combining TTV and RV analysis in a joint approach was exploited by Almenara et al. (2016) and Weiss et al. (2017), the former within a photodynamical framework. 4All the published measurements of planetary masses for WASP-47 are summarized in Table 1.It is worth mentioning that the measurement of the spin-orbit angle obtained through the Rossiter-McLaughlin effect is consistent with an aligned configuration (λ = 0 • ± 24 • ) for WASP-47b (Sanchis-Ojeda et al. 2015) and that a phasecurve re-analysis of K2 photometry constrained the WASP-47b albedo to values significantly lower than the average of HJs (Kane et al. 2020).Lastly, a very recent analysis by Bryant & Bayliss (2022) focusing on the characterization of planet -e has made use of both TESS and ESPRESSO data for the first time.
Despite the extensive amount of work carried out on WASP-47 so far, a convincing and complete description of the evolutionary history of this system is still elusive.Our current understanding of this system will change once the James Webb Space Telescope (JWST) enables a detailed atmospheric characterization of WASP-47's transiting planets, especially of the small ones: -e and -d (Bryant & Bayliss 2022).Unfortunately, a reasonably accurate prediction of the future transits of planet -d is difficult, as even the best available ephemerides based on a combination of K2 and RV data will have drifted by about one hour (1-σ) at epoch 2023.0.The photometric detection of WASP-47d is indeed beyond the reach of ground-based facilities, and TESS has been unable to significantly detect the only transit predicted to fall in Sector 42 (and will not re-observe it at least until the end of Cycle 6 in 2024).The need to recover a new set of reliable ephemerides is one of the initial motivations of our investigation.In addition, the mass measurements of planet -d and -e have historically been slightly inconsistent with each other due to a small statistical tension between the most precise estimate through a photodynamical approach and studies wholly based on RV data (Table 1).Although the most recent ESPRESSO measurements (Bryant & Bayliss 2022) appear to have solved this tension for planet -d (but still not for -e), WASP-47 remains one of the very few systems for which a precise mass measurement can be achieved by either TTVs or RVs, potentially shedding some light on an old debate about a possible systematic discrepancy between the two techniques (Mills & Mazeh 2017;Petigura et al. 2018).
In this paper, we present new data from the CHaracterising ExOPlanet Satellite (CHEOPS) and new ground-based photometric data for WASP-47b and WASP-47d, and we perform a global dynamical re-analysis of all the available data, including the latest TESS and ESPRESSO data sets.We focus particularly on the determination of the orbital and physical parameters of planet -d and -e in order to address the aforementioned issues.In Section 2, we describe all the photometric and spectroscopic data analyzed in this work.In Section 3, we deal with the light curve fitting and the dynamical modeling.Finally, we present and discuss our results, including a comparison with the literature, in Section 4.

Observations
For our dynamical analysis, we gathered all the publicly available RV and photometric data of WASP-47.This included both public data (from K2, TESS, and several RV surveys) and new proprietary data from CHEOPS and ground-based telescopes, described below.A log summarizing all the analyzed photometric data is reported in the Appendix.The time stamps associated with all the measurements were consistently converted to the BJD TDB standard, as prescribed by Eastman et al. (2010).

CHEOPS photometry
Launched in 2019, CHEOPS (Benz et al. 2021) is an ESA Sclass mission currently carrying out its 3.5-year nominal observing program.The scientific instrument of CHEOPS is a 32-cm reflecting telescope designed to deliver defocused images for the performance of ultra high-precision photometric measurements of bright stars.The high-performance capabilities of CHEOPS when working with transit timing in particular have been demonstrated by Borsato et al. (2021).
Eleven visits of WASP-47 were scheduled with CHEOPS in 2020 and 2021 within the Guaranteed Time Observations (GTO) programs CH_PR100017 and CH_PR100025 (see Table A.4 for a list of the CHEOPS data sets).Ten visits were centered on the transits of WASP-47b (Table A.1), and one was centered on WASP-47d (Table A.3).All the light curves were extracted from the raw satellite data by the official CHEOPS Data Reduction Pipeline v13.1 (DRP; Hoyer et al. 2020), and their plots are shown in the left panel of Fig. 1 (planet -b) and in Fig. 3 (planetd).The exposure time was set to 60 s and the minimum efficiency to 50%, resulting in some of the light curves showing noticeable gaps every 98.77 min (the orbital period of CHEOPS) due to the spacecraft crossings of the South Atlantic Anomaly as well as Earth-related constraints (blockage of the line-of-sight by the Earth and high levels of stray light).

K2 and TESS photometry
During its Campaign 3, K2 monitored WASP-47 in shortcadence mode (net cadence: 58.3 s) for 69 nearly uninterrupted days, from November 14, 2014, to January 23, 2015.These data sets were used to study the transits of planets -d and -e for the first time (Becker et al. 2015).We also use their light curves in our analysis, as they are already corrected for jitter-induced systematic errors by fitting splines as a function of spatial drift, following the procedure by Vanderburg & Johnson (2014).Overall, 108 transit light curves from K2 were ingested in our analysis: 16 of planet -b, seven of planet -d, and 85 of planet -e.These are summarized in Table A .1, A.3, and A.5-A.6, respectively.Almost seven years later, WASP-47 was observed by TESS for the first time (and the only time so far) in Sector 42 (camera 1, CCD 4), from August 20 to September 16, 2021.Being a known planet host, it was designated as a 120-s cadence target.The resulting data set is the same investigated by Bryant & Bayliss (2022), who used of the Pre-search Data Conditioned Simple Aperture Photometry (PDCSAP) light curve made available by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016).Unfortunately, the PDCSAP light curve is missing a large one-week section, including two transits of WASP-47b (on August 27 and 31), due to a "scattered light" quality flag.Unlike Bryant & Bayliss (2022), we based our analysis on Simple Aperture Photometry (SAP) data, where the two missing transits are preserved and do not appear to be impacted by a significant amount of systematic noise, as is demonstrated in Section 3.1.While the SAP algorithm does not correct the light curves for photometric dilution due to contaminating sources, we emphasize that no Gaia DR3 source (limiting magnitude G 21) is present at all within a radius of 30 from WASP-47, and the brightest one within 60 has G = 17.50, implying a negligible dilution factor.We also note that a constant dilution factor does not have any systematic effect on the measurement of transit times because it does not alter the symmetry of the transit signal.
Overall, six TESS light curves of planet -b (Table A.1; right panel of Fig. 1) were ingested.As already noted by Bryant & Bayliss (2022), individual transits from planet -d and -e are hidden in the TESS photometric scatter, and their detection is at most marginal.Our preliminary fitting tests confirmed this, and we did not attempt to retrieve their timings.The next visit to WASP-47 by TESS is not expected until at least the end of Cycle 6 (October 2024).

Ground-based photometry
One transit of WASP-47b was observed on September 24, 2021, by the 60-cm Rapid-Eye Mount telescope (REM; Zerbi et al. 2001) through the ROS2 instrument (Molinari et al. 2014), a charge-coupled device (CCD) camera able to observe the same field of view (9.1 × 9.1 ) simultaneously in four photometric bands (Sloan g , r , i , and z ).These observations (PI: Nascimbeni) were carried out as part of the TASTE program, a groundbased multi-site long-term TTV campaign to monitor transiting planets (Nascimbeni et al. 2011).Two bright reference stars (TYC 5805-338-1 and TYC 5805-739-1) were always imaged along with the target, allowing us to perform differential photometry.In addition, in order to mitigate systematic errors from guiding drifts, pixel-to-pixel non-homogeneity, and flat-field errors, a generous amount of defocus was applied (FWHM 10 pix).
The four transit light curves from the observation were extracted by running STARSKY (Nascimbeni et al. 2013), the photometric pipeline developed for TASTE.The STARSKY pipeline automatically selects the best combination of photometric aperture radii and weighting scheme to minimize the final out-of-transit scatter.No linear detrending against the external parameters automatically tested by STARSKY proved to be effective (i.e., linear baseline, airmass, x-y position, background level, stellar FWHM).Guiding drifts unfortunately turned out to be much larger than anticipated, totaling more than 30 pixels throughout the whole series.The resulting systematic errors impacted the g , i , and z data to an unacceptable extent, while the r light curve was mostly spared, thanks to its much better flat-field correction and its structure being more homogeneous at small spatial scales.For this reason, we only considered the r light curve (right panel of Fig. 2) for our subsequent analysis.
Seven additional transit light curves of WASP-47b (of six distinct transit events) were collected from 2020 to 2021 by Y. Jongen and A. Wünsche and are available on the public archive of the Exoplanet Transit Database (ETD; Poddaný et al. 2010).They were gathered through an R filter with a PlaneWave CDK 17" telescope coupled with a Moravian G4 16k CCD camera at the Deep Sky Chile facility (Jongen) and with an unfiltered PlaneWave CDK 20" telescope coupled with the same camera at the El Sauce Observatory in Chile (Wünsche).The exposure time was set to 120 s by both observers.After visually checking their quality and converting the time stamps to BJD TDB , the transit light curves were included in our analysis as well.These light  curves are detailed in Table A.2 and plotted in the left panel of Fig. 2.

Light curve fitting
All the photometric data described in the previous section were homogeneously (re-)analyzed through the PyORBIT (Malavolta et al. 2016(Malavolta et al. , 2018) ) code.The CHEOPS light curves also went through an additional detrending stage using cheope5 , an optimized python tool (of which pycheops is the back-end; Maxted et al. 2022) specifically developed to filter, correct, and fit CHEOPS data.The detrending model of cheope is defined as a linear combination of terms, including a quadratic baseline f 0 + d f /dt + d 2 f /dt 2 , the first and second order derivative of the centroid offset in x and y pixel coordinates ), and the first three harmonics of the spacecraft roll angle (in cos φ and sin φ).There is also an additional term called "glint" that empirically models the internal reflections as a spline-based smooth function (Borsato et al. 2021).The combination of terms to be adopted as the best detrending model was selected automatically by cheope according to its Bayes factor, first through an lmfit optimization (Newville et al. 2014) to find a reasonable starting point and then with an emcee fit (Foreman-Mackey et al. 2013).
For both CHEOPS and all the other data sets, the stellar activity signal and any residual instrumental trend were subtracted using wōtan (Hippke et al. 2019) after masking the transits of planets -b, -d, and -e to get the normalized light curves.The filter adopted was a Tukey's bi-weight with a window filter duration set equal to the transit duration plus its 1-σ uncertainty.All the detrended transits from ground-and space-based telescopes are shown in Figs. 1, 2, and 3.As for the TESS Sector 42 light curves, we focused our analysis on the extraction of six transits of planet -b only, after masking the transits of planets -d and -e, which were only marginally detectable in this data set and hence not usable for an individual T 0 fit.
In order to extract the central transit times T 0 , we applied PyORBIT independently on each of the five data sets (i.e., CHEOPS, TESS, K2, REM, and ETD), performing a modeling of each light curve using PyDE+emcee (Parviainen et al. 2016;Foreman-Mackey et al. 2013).We assumed Keplerian orbits for all three fitted planets, -b, -e, and -d.We fixed the a/R parameter of each planet at the values found on the highest S/N data set (K2; Vanderburg et al. 2017), leaving the individual T 0 values, the planetary radius R p /R , and the orbital inclination i for each planet as the only free parameters in order to take into account changes in the impact parameter as a function of time or different transit depths through different instruments and filters.The stellar parameters were derived by the CHEOPS Stellar Characterization team following the procedure described by Borsato et al. (2022), Section 3.2.1.The two parameters that are crucial for our dynamical analysis, the stellar radius (R ) and mass (M ), are listed in Table 2.We adopted a quadratic law u 1 , u 2 for the limb darkening coefficients, but we re-parameterized them as fitting parameters q 1 and q 2 , as done by Kipping (2013).The priors set on the coefficients u 1 and u 2 are reported in Table 2.We assumed a Gaussian prior distribution over the coefficients for all the instruments using a bi-linear interpolation of the limbdarkening profile defined in Claret (2017Claret ( , 2021)).We used PyDE to find the best starting point for the MCMC sampler.For each data set, PyDE used a population size of 4 × N par , where N par is the number of free parameters, and let it evolve for 64 000 generations, and for the Bayesian sampler, we set 100 chains for 150 000 steps, discarding the first 40 000 by adopting a thinning factor of 100.All the final T 0 values of this analysis are summarized in Tables A.1 to A.6 of Appendix 4, while the physical best-fit values of the common parameters among individual data sets (R p , i) are shown in Table 3.As anticipated in Section 2, we directly compared the transit timings we extracted from TESS SAP photometry with those obtained from the PDCSAP light curves by Bryant & Bayliss (2022) and listed in their Table 4 to check for any systematic offset.The average ∆T 0 = T 0,PDCSAP − T 0,SAP evaluated over the four transits in common yielded 14 s (with average error bars of ∼ 66 s), and no pair disagreed by more than 0.5 σ.In other words, timings from SAP and PDCSAP are in perfect statistical agreement.

Dynamical modeling
We employed the TRADES dynamical integrator (Borsato et al. 2014), which has already been successfully applied to Kepler/K2 (Borsato et al. 2019) and CHEOPS data (Borsato et al. 2021) as well as to simulate the Ariel TTV science case (Borsato et al. 2022).The TRADES dynamical integrator allowed us to fit transit times (T 0 ) and RVs simultaneously during the integration of the orbits.We fitted the mass ratios (M p /M ), the periods (P), and the mean longitude (λ) for all the planets.The eccentricity (e) and the argument of periastron (ω) were fitted as ( √ e cos ω, √ e sin ω) for planets -b, -c, and -d, while we used the (initial) circular orbit (e = 0, ω = 90 • ) for planet -e due to its extremely short tidal circularization timescale (Vanderburg et al. 2017).We adopted our initial values for M b , M c , M d , M e , e and ω from Almenara et al. (2016).We fixed the longitude of the ascending node (Ω) of -b, -c, and -e to 180 • .For the orbital inclination i, we took the weighted average of the fitted values from Section 3.1 and Table 3 (for instance, i b = 88.88 ± 0.14 • ).We found that by assuming i d < 90 • , one of the most recent transits of planet -b was missed by our integration (i.e., the impact parameter for that transit was greater than one).Almenara et al. (2016) reported a best-fit inclination higher than 90 • , but the authors specified that the supplementary angle is equally probable.If both planets -b and -d transit in the same stellar hemisphere (with respect to the observer's line of sight), then their relative distance would be shorter with respect to the case with planet -d on the opposite hemisphere.The mutual gravitational interaction would of course be much stronger in the former case, moving planet -b on some occasions to a non-transiting configuration that is not observed in our data.Therefore, we assumed that planet -d transits the opposite hemisphere with respect to -b (i b < 90 • ) and adopted the value 180 • − i d as the initial parameter for TRADES.We also fitted the i and Ω of -d, and we set i c to 90 • .All the orbital parameters are astrocentric and defined at the reference time 2456979.5 BJD TDB , the same as in Almenara et al. (2016).The initial periods were set to the linear ephemeris fitted to T 0 s determined in Section 3.1.For each RV data set, we fitted a log 2 -based jitter and an RV offset (γ).However, we found that changing the starting point of the jitter (e.g., log 2 of 0.5 and < 1 × 10 −6 ms −1 ) did not affect the final distribution of the chains and, therefore, did not change the best-fit solution.We split the CORALIE data set into three subsets, as done in Almenara et al. (2016), for a total of seven RV data sets.In all, we fitted 34 parameters, all with uniform-uninformative priors in the parameter space.
In order to save computational time, we first ran a local minimizer 6 and used the result as a starting point for the python package emcee (Foreman- Mackey et al. 2013Mackey et al. , 2019)), for which we initialized 68 walkers with a tight Gaussian.We adopted as a sampler algorithm a mix of the differential evolution proposal (80% of the walkers; Nelson et al. 2014) and the snooker differential evolution proposal (20% of the walkers; ter Braak & Vrugt 2008).We let the code run for 110 000 steps and discarded the first 100 000 steps as burn-in after checking the convergence of the chains by means of visual inspection as well as through diagnostics from Geweke 1991 (within 2-3σ) and Gelman-Rubin statistics7 (Gelman & Rubin 1992).Our uncertainties are computed as the high density interval (HDI) at 68.27% from the posterior distribution as the equivalent of the credible intervals at the 16th and 84th percentile.We defined our best-fit solution as the maximum likelihood estimator (MLE), that is, the parameter sample that maximizes the log-likelihood (log L) of the posterior distribution within the HDI.We also computed a symmetric uncertainty (σ) as the 68.27th percentile of the (sorted) absolute residuals of the posterior with respect to the best-fit solution.(See Table 4 for a summary of the fitted and physical parameters of the system determined with TRADES.See the O-C diagrams of planet -b, -d, and -e in Figs. 4, 5, and 6 and the RV plot in Fig 7 .)

Discussion and conclusions
The global analysis presented in this work is the most comprehensive dynamical modeling of the WASP-47 planetary system carried out so far.We exploited all the available RV and TTV data sets and merged them with 19 unpublished light curves, Notes.The columns give the name of the parameter, its best-fit value with its associated 68.27% HDI and σ, and the priors adopted.Astrocentric orbital parameters are defined at 2456979.5 BJD TDB .M is the mean anomaly computed as M = λ − ω − Ω.We fixed the value of i c = 90 • , so M p of planet -c is the minimum mass.Planet -e has a circular orbit (e = 0, ω = 90 • ) at the reference time, so e amd ω have not been fitted with TRADES.See Section 3.2 for details.
Article number, page 7 of 17  4) can be compared with the values reported in the literature and listed in Table 1.
For the two low-mass planets (-d, -e), we found best-fit masses (M d /M ⊕ = 15.5 ± 0.8, M e /M ⊕ = 9.0 ± 0.5) with rela-tive errors of about 5%, which is smaller than all the previous relative errors and is now limited mostly by our knowledge of the stellar mass (relative error on M : ∼4%).In particular, for the innermost planet, our solution is perfectly consistent with the upper extreme M e /M ⊕ = 9.1 ± 1.0 set by Weiss et al. 2017 (joint RV plus TTV analysis), but it is not consistent with the lower extreme M e /M ⊕ = 6.8 ± 0.6 set by Bryant & Bayliss 2022 using RVs only, apparently confirming, at 2.8σ significance, some kind of systematic offset between the two techniques (Mills & Mazeh 2017;Petigura et al. 2018).We did not see a similar behavior for planet -d, whose best-fit mass is fully consistent with all the previous measurements but more precise.By calculating the bulk densities as derived quantities from Table 4 and propagating the errors, we get ρ d = 1.69 ± 0.22 g cm −3 and ρ e = 8.1 ± 0.5 g cm −3 , thus confirming a Neptune-like density and composition for planet -d (implying an extended volatile envelope) but also moving planet -e very close to what planetary structure models predict to be an Earth-like composition rather than a pure-silicate one (in contrast to the claims by Vanderburg et al. 2017 andBryant &Bayliss 2022; see Fig. 8).If that conclusion were true, then planet -e would not appear as an outlier in the planetary density versus stellar irradiation diagram anymore and there would be no reason to hypothesize a different formation path with respect to other ultra-short period planets at similar irradiation levels, such as K2-141b and HD 213885b (Malavolta et al. 2018;Espinoza et al. 2020).Moreover, WASP-47 is a rather metal-rich star at [Fe/H] = +0.36 ± 0.05, according to Mortier et al. (2013).The high density of WASP-47e, which fits into the "super-Mercury" class, would confirm the finding by Adibekyan et al. (2021) that the bulk density of super-Earths correlates with stellar iron fraction.On the other hand, the relatively low density of planet -d would also confirm the trend discovered by Wilson et al. (2022), who found that sub-Neptunes around metal-rich stars have, on average, lower densities.Taken together, our findings would support the general idea that a higher stellar metallicity leads to forming planets with larger cores and, hence, generates denser hot super-Earths (such as -e) and warm sub-Neptunes that are able to retain more extended atmospheric envelopes (-d).
We caution, however, that the existing tension on M e could be due to subtle methodological biases in RV and/or RV plus TTV analyses that are yet to be fully explored and a deeper un-derstanding of the effects at play in such dynamically complex systems is needed.This is true of course for not only WASP-47 but for all planetary systems for which both techniques can be applied and that are steadily growing in number, especially after the launch of K2 and TESS.Two factors are most frequently discussed when dealing with such issues: the impact of stellar activity and the way RV and TTV information is merged at the analysis stage.Notably, WASP-47 is not a particularly active star, and all the studies up to and including Vanderburg et al. (2017) explicitly neglected any contribution from stellar activity on both photometric and spectroscopic data.Bryant & Bayliss (2022), on the other hand, noticed an excess scatter in the ESPRESSO data and concluded that the inclusion of a Gaussian process kernel in their RV modeling was justified by a lower BIC value of their fit, even though the statistical significance of the peak on the periodogram of the RV residuals at the claimed rotational period is marginal.The treatment of stellar activity, however, does not appear to be the (only) explanation for the M e discrepancy we see since Bryant & Bayliss (2022) came to a value that is almost identical to that reported by Vanderburg et al. (2017) without modeling for it: M e = 6.8 ± 0.7 vs. 6.8 ± 0.6 M ⊕ .Beyond that, in an RV plus TTV study such ours, it is extremely difficult to model stellar activity in a consistent framework because the two techniques are effected in different ways and at different timescales (Boisse et al. 2011;Oshagh et al. 2013;Ioannidis et al. 2016); indeed, transit light curves are impacted more by the local rather than global distribution of star spots over the photosphere.With WASP-47, an additional issue prevents us from adopting a more complicated data model (i.e., including stellar activity and/or a photodynamical approach): computational time.Having to deal with four planets, 34 free parameters, and an extensive set of 133 transits and 212 RVs, each of the MCMC steps of the optimization process described in Section 3.2 took about 15 s on a medium-power computing workstation, and the whole fit required weeks to reach convergence.This sizable computing time is mostly due to an unfortunate combination of a very long observing baseline (15 years) and the very short period of the inner planet (P e 0.79 d), forcing the N-body integration time step to unusually small values.Of course, under these assumptions, the inclusion of a Gaussian process treatment of stellar activity or, even worse, the implementation of a photodynamical algorithm would increase the computational time up to an unreasonable amount, at least with ordinarily available hardware.We emphasize, however, that our result agrees with previous RV plus TTV studies both based (Almenara et al. 2016) and not based (Weiss et al. 2017) on a photodynamical approach, so fitting our light curves and performing the dynamical modeling at separate stages cannot be the sole reason for the discrepancy discussed above.
As for the giant planets (-b and -c), it is worth noting that our newly derived masses are compatible with those presented by Vanderburg et al. (2017) and Bryant & Bayliss (2022) and statistically consistent with all the previous measurements reported in Table 1.The precision of our values M b /M ⊕ = 374 ± 17, M c /M ⊕ = 447 ± 20 is larger by a factor of about two with respect to the most recent RV works.This is easily understood if we consider that our error on M (4.4%) is larger than theirs (3%) and our dynamical analysis probes are M p /M rather than M p /M 2/3 , the latter being true for a pure RV analysis.In other words, our study is more sensitive to the uncertainty on stellar mass.The same is true for the derived planetary density ρ b = 0.98 ± 0.09 g cm −3 .An important byproduct of our analysis is the derivation of new, updated mean ephemerides8 for all the WASP-47 transiting planets, as this data can be exploited by any future follow-up study.Our best-fit relations are: Observations of WASP-47 are not yet scheduled with JWST, at least during Cycle 1, even though its transiting planets -e, -b, and -d have been recognized as suitable targets for transmission spectroscopy (Bryant & Bayliss 2022) and will very likely be proposed in the next cycles.Such observations will be crucial to test our planetary formation theories.It is widely accepted that three main channels for HJ migration exist (Dawson & Johnson 2018): one is through the protoplanetary disk at low eccentricity, another is through high-eccentricity dynamical migration with tidal damping at later times, and the third is through in situ formation close to or at the final orbit.A planet migrating through the disc or having formed in situ would accrete a large amount of gas from inside the ice line, whereas an HJ that migrated dynamically may have originated from beyond the ice line with a different composition (Öberg et al. 2011;Madhusudhan et al. 2017;Knierim et al. 2022).Simulations show that the dynamical migration channel would invariably destroy any close-in planets, as the orbit of the HJ would intersect with that of the close-in planets prior to circularization (Mustill et al. 2015).Thus, unlike ordinary HJs, for very rare systems such as WASP-47, we can rule out one of the main scenarios and test the other models through detailed atmospheric characterization.
Other space-based follow-up opportunities include, at least in principle, TESS, PLATO (Rauer et al. 2014), and ARIEL (Tinetti et al. 2018).As for TESS, we mentioned in the introduction that no other sector will include WASP-47, at least up to the end of Cycle 6, implying that this target will not be reobserved until fall 2024, at best.Unfortunately, individual transits of planets -e and -d are essentially undetectable by TESS due to its noise properties, and the time baseline of each sector is limited to about four weeks.Future sectors, if any, will thus not be effective in significantly extending the TTV analysis carried out so far.
On the other hand, PLATO will be able to deliver highprecision transit light curves of planets -e and -d since the predicted noise level ranges from 90 to 180 ppm in one hour, according to the number of co-pointing PLATO cameras involved (Montalto et al. 2021).Due to its very low ecliptic latitude, WASP-47 will not lie within the proposed Long-duration Observing Phase fields of PLATO (Nascimbeni et al. 2022), where targets will be continuously observed for at least two years.However, WASP-47 could be selected for the Short-duration Ob- Fig. 8. Mass-radius diagram with all the known planets with a relative error on bulk density better than 20% (circle points) color coded according to their equilibrium temperature.The colored lines represent models at different bulk compositions from Zeng et al. (2019).The two measurements of WASP-47e by us and by Bryant & Bayliss (2022) are plotted as squares.
serving Phase and monitored for two to three months without any significant interruption (see Fig. 1 of Heller et al. 2022).
As a closing note, ARIEL, though focused on the atmospheric characterization of exoplanets through emission and transmission spectroscopy (Tinetti et al. 2021), has also been proposed as a very powerful instrument to investigate TTVs (Borsato et al. 2022).WASP-47 is consistently included in the provisional ARIEL target list (Zingales et al. 2018;Edwards & Tinetti 2022) as a "Tier 2" object (aim: in-depth atmospheric characterization through repeated observations), but these observations are only devoted to observing eclipses of planet -b.

Fig. 1 .-
Fig.1.Space-based light curves of WASP-47b from CHEOPS and TESS analyzed for the present work after de-trending.For each light curve, the corresponding label shows the acquisition date in UT and the timing error σ T 0 (in seconds).Arbitrary vertical offsets of 0.0165 and 0.03 were added, respectively, to both sets for visualization purposes.

Fig. 4 .Fig. 5 .
Fig. 4. Transit times of WASP-47b.Upper panel: Observed -Calculated (O − C) diagram.The O − C was calculated by subtracting the predicted T 0 calculated from the linear ephemeris to the observed transit times.The O − C values computed from the observed T 0 s are plotted as open orange circles, while the O −C of the TRADES-simulated T 0 s from the best-fit model are plotted as blue circles.Samples drawn from the posterior distribution within HDI are shown as grey lines.Lower panel: Residuals computed as the difference between observed and simulated T 0 s.

Fig. 6 .
Fig. 6.Transit times of WASP-47e.Elements are the same as in Figs. 4 and 5, but for planet -e.

Fig. 7 .
Fig. 7. Radial velocities of WASP-47.Upper panel: RV plot.There is a different marker and color for each data set.The TRADES best-fit RV model is plotted with open black shapes.The maximum-likelihood orbital solution is shown by a light gray line.The CORALIE data set was divided into three data sets, as described in Almenara et al. (2016).Lower panel: RV residuals with respect to the TRADES RV best-fit model.The corresponding jitter determined from the best-fit model has been added in quadrature to the measured uncertainty of each data point.

T
0,b = 2459407.761987± 0.0000002 BJD TDB (1) +N × (4.159151 ± 0.0000004), T 0,d = 2457024.430107± 0.0000003 BJD TDB(2) +N × (9.030501 ± 0.000013), T 0,e = 2457012.13666± 0.00003 BJD TDB(3) +N × (0.78961 ± 0.00003), where N is an integer number commonly known as "epoch" and set arbitrarily to zero at our reference transit time T ref .We highlight that the 1σ uncertainty on transit prediction at 2023.0 for planet -d has now been reduced from about one hour to just six minutes.Of course this is crucial for planning space-based observations, where the cost of the invested observing time and the time-sensitive nature of the observations makes pinpointing the transit window as precisely as possible critical.

Table 1 .
Mass measurements of the WASP-47 planets published in the literature.

Table 2 .
Stellar and limb darkening parameters adopted in the transit fitting process (Section 3.1) and subsequent dynamical analysis (Section 3.2).

Table 3 .
Best-fit values of planetary radius R p (derived from R p /R ) and orbital inclination i.All these parameters were fitted as common parameters among individual data sets.
Notes.See Section 3.1 for details about the transit fit.

Table 4 .
Summary of the best-fit MLE parameters for the WASP-47 system determined with TRADES.