Discovery and characterization of the exoplanets WASP-148b and c. A transiting system with two interacting giant planets

We present the discovery and characterization of WASP-148, a new extrasolar system including at least two giant planets. The host star is a slowly rotating, inactive late-G dwarf with a V=12 magnitude. The planet WASP-148b is a hot Jupiter of 0.72 RJup and 0.29 MJup transiting its host with an orbital period of 8.80 days. We found the planetary candidate with the SuperWASP photometric survey then characterized it with the SOPHIE spectrograph. Our radial velocity measurements subsequently revealed the presence of a second planet in the system, WASP-148c, with an orbital period of 34.5 days and a minimum mass of 0.40 MJup. No transits of that outer planet were detected. The orbits of both planets are eccentric and fall near the 4:1 mean-motion resonances. That configuration is stable on long time scales but induces dynamical interactions so the orbits slightly differ from purely Keplerian orbits. In particular, WASP-148b shows transit timing variations of typically 15 minutes making it one of the few cases with TTVs detected on ground-based light curves. We establish the orbital plane of both planets could not have a mutual inclination larger than 35 degrees, and the true mass of WASP-148c is below 0.60 MJup. We present photometric and spectroscopic observations of that system covering a time span of ten years, and their Keplerian and Newtonian analyses; they should be significantly improved thanks to future TESS observations.


Introduction
Extrasolar planets which transit their host stars are especially interesting. When they are characterized both in photometry and spectroscopy they allow numerous studies and the determination of many of their parameters, including their radius and mass. About 700 exoplanets have such a double characterization according the exoplanet archives 1 . A few were first discovered from radial-velocity (RV) surveys then photometric follow up subsequently revealed their transits (e.g. Charbonneau et al. 2000, Motalebi et al. 2015. But the vast majority were first identified from photometric surveys then characterized with RV follow up. Spectroscopic observations of planetary candidates revealed by photometry are used to establish or reject their planetary nature, in particular by measuring the mass of the transiting bodies using the RV method (e.g. Hébrard et al. 2014, Cooke et al. 2020. They are also used to measure the orbits' eccentricity and obliquity and characterize the host stars. Long-term RV follow up could also reveal the presence of additional, nontransiting planets in the system (e.g. Christiansen 2017, Rey et al. 2018. Such multi-planetary systems are particularly interesting for the studies of their dynamics. Dynamics can also be studied in multi-planetary systems when transit timing variations (TTVs) are detected. Indeed, whereas in a purely Keplerian orbit the epoch of a transit is exactly periodic, several gravitational perturbations can produce small deviations of the transit epochs with respect to a perfect periodicity. These Newtonian orbits could be caused by orbital decay due to tides (e.g. Birkby et al. 2014) or to gravitational interactions between bodies in multiple systems (e.g. Holman & Murray 2005, Agol et al. 2005. In the case of multi-planetary systems, TTVs have larger amplitudes when the planets have orbital periods in near-commensurability, i.e. in mean-motion resonances (MMR). This makes TTV analysis a powerful technique for characterizing such systems, and in particular measuring planetary masses and eccentricities, or even detecting additional, perturbing planets.
For years, several attempts have been made to detect TTVs of transiting planets with ground-based photometry (e.g. Díaz et al. 2008, Maciejewski et al. 2010) but most of them were unconfirmed hereafter (e.g. Petrucci et al. 2020). TTVs indeed are difficult to identify as transit timing strongly depends on the steepest portions of the light curves (the planetary ingress and egress) which are short-duration events easily subject to systematics, especially from the ground. The high-quality, long-duration light curves of the Kepler space telescope finally allowed the first detection of TTVs by Holman et al. (2010) with the famous case of Kepler-9. That star is transited by two giant planets each ∼ 19.2 and ∼ 38.9 days. That near 2:1 MMR causes TTVs with an amplitude of about one day clearly detected with Kepler.
Now, a few tens of exoplanets have been detected or characterized thanks to their TTVs. Most of those detections were made from Kepler photometry. This includes multi-planetary transiting systems as well as single transiting planets showing TTVs allowing the detection and characterization of additional, non-transiting planets. Notable cases include KOI-142b, a P = 10.95-d planet showing TTVs up to one day caused by a non-transiting, giant planet with a period two times longer (Nesvorny et al. 2013) eventually detected in RVs (Barros et al. 2014), or the seven Earth-size planets orbiting the star Trappist-1 with periods between 1.5 and 18.7 days and TTVs up to a few tens minutes (Gillon et al. 2017).
Most of TTV systems imply Earth-or Neptune-size planets. Hot Jupiters presenting TTVs remain rare, which is one of the reasons why TTVs were late to be detected as the first researches were mainly attempted on such kind of systems. Only 1 exoplanet.eu, exoplanetarchive.ipac.caltech.edu. three are confirmed today. WASP-4b and WASP-12b are hot Jupiters of 1.3-d and 1.1-d periods, respectively; both show long-term deviations from a purely periodic orbit (Bouma et al. 2019, Southworth et al. 2019, Maciejewski et al. 2016, Yee et al. 2020. Several scenarios have been proposed to explain those TTVs, whose amplitudes are below two minutes. They include stellar activity, tide-caused orbital decay, or additional companions. In the case of WASP-4, Bouma et al. (2020) recently showed that it could be mostly or entirely produced by the line-of-sight acceleration of the system (see, however, Baluev et al. 2020). The third case is WASP-47b, a hot Jupiter on a 1.1-d period showing TTVs of half a minute explained by the presence of two smaller, short-period planets (Becker et al. 2015, Weiss et al. 2017.
Here we present the new planetary system WASP-148, a fourth case with a hot Jupiter showing TTVs. WASP-148b was first identified as a promising transiting planet candidate by the SuperWASP photometric survey with an orbital period of 8.80 days. The RV follow up with the SOPHIE spectrograph established the planetary nature of the transiting object, and revealed the presence of a second, outer, giant planet. WASP-148c orbits the host star with a period of 34.5 days which is near the 4:1 MMR, apparently without transiting it. The few photometric transits of WASP-148b observed with ground-based telescopes reveal significant deviations from a constant orbital period, with TTV amplitudes of the orders of a few minutes. They are likely to be due to WASP-148c, at least partially.
We present the observations of WASP-148 in Sect. 2, determine the properties of the host star in Sect. 3, and assess the evidence for the presence of the two planets in Sect. 4. Section 5 describes the Keplerian fit to the data whereas Sect. 6 discusses dynamic analyses of the system before we conclude in Sect. 7.  . It observes with a broad-band filter (400-700 nm), and secured thousands of photometric points over several seasons per star. Periodic signatures of possible planetary transits are identified in these light curves using the algorithms presented by Collier . Using this facility and procedure, WASP-148 was identified as the host star of a promising candidate for a transiting planet .  Three similar planetary-transit-like features were observed on  2008 June 20, 2010 June 04, and 2011 May 31, with a depth of ∼ 0.0070 mag, a duration of ∼ 3 hr, and a possible periodicity of 8.80 days. The catalog IDs, coordinates, magnitudes, and distance of the star WASP-148 are reported in Table 1. The SuperWASP light curves are shown in the three first plots of Fig. 1; the transit-like features are not obvious in those initial data, which illustrates how sensitive the candidate detection algorithms should be.
We note that the SuperWASP total light curve also shows a hint of a long-term sinusoidal modulation with a period near 24 days, which could be due to the stellar rotation. Its lowamplitude is not visible on the light curves plotted in Fig. 1 which just focus around short-duration transits.

Radial-velocity follow-up with SOPHIE
After its identification from SuperWASP photometry, we started a RV follow up of WASP-148 with the SOPHIE spectrograph at the 1.93-m telescope of Observatoire Haute-Provence, France. The goal was first to establish the putative planetary nature of the transiting candidate, and then in case of positive detection to characterize the planet by measuring notably its mass and orbital eccentricity, as well as to search for potential additional bodies in the system. SOPHIE is a stabilizedéchelle spectrograph dedicated to high-precision RV measurements (Perruchot et al. 2008, Bouchy et al. 2009a. Here we used its High-Efficiency mode with a resolving power R = 40 000 and slow readout mode to increase the throughput for this faint star. The first season of observation (35 SOPHIE measurements over six months) revealed RV variations in phase with the 8.80day signal detected with SuperWASP, no significant variations in the spectral line profiles, and an amplitude of the RV variations of the order of 30 m/s. This showed that the transiting body is a planet slightly more massive than Saturn, designated as WASP-148b hereafter. However, the residuals of the one-planet Keplerian fit exhibited a dispersion significantly larger than expected, with a possible periodicity.
The second season allowed us to secure a total of 75 measurements over 18 months. The dataset clearly showed that a second periodic signal was present in the RVs, with a period around 35 days, here again with no significant variations in the spectral line profiles and an amplitude of the RV variations of the order of 30 m/s. This indicated the star hosts a second, outer planet slightly more massive than Saturn, hereafter designated as WASP-148c The final dataset we present here includes 116 SOPHIE measurements secured between April 2014 and June 2018. Depending on weather conditions, exposure times ranged from 200 to 2000 seconds with a typical value of 1200 seconds. This allowed a nearly constant signal-to-noise ratio SNR = 27 ± 3 to be reached on each exposure in order to reduce CCD charge transfer inefficiency. Table 2 shows the log of the observations and the corresponding barycentric RVs which were obtained as follow. The spectra were extracted using the SOPHIE pipeline (Bouchy et al. 2009a) and the RVs were measured from the weighted cross-correlation with a G2-type numerical mask (Baranne et al. 1996, Pepe et al. 2002. We excluded the 15 bluer SOPHIE spectral orders from the cross-correlation as they were particularly noisy. Spectra were corrected for CCD charge trans- fer inefficiency (Bouchy et al. 2009b) and RV error bars were computed from the cross-correlation function (CCF) using the method presented by Boisse et al. (2010). The resulting CCFs have full width at half maximum (FWHM) of 10.1 ± 0.1 km/s and contrast that represents ∼ 35 % of the continuum. The final 116-point dataset considered here has been cleared from measurements secured near transit of WASP-148b to avoid any possible deviation due to the Rossiter-McLaughlin effect, as well as a few too low-SNR measurements (RVs less accurate than ±15 m/s have been removed). The HE mode of SOPHIE is known to present possible instrumental drifts (see, e.g., Hébrard et al. 2013). The causes for these drifts are not well understood nor identified but might be due to thermal effects. Following the procedure discussed by Santerne et al (2016), we used the constant star HD 185144 observed the same nights with SOPHIE in HE mode to correct for those potential drifts. The RV dispersion of HD 185144 is 8.5 m/s over the nights where WASP-148 was observed, with a maximum amplitude of 40 m/s observed on a one-month scale in Summer 2015. So that correction allows our results to be significantly improved.
Following the method described e.g. in Pollacco et al. (2008) and Hébrard et al. (2008), we estimated and corrected for the moonlight contamination by using the second SOPHIE fiber aperture, which is targeted on the sky while the first aperture points toward the star. We estimated that 28 spectra of the 116 were significantly polluted by moonlight. On each of them, the moonlight correction ranged from a few to 150 m/s, with a typical value around 30 m/s. Removing these points does not significantly modify the orbital solution.
Our final SOPHIE dataset thus includes 116 RVs having precisions ranging from 6.5 to 14.9 m/s depending on the SNR, with a median value of 9.2 m/s. They are reported in Table 2 and displayed in the lower panel of Fig. 1. The observed 31.0-m/s dispersion is significantly higher than estimated error bars on the measurements, indicating the presence of variability.

Additional photometry
Once the RVs established that the of 8.80-day signal identified with SuperWASP actually was due to a planet, we obtained five extra transit light curves with four different larger ground-based telescopes. They allowed improved spatial and temporal resolutions as well as more precise time-series photometry during   Table 5). The two upper panels show photometric transit data of WASP-148 acquired by a series of observatories (see Table 3). The title of each panel shows the name of the observatory and the date of the observation. The relative flux is provided for each transit. Below each panel, the residuals of the MAP (maximum-a-posteriori) model are shown. For the transit observed by the Sánchez telescope, we present the full dataset in blue empty circles, and a binned version as purple points. The third, lower panel shows the SOPHIE radial velocities and their 1-σ error bars (Table 2). On left they are plotted as a function of time; the residuals to the MAP model are also plotted. The two panels on the right are the phase-folded RV curves for WASP-148b (P = 8.80 d) and WASP-148c (P = 34.5 d) after removing the effect of the other planet. For the transit panels and the phase-folded RVs, the solid black curve is the MAP model, and the gray thin curves are 100 models drawn randomly from the posterior distribution (see Sect. 5).
transits. The goal was to refine the determination of the parameters derived from photometry. These observations are briefly described below. Data reductions are standard and include bias and flat-field corrections, aperture photometry, comparison stars selected to minimise the scatter out of transit, and flux normalization. The images revealed no contamination on the star. The logs of the photometric observations are reported in Table 3, and the five corresponding transit light curves (four complete and one partial transit) are plotted in Fig. 1; they show obvious transit detections.
That table also shows the measured epochs of each transit measured below in Sect. 5. They exhibit significant deviations from the constant ephemerides. We double-checked that the times reported by the observers were correct. In addition, these telescopes and their clocks have been regularly used for other planetary transit studies without showing any timing issue (e.g. Hay et al. 2016, Spake et al. 2016, Demangeon et al. 2018. We thus concluded the WASP-148b transit light curves show hints for TTVs. We did not secure dedicated photometric observations to search for any possible transit of the outer, 34.5-day planet. An inspection of the SuperWASP light curves did not show any significant signature but their time coverage is poor. So we can just conclude here the planet WASP-148c does not show obvious signatures of transits. Table 3. Photometric observations, measured jitters, and transit timing variations of WASP-148b transits. The four last columns summarize the statistics for the marginal posteriors (Sect. 5.1). In particular, the TTVs below are obtained using the mean ephemeris derived in Sect. 5.4.1, the reported error on the TTV amplitude for each transit corresponding to the error on T 0 .

Sánchez
A full transit of WASP-148b was observed on 2015 June 18 at the 1.52-m Carlos Sánchez Telescope (TCS) at the Teide Observatory, Tenerife in the Canary Islands, Spain (Oscoz et al. 2008). We used the Wide-FastCam camera with a Johnson-R filter, and the telescope was manually guided. The mean FWHM through the night was 3.30 pixels and the best aperture for data reduction was 9 pixels.

Mars
One partial transit of WASP-148b was observed at Observatoire Hubert Reeves in Mars, France, without filter on 2017 April 17. A 0.6-m ( f /8) telescope was used with a 2750 × 2200-pixel Atik 460EX camera, and successive exposures of 80-sec durations. The observations and their reduction were done by amateur astronomers, using AudeLa and Muniwin softwares. The telescope has a field of view of 8.68 × 6.9 , and a pixel scale of 0.77 pixel −1 .

Rise
A full transit of WASP-148b was observed on 2017 July 22 with a V+R (og515+kg5) filter using the Rise instrument mounted on the robotic 2.0-m Liverpool Telescope at La Palma in the Canary Islands, Spain. Rise is equipped with a back-illuminated, frametransfer, 1024 × 1024 pixel CCD. The scale is 1.08 pixel −1 and a 2 × 2 binning was used. A total of 2720 6-sec exposures was secured in a row to have a good sampling of the transit. This is the best transit light curve of WASP-148b in the dataset presented here.

Lucky imaging
To further investigate the possibility of stellar contamination of our photometric light curves, on 2016 March 9 we carried out a lucky imaging search for additional companions around WASP-148A using FastCam (Oscoz et al. 2008) on the Carlos Sánchez Telescope at the Teide Observatory, Tenerife, Spain. FastCam has a field of view and pixel scale of 6 and 0.042 pixel −1 , respectively. The detector is an L3 electron-multiplying CCD (EMCCD) with rapid readout, and essentially zero readout noise. We obtained 10 data cubes for WASP-148, each with 1000 images of 50 ms exposure time, using no filter. The data in each cube were bias-subtracted, aligned, and stacked to increase the signal-to-noise ratio. No visible companion objects were found within 6 . So there are no hints here for contamination nor blend.

Stellar properties of WASP-148
The host-star SOPHIE spectra unpolluted by Moonlight were RV corrected and averaged to produce a single spectrum. It was used for spectral analysis using the methods described in Doyle et al. (2013). We used the Hα line to estimate the effective temperature (T eff ), and the Na i D and Mg i b lines as diagnostics of the surface gravity (log g). The iron abundances [Fe/H] were determined from equivalent-width measurements of several clean and unblended Fe i lines and are given relative to the Solar value presented in Asplund et al. (2009). The derived abundance errors include the uncertainties in T eff and log g, as well as the scatter due to measurement and atomic data uncertainties. The projected rotation velocity (v sin i * ) was determined to be around 2 km/s by fitting the profiles of the Fe i lines, after convolving with the SOPHIE instrumental resolution and a macroturbulent velocity of 2.8 ± 0.7 km/s adopted from the calibration of Doyle et al. (2014). We obtained T eff = 5460 ± 130 K, log g = 4.40 ± 0.15 (cgs), [Fe/H] = +0.11 ± 0.08, and logA(Li) < 0.5. Using the calibration of Torres et al. (2010), we derive a stellar mass and radius of 1.00 ± 0.08 M and 1.03 ± 0.20 R , respectively. These values agree with the ones from the Gaia DR2 catalog (Gaia Collaboration 2018).
The SOPHIE spectra show no chromospheric emission peaks in the Ca ii H+K lines, whereas emission would be a sign of stellar activity. We computed the activity index log R HK = −5.09 ± 0.11 from the spectra, and the projected rotational velocity v sin i * = 4.1 ± 1.0 km/s from the parameters of the CCF, both using the calibrations of Boisse et al. (2010). The analysis of spectral lines themselves providing v sin i * ∼ 2 km/s, we finally adopted the conservative value v sin i * = 3 ± 2 km/s. This agrees with the stellar rotation period tentatively detected around 24 days (Sect. 2.1) but the large error bars are not constraining here.
The log R HK value is consistent with the basal limit corresponding to a quiet star. The low RV and photometric jitters found below in Sect. 5.4.2 corroborate this. Some hot Jupiter hosts have apparently depressed values of log R HK , which can be explained by absorption in circumstellar gas ablated from the hot planets (Haswell et al. 2012, 2019, Fossati et al. 2013, Staab et al.2017, but this phenomenon does not appear to have significantly depressed the log R HK value of WASP-148.
From all these analyses, we conclude WASP-148 is a slowly rotating, inactive late-G dwarf.

Radial velocity periodograms
In order to study the periodic signals possibly present in our final SOPHIE dataset, Fig. 2 presents their Lomb-Scargle periodograms (Press et al. 1992) in four different cases, as well as the limits corresponding to false-alarm probabilities of 1 × 10 −3 . The upper panel shows the periodogram of the WASP-148 RVs. Two main significant peaks are clearly detected at periods ∼ 8.80 days and ∼ 34.5 days corresponding to the two planets reported above, together with their weaker aliases around one day.
The second panel shows the periodogram of the RV residuals to a fit including the inner, transiting planet only. The standard deviation of the residuals of this fit is σ O−C = 21.3 m/s, indicating an additional signal is present. On this periodogram, as expected the peak at 8.80 days and its aliases are no longer visible. The main remaining peak is the one at 34.5 days, together with the two fainter double-peaks near 0.97 and 1.03 d corresponding to its aliases with one synodic and one sidereal day due to ground-based sampling. Similarly, the third panel of Fig. 2 shows the periodogram of the RV residuals to a fit including the outer planet only (here again with a large residuals standard deviation σ O−C = 25.3 m/s). The peak at 34.5 days and its aliases are removed, and the 8.80-day signal remains together with the two alias, double-peaks near 0.90 and 1.30 d.
Finally, the bottom panel of Fig. 2 shows the periodogram of the residuals after the two-planet fit shown in Fig. 1. Only lower peaks remain, with false-alarm probabilities below 1 × 10 −3 . We note that in addition to the peaks around one day (which is due to the aliases of all the detected signals), all the four panels in Fig. 2 show possible peaks at longer periods, in particular near 150 days. This could be the signature of a third, outer planet. It would have an RV semi-amplitude around 10 m/s corresponding to a sky-projected mass of ∼ 0.25 M Jup . However, that longperiod signal is not strong enough to claim any detection with the available data. Further observations of this star on a longer time baseline are necessary to establish or not the presence of a third planet in this system.
We also note that none of the periodograms shows any significant power near 24 days, as the one seen in SuperWASP photometry and possibly linked to stellar rotation (Sect. 2.1). Possible RV signals caused by stellar effects are discussed in more details in the following subsection.
In conclusion, this analysis clearly shows that the SOPHIE RVs include a periodic signal at 8.80 days corresponding to the one detected in photometry, together with an additional one at 34.5 days.

Validation of the planetary nature of the RV signals
Here we argue that both RV periods are caused by Doppler shifts due to planets orbiting WASP-148, and not to spectral line profile variations due to stellar activity nor blended binaries.
It is well known that stellar blended configurations can mimic planetary transits, including undiluted eclipsing binaries with low-mass stellar companions or diluted eclipsing binaries (e.g., Almenara et al. 2009). The astrometric excess noise of WASP-148 measured by Gaia is 0.7 mas which is below its detection limit (Kiefer et al. 2019), thus showing no signatures for contamination nor blend caused by a possible massive companion. This agrees with the lack of companion detection in the images obtained for photometry (Sect. 2.3) as well as on our highresolution imaging (Sect. 2.4). In addition, RVs measured using different stellar masks (F0, K0, or K5) produce variations with similar amplitudes to those obtained with the G2 mask, so it is unlikely that these variations are produced by blend scenarios composed of stars of different spectral types.
Similarly, the measured CCF bisector spans (Table 2) quantify possible shape variations of the spectral lines. Indeed, a cor- relation between the bisector and the RV could be the signature of RV variations induced by blend configurations or stellar activity (see, e.g., Queloz et al. 2001, Boisse et al. 2009). Here, bisector spans show a dispersion of 21.9 m/s which is 1.4 times smaller than the RV dispersion, whereas each bisector span is roughly two times less precise than the corresponding RV measurement. Therefore, whereas the RV dispersion is caused by the RV periodic signals, the bisector spans show no significant variations. Moreover, they show no correlations with the RVs. The linear correlation parameter is 0.00 ± 0.07 and the Pearson and Spearman rank correlation factors have low values of -0.01 and -0.06, respectively. This agrees with the conclusion that the RV variations are caused by planetary signals and not by spectralline profile changes attributable to blends or stellar activity.
We made the same tests between bisector spans and RV residuals after fitting WASP-148b alone, WASP-148c alone, or both planets b and c; in none of those cases we detected any correlation. In addition we made the same tests with the FWHMs of the CCF, which as the bisector spans quantify the shape of the spectral lines; here also, we found no correlation between that parameter and the RV nor their residuals. Finally we computed the Lomb-Scargle periodograms of bisector spans and FWHMs; they present no significant periodicities, in particular at the periods of the two signals seen in RVs.
All of this strengthens the inference that the RV variations are not caused by spectral-line profile changes attributable to blends nor stellar activity. We conclude they are caused by exoplanets with orbital periods of 8.80 and 34.5 days. The detected mutual interactions between planets provide an additional argument to this interpretation (see below in Sect.6.4).

Keplerian characterization of the system
Here we present a global Keplerian analysis of our photometric and spectroscopic data. It allows a good fit of them and reliable parameters to be derived for the WASP-148 system. There are compared in Sect. 6 with those obtained from a Newtonian analysis taking into account for planets mutual interactions.

Models and parametrisation
The data are fitted with the models implemented in the pastis package , originally developed to perform statistical validation of transiting candidates. In brief, the RV time series is modeled as a sum of non-interacting Keplerian curves, one per planet. This means that the mutual interactions between planets are neglected. To model the transit light curves pastis implements the JKTEBOP code (Southworth 2011) based on EBOP (Nelson & Davis 1972, Etzel 1981, Popper & Etzel 1981. The Keplerian curves are parametrized by their period P, the semi-amplitude K 1 , and two parameters involving the eccentricity and longitude or pericenter, √ e cos(ω) and √ e sin(ω). For the fifth parameter of the curve, we used the time of transit T c for the inner planet and the mean longitude at epoch BJD = 2 455 500, λ 0 , for the outer one. Finally we included a systemic RV offset, γ. We assumed that the RV residuals are independent and normally distributed around zero. The variance of that distribution is equal to the quadratic sum of the inferred uncertainty for each observation and an additional term, identical for all data points; this constitutes an additional (hyper)parameter of our model, σ 2 RV . The transit light curves are parametrized by the radius ratio k = R p /R s , the stellar density, ρ * , and the impact parameter, b. We chose a linear law for the stellar limb-darkening (LD); this requires an additional wavelength-dependent parameter, the LD coefficient u j , where j = 1, 2, 3 corresponds to the three bandpasses used for the observations (Johnson-R, og515+kg5, and clear). The fluxes of each light curve outside of transit are additional parameters of the model to allow normalization to be adjusted. As for the RVs, the residuals are assumed to be distributed as the sum of two normal distributions, one with width corresponding to the uncertainties of the measurements, and the other one with a variance σ 2 LC i different for each light curve i. The data were modeled assuming the presence of two planets only. Initially we uniquely considered the light curve data, assuming a constant period for the inner planet and that the outer planet did not transit. However, as reported in Sect. 2.3 this did not allow us to explain the ensemble of the light curves which show TTVs. So we implemented a new model in pastis that includes a time shift for each light curve. This means there is an additional model parameter δT c for each light curve. This model accommodated the observations and permitted measuring significant departures from transit times derived from a constant period.
We then added the RVs, modeled as described above, and kept the possibility that each transit light curve presented timing difference to the expected value based on a constant ephemeris. This means that the model has an internal inconsistency, as the departure from perfect Keplerian motion for the planets is assumed to only be reflected in the photometric data. This is justified by the small amplitude of mutual interactions expected for the RVs (see Sect. 5.4.1), as confirmed by our n-body analysis (see Sect. 6.1).
In these fits, the orbital inclination, i p , is not known nor constrained for WASP-148c, as well as the longitude of the ascending node, Ω, for both planets (which is put to 0 here). Some dynamical constraints are put on these parameters however in Sect. 6.3.

Priors
The priors chosen for each parameter are presented in Table 4. These are mostly uninformative priors with some reasonable bounds. Note that as we are not performing model comparison, the exact choice of bounds is not critical here. The parameters that use informative priors are the ephemeris parameters for the transiting planet for which we employed the knowledge from the transit analysis, the normalising flux out of transit for which we chose a normal distribution centered around 1, and the timing offsets used to account for possible timing variations, δT c . For these parameters, we chose a normal distribution centered around zero and with a width of 0.05 day. This choice is equivalent to adding a regularisation term in the likelihood function (Bishop 2007) hindering the values the δT c parameters from becoming extremely large by changing the value of the period or the nominal transit time accordingly.

Posterior sampling
The posterior distribution was sampled using Markov Chain Monte Carlo (MCMC) algorithms. For the first analyses we employed the MCMC algorithm implemented in pastis and described by . This algorithm automatically chooses the parametrisation that minimises correlations between the parameters. To do this, the eigenvectors of the empirical covariance matrix of the parameters are used to define the directions in which the Markov Chain moves in parameter space. [m/s] U(0, 500) U(0, 100) √ e cos ω U(−0.9, 0.9) U(−0.9, 0.9) √ e sin ω U(−0.9, 0.9) U(−0.9, 0.9) Impact parameter, b U(0, 1) −− Radius ratio, R p /R * J(0.01, 0.5) −− Stellar density, ρ * We run the emcee algorithm with 150 walkers for 2 × 10 5 interactions, and we thinned by a factor of 100 because of performance (memory) limitations. We removed the first 10 000 steps, assumed to be the burn-in period. The walkers exhibited adequate mixing and the acceptance rate was centered around 0.15. All the walkers reside in the same region of parameter space, around a maximum of the posterior density. We computed the Geweke (1992) statistics for every parameter and walker, comparing the first 20% of the chain, with successive fractions of the same size. The results are distributed like a normal centered in zero, with a width of around 0.28. All of this indicate that the algorithm is likely converged.
Additionally, we computed the autocorrelation function for the model parameters and selected functions, such as the time of inferior conjunction. We did this for each walker individually in an attempt to identify chains that presented issues. We did not find any problematic walker, and the mean correlation lengths over walkers, defined as the smallest lag for which the autocorrelation is below 1/e, are below 20 steps -i.e. 2000 steps from the unthinned chain -for all parameters (mean: 13.3, median: 12.3), except for √ e b sin ω b and the stellar density, ρ * , for which the mean correlation lengths are around 23.5 and 20.8, respectively. We therefore have a minimum of 6000 independent samples to perform the parameter inference.

Results
The samples obtained from the posterior distribution using the MCMC algorithm allows us to derive the maximum-a-posteriori (MAP) estimate for each parameter, and compute the corresponding model. The MAP model is overplotted on the data in Fig. 1, together with other sample models. The derived parameters are reported in Table 5.

Transit timing variations of WASP-148b
As reported above, the model employed to describe the data allows for a timing offset for each light curve. Based on the posterior samples of the nominal period and transit time, and the individual timing offsets, δT c , we computed the time of inferior conjunction (transit) of WASP-148b for each light curve epoch. Summary statistics for the marginal posterior distributions are reported in Table 3. Figure 3 shows the shape of the marginal posterior distributions. With the exception of the transit times of the SuperWASP light curves, which are the least precise ones (Table 3), the inferred marginal posterior density functions resemble closely a normal distribution. This strengthens the reliability of their timings. With this in mind, we fitted a straight line to the transit times, assuming normal errors equal to the standard deviations reported above. We quadratically added to the error on each transit time a supplementary error representing the typical amplitude of the possible TTVs, assuming they would present sinusoidal variations with time; we show below in Sect. 6.4 and Fig. 9 this actually is the case. We adjusted that supplementary error in order to have a reduced χ 2 = 1 to the straight line and found it is equal to 15.1 min. Prob. Density MARS;2017-04-17 NITES;2014-08-13 NITES;2016-06-13 RISE;2017-07-22 SANCHEZ;2015-06-18 WASP;2008-06-20 WASP;2010-06-04 WASP;2011-05-31 Fig. 3. Marginal posterior distribution for the WASP-148b transit times for each dataset, offset to the mean value of the distribution and normalised using the standard deviation. In most cases, the distributions closely resemble the standard normal distribution, plotted as a black curve.
This leads to a mean ephemeris: where the number between parenthesis represent the error in the parameters and their covariance is cov T (0) c , P = 1.76 × 10 −7 . The residuals of the fit are shown in Fig. 4. This confirms there are significant variations with respect to a constant period. Even removing from the regression the SuperWASP transits (which have large error bars) or the Mars transit (which is partial), the remaining timing measurements exhibit significant TTVs. So the TTVs are significantly detected here.
We consider the mean ephemeris in Equation 1 as our final ones for WASP-148b. There are in agreement with the ephemeris obtained by the fit above considering Keplerian orbits and a time shifts δT c for each light curve.

System parameters
The results from the sampling of the posterior distribution are summarized in Table 5 for some of the model parameters and a number of derived quantities. We also provide the 95%-Highest Density Interval (HDI), defined as the interval containing 95% of the marginal distribution mass, such that all points in the interval have probability densities larger than any point outside. The reported orbital period and time of transit for WASP-148b are the averaged ones obtained from the precedent section. The times for inferior conjunction or transit, T c , and pericenter passage, T p , reported in Table 5 WASP;2008-06-20 WASP;2010-06-04 WASP;2011-05-31 NITES;2014-08-13 SANCHEZ;2015-06-18 NITES;2016-06-13 MARS;2017-04-17 RISE;2017-07-22 100 0 10 0 Fig. 4. Timing residuals of a linear regression model to the inferred WASP-148b transit times. The corresponding constant orbital period is P = 8.803810 ± 0.000043 days as reported in Table 3. The inset is a zoom to the region of the last four transits.
The results for the nuisance model parameters related to the data sets (flux jitter amplitudes) are reported on a separate table (Table 3). The smallest one is obtained for the Rise data and is 70 ± 54 ppm, which could be considered as an upper limit for the flux intrinsic variations of WASP-148. The derived limb-darkening coefficients are 0.93 ± 0.05, 0.75 ± 0.07, and 0.53 ± 0.05 for clear, Johnson-R, and og515+kg5, respectively. This agrees with expected values in those band passes (e.g. Claret et al. 2013). The residuals dispersion of the RVs is 13.9 m/s, which is slightly larger than the typical estimated RV error bars (Sect. 2.2). The RV jitter we fitted to account for that is 11.1 ± 1.4 m/s.
For the derived quantities that required the input from the stellar models, such as the mass of the star for the semi-major axis of the orbits, we sampled values of the required parameter from normal distributions with mean and width corresponding to the values reported in Sect. 3. The transit of WASP-148b allows stellar radius and density to be directly measured to R * = 0.92±0.07 R and ρ * = 1.89±0.48 g/cm 3 . Those values are more accurate but agree with those obtained from stellar analysis in Sect. 3 (R * = 1.03±0.20 R and ρ * = 1.3 +1.2 −0.5 g/cm 3 ). Similarly, the a/R * ratio computed from WASP-148b transits translates to a b = 0.103 ± 0.021 AU, in agreement with the more accurate value a b = 0.0845 ± 0.0022 obtained from M * and the third Kepler law. These agreements reflect the coherence of the fit and its results.
WASP-148b is a hot Jupiter with a mass M b = 0.291 ± 0.025 M Jup and a radius R b = 0.722 ± 0.055 R Jup , translating into a bulk density of 0.95 g/cm 3 . WASP-148c has a sky-projected Table 5. WASP-148 system parameters. The stellar parameters are obtained from Sect. 3. The transit and orbital parameters are obtained from the Keplerian fit statistics of the joint posterior sample (Sect. 5). In particular, the mean ephemeris of WASP-148b are obtained from the TTV fit presented in see Sect. 5.4.1. The maximum a posteriori probability estimates are given together with the standard deviations of the marginal distribution following the ± sign. For some parameters, the extremes of the 95%-Highest Density Interval are also given in brackets. Except the semi-amplitude K 1 which refers to the star, the other orbital parameters (ω, λ 0 , T 14 , b, a...) refer to the orbit of the planets.  i.e. ω b − ω c = 0 is within the 95%-HDI, but the uncertainties remain too large for a significant determination. Figure 5 presents the histograms of the marginal distributions of the eccentricities and of ω b − ω c . We also report in Table 5  perature T eq of both planets. It is computed while the planet is at at the semi-major axis of the star, assuming black body, a Bond albedo of 0.1, and uniform heat redistribution to the planetary night side. We find T eq = 940 ± 80 K and 590 ± 50 K for WASP-148b and WASP-148c, respectively. The 13.9-m/s dispersion of the RV residuals might include the signature of additional planets, as the possible third planet discussed in Sect. 4.1. A three-planet fit does not significantly modify the parameters of WASP-148b and WASP-148c.
Several studies have shown that the Keplerian signature of a single eccentric planet might also be fit with a model including two planets on circular (or nearly circular) orbits in 2:1 resonance, depending on the RV accuracy and their time sampling (e.g., Anglada-Escudé et al. 2010, Wittenmyer et al. 2013, Kürster et al. 2015. So we attempted to fit our dataset with threeplanet models, thus including WASP-148b and two outer planets on circular orbits around 17.3 and 34.5 days instead of the eccentric orbit of WASP-148c presented above. The quality of that fit is similar to the two-planet fit. Except its eccentricity, the resulting 34.5-d planet as properties similar to those of WASP-148c presented above, in particular its mass. The fitted, circular, 17.3d planet produces a semi-amplitude of 9±2 m/s corresponding to a Neptune mass. Such a signal is at the limit of detection according the accuracy of our RV data, and the periodograms presented in Fig. 2 show no hints for a signal around 17 days. So the presence of the planet WASP-148c on a 34.5-d orbit is not put into doubt, but the possibility that its orbit is circular and another planet is present on a 17.3-d, circular orbit between WASP-148b and WASP-148c could not be totally excluded from the available RVs. However, we favor here the solution with two planets only, both on eccentric orbits. Section 6 below and the presence of TTVs give additional arguments in favor of eccentric orbits.

Dynamic analyses
The orbital solution given in Table 5  In addition, the ratio between the orbital periods, P c /P b = 3.92, is close to a 4:1 MMR. As a consequence, mutual gravitational interactions between planets are likely to be significant. We study and quantify those dynamical aspects here, with particular focus on the instabilities which may rise, and their effects on the TTVs.

N-body characterization of the planetary system
The solution given in Table 5 was obtained assuming noninteracting, Keplerian orbits. We first performed an n-body Newtonian fit to the RV data taking into account planets mutual interactions. We performed that fit using the approach presented by Correia et al. (2010) and fixing the epoch of WASP-148b transits from Table 5.
The dispersion of the resulting RV residuals is 14.1 m/s. This is similar to the dispersion of the resulting RV residuals of the Keplerian fit (13.9 m/s, see Sect. 5.4.2), so both fits are equivalent on this point of view. The obtained orbital parameters show no significant differences with those obtained to the Keplerian ones. This justifies the hypotheses made in Sect. 5 and the reliability of the parameters reported in Table 5 over time scales of a few years, i.e. that of our datasets.
We study below the dynamical effects on longer time scales.

Stability analysis
We performed a global frequency analysis (Laskar 1990(Laskar , 1993 in the vicinity of the best fit (Table 5), in the same way as achieved for other planetary systems (e.g. Correia et al. 2005. This allows us to analyse and estimate the stability of the orbital solution, Fig. 6. Stability analysis of the WASP-148 planetary system, assuming coplanar orbits. For fixed initial conditions, the phase space of the system is explored by varying the semi-major axis a c (x-axis) and eccentricity e c (y-axis) of the outer planet WASP-148c. The step size is 10 −2 in eccentricity and 10 −3 in semimajor axis. For each initial conditions, the system is integrated over 50 kyr and a stability criterion is derived with the frequency analysis of the mean longitude. The chaotic diffusion is measured by the variation in the frequencies. The red zone corresponds to highly unstable orbits, while the dark-blue region can be assumed to be stable on a billion-years timescale. The reduced-χ 2 level curves of the Newtonian fit are also plotted.
The system is integrated on a regular 2D mesh of initial conditions, with varying semi-major axis and eccentricity of WASP-148c, while the other parameters are retained at their nominal values. We used the symplectic integrator SABA1064 of Farrès et al. (2013) with a step size of 5 × 10 −3 yr and general relativity corrections. Each initial condition is integrated over 50 kyr, and a stability indicator is computed to be the variation in the measured mean motion over the two consecutive 25 kyr intervals of time (for more details, see Couetdic et al. 2010). For regular motion, there is no significant variation in the mean motion along the trajectory, while it can vary significantly for chaotic trajectories. Figure 6 shows the wide vicinity of the nominal solution, together with the reduced-χ 2 level curves (whose minimum gives the best fitted solution of the Newtonian fit). The stability indicator is reported using a color index, where the red zones represent the strongly chaotic trajectories, and the dark-blue zones the extremely stable ones. We observe the presence of the large 4:1 MMR and its chaotic separatrix. We can see that the present system is outside this resonance, in a more stable area (dark region). We hence conclude that the WASP-148 planetary system is stable.
We also directly tested the stability of the MAP solution from Table 5 by performing a numerical integration over 1 Gyr. As expected, the orbits evolve in a regular way, and remain stable throughout the simulation. Nevertheless, the nominal solution is close to an unstable region (Fig. 6) due to the high eccentricity of the outer planet. This suggests the WASP-148c eccentricity might be slightly smaller to bring the system to an even more stable region. Overall, that analysis allows further constrains to be put on the planetary parameters, by reducing the region of parameter space where the orbits are stable.
In addition we performed a frequency analysis of the orbital solution computed over 100 kyr, assuming coplanar orbits. The fundamental frequencies of the systems are the mean motions n b and n c , and the two secular frequencies of the pericenters g 1 and g 2 ( Table 6). Because of the proximity of the two innermost orbits, there is a strong coupling within the secular system (see Laskar 1990). Both planets WASP-148b and c precess with the same frequency, g 1 . The two pericenters are thus locked and ∆ω = ω b − ω c oscillates around 0 • (aligned ellipses), with a maximal amplitude of about 45 • (Fig. 7). The secular period for the eccentricity and ∆ω oscillations is 2π/(g 2 − g 1 ) = 881 yr. Thus, on this time scale the eccentricities of WASP-148b and c roughly oscillate between 0.1 and 0.45, and 0.25 and 0.4, respectively. Whereas most hot Jupiters have circularized orbits due to tidal dissipation (e.g. Dawson & Johnson 2018), the eccentricity of WASP-148b could be at least partially explained by its interactions with WASP-148c.  Fig. 7. Secular evolution of the WASP-148 system, assuming coplanar orbits (Table 5). We show the eccentricity of WASP-148b (top) and WASP-148c (middle), and the angle ∆ω = ω c − ω b (bottom). Those three parameters oscillate with a 881-yr period. Table 6. Fundamental frequencies for the nominal orbital solution in Table 5. n b and n c are the mean motions, and g 1 and g 2 are the secular frequencies of the pericenters. We assumed here coplanar orbits. ig. 8. Global stability analysis of the WASP-148 planetary system. We fix all orbital parameters of the solution shown in Table 5 and vary the only two unconstrained parameters, both from the outer planet WASP148c: the longitude of its ascending node, Ω c (in degrees, on the x-axis), and its inclination, i c (in degrees, on the y-axis). The upper panel shows the whole (i c , Ω c ) domain whereas the lower panel zooms on the most stable regions. The step size is 0.2 • in the node and 0.5 • in the inclination. For each initial condition, the system is integrated over 50 kyr and a stability criterion is derived with the frequency analysis of the mean longitude. Due to dynamical invariant, the figure is symmetric with respect to the (i c = 90 • , Ω c = 0) center. White, dashed curves give the isolines of constant mutual inclination I = 10 • , 20 • , 30 • , 40 • . The red zones correspond to highly unstable orbits, while the dark-blue region can be assumed to be stable on a billion-years timescale. The reduced-χ 2 level curves of the Newtonian fit are also plotted in white.

Additional constraints
Since WASP-148b transits its host star, we are able to determine its inclination to the line of sight i b = 89.80 • ± 0.27 • ( Table 5). As a result, the system is left with only two undetermined parameters: the orbital inclination of the outer planet, i c , and the difference between the longitude of the ascending nodes, ∆Ω = Ω c − Ω b . The longitude of the ascending node of WASP-148b can be fixed at any value, so for simplicity we fix Ω b = 0 • , and ∆Ω = Ω c . We can thus build a 2D stability map for the two unknown parameters, to see how dynamics can constrain their possible values. Figure 8 explores the stability in the (i c , Ω c ) domain, by keeping the remaining parameters fixed at their values shown in Table 5. We also show the reduced-χ 2 level curves. They present a minimum for i c = 73 • , Ω c = 26 • (and by symmetry i c = 107 • , Ω c = −26 • ), but the contour levels do not put strong constraints on the determination of this minimum with the available data. This is why we are unable to determine at present these parameters from our Newtonian fit (see, e.g., . However, there is only a subset of (i c , Ω c ) values for which the system can be stable: 55 • i c ≤ 125 • and |Ω c | 35 • (Fig. 8, lower panel). We note there exists also a stable zone for retrograde orbits (around Ω c = 180 • ; Fig. 8, upper panel) but it is excluded as RVs indicate a prograde orbit of WASP-148c by comparison to WASP-148b. These dynamical constraints have direct consequences on the determination of the true mass of WASP-148c, M c , despite the absence of transit detection. Using the 95%-HDI upper limit M c sin i c < 0.49 M Jup (Table 5), we conclude that M c ≤ 0.60 M Jup . This is a stringent upper limit on the true mass of an exoplanet detected from RV only.
Another constraint can be derived for the mutual inclination between orbital planes, I. Assuming i b = 90 • , we have cos I = sin i c cos Ω c . (2) In Fig. 8 we plot the lines of constant mutual inclination, which describe circles around (i c = 0 • , Ω c = 0 • ). We conclude that all stable areas correspond to I 35 • , so the orbital plan of both planets could not have a mutual inclination larger than this value.

Transit-timing variations
As shown in Sects. 2 and 5.4.1, apparently only WASP-148b does transit in front of the host star, showing significant TTVs (Sect. 5.4.1). These TTVs are likely to be produced by gravitational interactions with WASP-148c and can be used to constrain its orbit, in particular the eccentricity or the true mass (e.g. Lissauer at al. 2011). We model here those TTVs.
Using the MAP solution from Table 5 and fixing Ω c = 0 • , we generated the TTVs for two different configurations: one with i c = 90 • , corresponding to a coplanar system, and another with i c = 60 • , corresponding to a mutual inclination I = 30 • . In the upper panel of Fig. 9 we show the variations corresponding to each solution. Over our ten-year-span observations, they present a sinusoidal shape with a period of about 460 days and an amplitude of about two hours. This is the expected shape in such configuration of two planets in MMR, and agrees with the assumption made in Sect. 5.4.1. Figure 9 also shows the measured TTVs as reported in Table 3. The two-hour amplitude of the computed TTVs is somewhat larger than the observations. The modest agreement could be partially due to the fact that we do not fit the TTVs simultaneously with the other data, but instead measure the TTVs then fit them. It could also be partially due to the uncertainty in some of the orbital parameters, in particular the eccentricities and the arguments of the pericenter. The lower panel of Fig. 9 presents the predicted TTVs assuming lower eccentricities, within the 95%-HDI. Here they have amplitudes in better agreement with the observations, in particular the model assuming eccentricities of 0.11 and 0.21 for WASP-148b and c, respectively. The precision and the number of photometric measurements currently available for the WASP-148 system, together with the fact TESS will soon observe it (Sect. 7), do not justify to run an exhaustive search for a best fit solution, and hence reduce the uncertainty in the parameters (e, ω) and constrain the unknown parameters (i c , Ω c ).
Finally, we note that if we assume both orbits are circular, the amplitude of the WASP-148b TTVs would be negligible, on the order of a few seconds or below. So the detection of TTVs supports the fact that the orbits are eccentric (Sect. 5.4.2). And if WASP-148c eventually is discovered as transiting its host star, it should also present TTVs. They would be anti-correlated with  Table 5 with i c = 90 • , Ω c = 0 • (coplanar orbits, in red), or with i c = 60 • , Ω c = 0 • (mutual inclination I = 30 • , in blue). The dots correspond to the observed TTVs (Table 3). For comparison, the lower panel presents curves corresponding to coplanar solutions with lower eccentricities, within our 95%-HDI: e b = 0.18, e c = 0.29 (in orange) or e b = 0.11, e c = 0.21 (in green).
those of WASP-148b but with a larger amplitude due its longer orbital period, as observed e.g. for Kepler-9.

Conclusions
We have presented the discovery and characterization of the WASP-148 exoplanetary system. This is based on ten years of photometric and spectroscopic observations and their Keplerian and Newtonian analyses. The system includes a 0.29-M Jup , 0.72-R Jup hot Jupiter transiting its star each 8.80 days, and an outer planet apparently not transiting with a period of 34.5 days and a sky-projected mass of 0.40-M Jup (true mass below 0.60 M Jup ). The planetary equilibrium temperatures are 940 K and 590 K, respectively. The orbits of both planets are eccentric and have a mutual inclination below 35 • . They present significant gravitational interactions due to their period ratio near the 4:1 meanmotion resonances. That orbital configuration is stable but shows significant deviations from purely Keplerian orbits. In particular, the inner planet exhibits transit timing variations of about 15 minutes. That configuration makes WASP-148 a unique case. Among systems with the most similarities with WASP-148, but still with significant differences, one can cite Kepler-9, Kepler-277, or TOI-216. Kepler-9 (Holman et al. 2010) is briefly described above in Sect. 1. Kepler-277 hosts two transiting planets with radii about three times the Earth radius and orbital periods of 17.3 and 33.0 days (Wu & Lithwick 2013, Xie 2014). That configuration causes TTVs of a few minutes implying masses similar to those of WASP-148 planets. The Kepler-277 planets are thus particularly dense but their masses remain poorly known as the eccentricities of the orbits are not measured. TOI-216 hosts two transiting giant planets with periods of 17.1 and 34.6 days (Kipping et al. 2019, Dawson et al. 2019). TTVs of a few minutes are detected and allow masses of 0.56 and 0.08 M Jup to be respectively evaluated. Here the eccentricities and the planetary densities are low. Those three systems show TTVs and are near the 2:1 MMR. Their planets are located in the "period valley" a domain of orbital periods between 10 and 100 days known to be sparse in giant planets (e.g. Udry et al. 2003, Santerne et al 2016. Their MMR configurations might be linked to the fact they are present in the valley.
WASP-148 is in a similar configuration but here close to the 4:1 MMR and WASP-148b share properties with standard hot Jupiters. A significant difference between WASP-148 and the three systems above is the eccentricity of its planets, which here are significantly different from zero. This is also a significant difference of WASP-148b with other hot Jupiters, whose orbits are circular for most of them.
The NASA Transiting Exoplanet Survey Satellite (TESS) secures an all-sky photometric survey to detect planetary transits in front of bright stars (Ricker et al. 2015). Its camera #2 will observe WASP-148 on its Sectors #24, #25, and #26 from 2020 April 16 to July 4. It will thus observe nine consecutive transits of WASP-148b, i.e. significantly increase the number of available transits and dramatically improve their accuracy. This will improve the transit parameters, and more importantly, it should allow the TTVs to be confirmed and refined. Together with the fact that they will be obtained consecutively and with a single instrument, the higher quality of those new transit light curves will allow the dynamical model of the WASP-148 system to be refined, and in particular a full TTV analysis to be done. In addition, TESS will cover three inferior conjunctions of WASP-148c (on 2020 April 23, May 27, and July 1). This could reveal this outer planet does transit its host star (and also present TTVs), or not. If WASP-148b and c are coplanar, there is a 99.4-% probability for WASP-148c to transit. If the mutual inclination between both orbits is I = 1 • or 5 • , that probability falls to 63 % and 4 %, respectively.
The future TESS 79-day continuous observation may also reveal additional planets in the system at short or long periods, in particular small-size planets which now could not be significantly detected in our RV data but could be found in such accurate, space-based light curve. We will also continue our RV follow-up of the system in order to refine the system parameters, in particular the planet eccentricities, and to possibly confirm and characterize the possible long-period planet that shows a detection hint in our current RV dataset. Spectroscopic observations of a transit of WASP-148b could also be secured to measure its obliquity (e.g. Winn et al. 2009, Hébrard et al. 2011. The amplitude of the Rossiter-McLaughlin anomaly could be around 10 m/s or larger, so reachable by several spectrographs.