© G. Hébrard et al. 2020

1 Introduction

Extrasolar planets that transit their host stars are especially interesting. When they are characterized 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 to the exoplanet archives¹. A few were first discovered from radial-velocity (RV) surveys, and photometric follow-up subsequently revealed their transits (e.g., Charbonneau et al. 2000; Motalebi et al. 2015). The vast majority were first identified from photometric surveys and were then characterized with RV follow-up, however. 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 orbital eccentricity and obliquity and to characterize the host stars. Long-term RV follow-up could also reveal additional nontransiting planets in the system (e.g., Christiansen 2017; Rey et al. 2018). Multiplanetary systems like this are particularly interesting for the studies of their dynamics.

Dynamics can also be studied in multiplanetary systems when transit-timing variations (TTVs) are detected. 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 might 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 multiplanetary systems, TTVs have larger amplitudes when the orbital periods of the planets are nearly commensurable, that is, in or near mean-motion resonances (MMR). This makes the TTV analysis a powerful technique for characterizing such systems, and in particular, for measuring planetary masses and eccentricities, or even detecting additional perturbing planets (e.g., Nesvorny et al. 2012).

For years, several attempts have been made to detect TTVs of transiting planets with ground-based photometry (e.g., Díaz et al. 2014a; Maciejewski et al. 2010), but most of them later were not confirmed (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 that are 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. This star is transited by two giant planets on ~ 19.2 and ~ 38.9 days each. This almost 2:1 MMR causes TTVs with an amplitude of about one day that are clearly detected with Kepler.

Today, a few dozen exoplanets have been detected or characterized based on their TTVs. Most of these detections were made based on Kepler photometry. This includes multiplanetary transiting systems as well as single-transiting planets showing TTVs that allow the detection and characterization of additional nontransiting planets. Notable cases include KOI-142b, a P = 10.95 d planet showing TTVs of up to one day that are caused by a nontransiting giant planet with a period that is twice longer (Nesvorny et al. 2013) and was eventually detected in RVs (Barros et al. 2014), or the seven Earth-size planets orbiting the star Trappist-1, which have periods between 1.5 and 18.7 days and TTVs up to a few dozen minutes (Gillon et al. 2017).

Most 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: the first researches were mainly attempted on this type of system. Only 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; Patra et al. 2017a,b; Yee et al. 2020). Several scenarios have been proposed to explain these TTVs, whose amplitudes are lower than 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 that are explained by 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 of a hot Jupiter showing TTVs, and the first interacting system with TTVs detected from the ground. 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 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 that have been observed with ground-based telescopes reveal significant deviations from a constant orbital period, with TTV amplitudes of about 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, and Sect. 6 discusses dynamic analyses of the system before we conclude in Sect. 7.

2 Observations and data reduction

2.1 Photometric identification with SuperWASP

Located on La Palma in the Canary Islands, Spain, SuperWASP-North consists of eight Canon 200 mm f∕1.8 focal lenses coupled to e2v 2048 × 2048 pixel CCDs with 13.7′′ pixels and a field of view of 7.8° × 7.8°, associated with a custom-built photometric reduction pipeline (Pollacco et al. 2006). It observes with a broadband 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 (Cameron et al. 2006).

With 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 h, 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 SuperWASPlight curves are shown in the three first plots of Fig. 1; the transit-like features are not obvious in these initial data, which illustrates how sensitive the candidate detection algorithms should be.

2.2 Radial-velocity follow-up with SOPHIE

After its identification from SuperWASP photometry, we started an RV follow-up of WASP-148 with the SOPHIE spectrograph at the 1.93 m telescope of the 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, 2013). 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 observation season (35 SOPHIE measurements over six months) revealed RV variations in phase with the 8.80-day signal detected with SuperWASP, no significant variations in the spectral line profiles, and an amplitude of the RV variations of about 30 m s⁻¹. This showed that the transiting body is a planet slightly more massive than Saturn, designated 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 of about 35 days, here again with no significant variations in the spectral line profiles and an RV variation amplitude of about 30 m s⁻¹. This indicated the star hosts a second outer planet that is slightly more massive than Saturn, hereafter designated 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 s with a typical value of 1200 s. This allowed us to reach a nearly constant signal-to-noise ratio (S∕N) = 27 ± 3 on each exposure in order to reduce CCD charge transfer inefficiency. Table 2 shows the observation log and the corresponding barycentric RVs, which were obtained as follows. 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-transfer 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 a full width at half maximum (FWHM) of 10.1 ± 0.1 km s⁻¹, and the contrast represents ~35% of the continuum. The final 116-point dataset considered here was cleared from measurements secured near the transit of WASP-148b to avoid any possible deviation due to the Rossiter-McLaughlin effect, as well as a few measurements whose S/N was too low (RVs less accurate than ± 15 m s⁻¹ were 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 or identified, but might be due to thermal effects. Following the procedure discussed by Santerne et al. (2016), we used the constant star HD 185144 that was observed on the same nights with SOPHIE in HE mode to correct for these potential drifts. The RV dispersion of HD 185144 is 8.5 m s⁻¹ on the nights where WASP-148 was observed, with a maximum amplitude of 40 m s⁻¹ observed on a one-month scale in July-August 2015. This correction allowed a significant improvement of our results.

Following the method described in Pollacco et al. (2008) and Hébrard et al. (2008), for example, we estimated and corrected for the moonlight contamination 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. In each of them, the moonlight correction ranged from a few to 150 m s⁻¹, with a typical value of about 30 m s⁻¹. Removing these points does not significantly modify the orbital solution.

Our final SOPHIE dataset thus includes 116 RVs with precisions ranging from 6.5 to 14.9 m s⁻¹ depending on the S/N, 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 the estimated error bars on the measurements, indicating variability.

Fig. 1 Photometric and RV observations of WASP-148, together with their two-planet Keplerian fit described in Sect. 5 (the parameters are reported in Table 6). The two upper panels show photometric transit dataof WASP-148 acquired by a series of observatories (see Table 3). The title of each panel shows the observatory name and observation date. The relative flux is provided for each transit. Below each panel, the residuals of the MAP model are shown. For the transit observed by the Sánchez telescope, we present the full dataset as blue empty circles, and a binned version as purple points. The lower set of panels shows the SOPHIE RVs and their 1σ error bars (Table 2). At left they are plotted as a function of time; the residuals to the MAP model are also plotted. The two panels at the right are the phase-folded RV curves for WASP-148b (P = 8.80 d) and WASP-148c (P = 34.5 d) after the effect of the other planet is removed. In 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).

2.3 Additional photometry

After the RVs established that the of 8.80-day signal identified with SuperWASP was indeed due to a planet, we obtained five additional 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 transits. The goal was to refine the determination of the parameters derived from photometry. These observations are briefly described below. Data reductions were standard and include bias and flat-field corrections, aperture photometry, comparison stars selected to minimize the scatter out of transit, and flux normalization. The images revealed no contamination on the star. The logs of the photometricobservations are reported in Table 3, and the five corresponding transit light curves (four complete transits and one partial transit) are plotted in Fig. 1; they show obvious transit detections.

Table 3 also shows the measured epochs of each transit measured below in Sect. 5. The epochs exhibit significant deviations from the constant ephemerides. We verified 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 problems (e.g., Hay et al. 2016; Spake et al. 2016; Demangeon et al. 2018). We therefore concluded the WASP-148b transit light curves show indications of 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. We can therefore only conclude here the planet WASP-148c does not show obvious signatures of transits.

2.3.1 Nites

Two full transits of WASP-148b were observed with the Near Infra-red Transiting ExoplanetS (NITES) Telescope: afirst transit without filter on 2014 August 13, and a second transit with a Johnson-R filter on 2016 June 13. NITES, located at La Palma in the Canary Islands, Spain, is a semirobotic, 0.4-m (f∕10) Meade LX200GPS Schmidt-Cassegrain telescope (McCormac et al. 2014). It is mounted with a Finger Lakes Instrumentation Proline 4710 camera and a 1024 × 1024 pixel deep-depleted CCD made by e2v. The telescope has a field of view of 11′ ×11′ and a pixel scale of 0.66′′ pixel⁻¹.

2.3.2 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.

2.3.3 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 s durations. The observations and their reduction were made by amateur astronomers, using the AudeLa and Muniwin softwares. The telescope has a field of view of 8.68′ ×6.9′, and a pixel scale of 0.77′′ pixel⁻¹.

2.3.4 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 frame-transfer 1024 × 1024 pixel CCD. The scale is 1.08′′ pixel⁻¹, and a 2 × 2 binning wasused. A total of 2720 6 s 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.

2.4 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-148 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⁻¹. The detector is an L3 electron-multiplying CCD (EMCCD) with rapid readout, and essentially zero readout noise. We obtained ten 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 S/N. The resulting data do not allow us to quantify themagnitude contrast well, but we can conclude that no visible companion objects were found within 6′′ of the star. This means that there are no indications here for contamination or blend.

3 Stellar properties of WASP-148

3.1 Rotation periods

We used the sine-wave fitting method described in Maxted et al. (2011) to search for quasi-periodic modulation in the SuperWASP light curves of WASP-148 caused by the combination of the stellar rotation and magnetic activity, that is, star spots. Variability due to star spots is not expected to be coherent on long timescales as a consequence of the finite lifetime of star spots and differential rotation in the photosphere. We therefore analyzed each season of data for WASP-148 separately. We also analyzed the data from each camera that was used to observe WASP-148 separately so that we could assess the reliability of the results. Only combinations of cameras and seasons with more than 2000 observations were included in this analysis. We removed the transit signal from the data prior to calculating the periodograms by subtracting a simple transit model from the light curve. We calculated periodograms over 8192 uniformly spaced frequencies from 0 to 0.5 cycle day⁻¹. The false-alarm probability was calculated using the bootstrap Monte Carlo method described in Maxted et al. (2011).

The results of this analysis are given in Table 4. The best-fit period for all the light curves obtained with camera 143 are consistent with a rotation period of about 26 days or its first harmonic. This supposed rotation period is not detected in the data obtained with camera 141. This may be because the data from this camera are more strongly affected by instrumental noise than the data from camera 143. The mean period from the data obtained with camera 143 is 26.2 ± 1.3 days. This is the value we adopted for the stellar rotation period P_rot. The amplitude of this long-term modulation is of a few millimagnitudes and therefore is not visible in the light curves plotted in Fig. 1, which focus on short-duration transits alone.

3.2 Spectral characterization

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. (2016). 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 that is due to measurement and atomic data uncertainties. The projected rotation velocity (v sin i_*) was determined to be about 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 derived a stellar mass and radius of 1.00± 0.08M_⊙ and 1.03± 0.20R_⊙, respectively. These values agree with the values of T_eff = 5350 ± 115 K and R_* = 0.97 ± 0.04R_⊙ from the Gaia DR2 catalog (Gaia Collaboration 2018). As a sanity check, T_eff was also obtained from the spectral energy distribution (SED). This was obtained using broad-band photometry from Table 1, except for the three from Gaia Collaboration (2018) in order to be independent of the DR2 T_eff value. The photometry was converted into fluxes and the best-fitting (Kurucz 1993) model flux distribution was determined, which gave a value of T_eff = 5430 ± 140 K, which again agrees well (see Fig. 2).

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 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 provided v sin i_* ~ 2 km s⁻¹, therefore we finally adopted the conservative value v sin i_* = 3 ± 2 km s⁻¹. This is consistent with the stellar rotation period and radius reported above, but the large error bars are not constraining here.

The 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. The values of of some hot-Jupiter hosts apparently lie below this basal limit, which is inconsistent with our understanding of the behavior of late-type stars. This might be explained by absorption in circumstellar gas that is ablated from the hot planets. This gas will absorb in cores of the Ca II H+K lines (Haswell et al. 2012, 2019; Fossati et al. 2013; Staab et al. 2017). This circumstellar absorption phenomenon may operate in the WASP-148 system, but it has not decreased the value below the basal limit. The intrinsic activity of WASP-148 might consequently be slightly above that indicated by taken at face value, but there is nothing to suggest that this effect is significant here as the star is quiet. Based on all these analyses, we conclude that WASP-148 is a slowly rotating inactive late-G dwarf.

Fig. 2 Fitted spectral energy distribution of the star WASP-148 (gray line) overplotted with measured magnitudes (black circles, Table 1).

4 Evidence for planets orbiting WASP-148

4.1 Radial velocity periodograms

In order to study the periodic signals that might be present in our final SOPHIE dataset, Fig. 3 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⁻³. 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 alone. The standard deviation of the residuals of this fit is σ_O-C = 21.3 m s⁻¹, indicatingthat an additional signal is present. In this periodogram, the peak at 8.80 days and its aliases are no longer visible, as expected. The main remaining peak is the peak at 34.5 days, together with the two fainter double-peaks near 0.97 and 1.03 d that correspond to its aliases with one synodic and one sidereal day due to ground-based sampling. Similarly, the third panel of Fig. 3 shows the periodogram of the RV residuals to a fit including the outer planet alone (here again with a large residuals standard deviation σ_O-C = 25.3 m s⁻¹). The peak at34.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. 3 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⁻³. We note that in addition to the peaks around one day (which are due to the aliases of all the detected signals), all the four panels in Fig. 3 show possible peaks at longer periods, in particular near 150 days. This might be the signature of a third outer planet. It would have an RV semiamplitude of aobut 10 m s⁻¹ that corresponds to a sky-projected mass of ~ 0.25 M_Jup. However, this long-period signal is not strong enough to claim any detection with the available data. Further observationsof this star on a longer time baseline are necessary to establish or discard the presence of a third planet in this system.

We also note that none of the periodograms shows any significant power near 26 days, as the one seen in SuperWASP photometry and possibly linked to stellar rotation (Sect. 3.1). Possible RV signals caused by stellar effects are discussed in more detail in the following subsection.

In conclusion, this analysis clearly shows that the SOPHIE RVs include a periodic signal at 8.80 days. This corresponds to the signal that was detected in photometry, together with an additional signal at 34.5 days.

Fig. 3 Lomb-Scargle periodograms of the SOPHIE RVs of WASP-148. Upper panel: periodogram computed based on the initial RVs without any fit removed. Second and third panels: periodograms computed based on the residuals of the fits including WASP-148b or WASP-148c, respectively. Bottom panel: periodogram after subtraction of the two-planet fit. The horizontal dashed lines correspond to the false-alarm probability of 1 × 10−3.

4.2 Validation in planetary nature of the RV signals

Here we argue that both RV periods are caused by Doppler shifts that are caused by planets orbiting WASP-148, and not by spectral line profile variations that are due to stellar activity or 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 or a 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 in our high-resolution 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, therefore 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. A correlation between the bisector and the RV might indeed 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 half as precise as the corresponding RV measurement. This means that whereas the RV dispersion is caused by the periodic RV 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 spectral-line profile changes that are 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 these cases did we detect 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 we found no correlation between this parameter and the RV or their residuals either. 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 or stellar activity. We conclude that they are caused by exoplanets with orbital periods of 8.80 and 34.5 days. The detected mutual interactions between the planets provide an additional argument for this interpretation (see below in Sect. 6.4).

5 Keplerian characterization of the system

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

5.1 Models and parameterization

The data were fit with the models that are implemented in the pastis package (Díaz et al. 2014b), which was originally developed to perform statistical validation of transiting candidates. In brief, the RV time series is modeled as a sum of noninteracting 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 parameterized by their period P, the semiamplitude K₁, and two parameters involving the eccentricity and longitude of pericenter, and . 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, λ₀, 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 this 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, .

The transit light curves were parameterized 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 the normalization to be adjusted. The RV residuals were assumed to be distributed as the sum of two normal distributions, one with a width corresponding to the measurement uncertainties, and the other one with a variance 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 that shows TTVs. We therefore implemented a new model in pastis that included 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 a 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 was set to 0 here). Some dynamical constraints are placed on these parameters in Sect. 6.3, however.

5.2 Priors

The priors chosen for each parameter are presented in Table 5. These are mostly uninformative priors with some reasonable bounds. Because we did not perform a 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 normalizing 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 regularizing term in the likelihood function (Bishop 2007) and prevents the values the δT_c parameters from becoming extremely high by changing the value of the period or the nominal transit time accordingly.

5.3 Posterior sampling

The posterior distribution was sampled using Markov chain Monte Carlo (MCMC) algorithms. For the first analyses we employed the MCMC algorithm that is implemented in pastis and described by Díaz et al. (2014b). This algorithm automatically chooses the parameterization that minimizes 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. For the final run with the full data set, we employed the ensemble sampler described by Goodman & Weare (2010) that is implemented in the emcee package by Foreman-Mackey et al. (2013).

We ran the emcee algorithm with 150 walkers for 2 × 10⁵ interactions, and we thinned by a factor of 100 because of performance (memory) limitations. We removed the first 10 000 steps, which we 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 and compared the first 20% of the chain, with successive fractions of the same size. The results are distributed like a normal centered on zero, with a width of around 0.28. All of this indicates that the algorithm likely has 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 problems. 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 and the stellar density, ρ_*, for which the mean correlation lengths are about 23.5 and 20.8, respectively. We therefore have a minimum of 6000 independent samples on which to perform the parameter inference.

5.4 Results

The samples obtained from the posterior distribution using the MCMC algorithm allowed 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 6.

5.4.1 Transit-timing variations of WASP-148b

As reported above, the model we 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 on 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 4 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 transit times (Table 3), the inferred marginal posterior density functions closely resemble a normal distribution. This strengthens the reliability of their timings.

With this in mind, we fit a straight line to the transit times, assuming normal errors equal to the standard deviations reported above. We quadratically added a supplementary error to the error on each transit time that represents the typical amplitude of the possible TTVs, assuming they would present sinusoidal variations with time; we show below in Sect. 6.4 and Fig. 10 that this is the case. We adjusted this supplementary error in order to have a reduced χ² = 1 to the straight line and found that it is equal to 15.1 min.

This leads to a mean ephemeris (1)

where the number between parentheses represents the error on the parameters, and their covariance is .

The residuals of the fit are shown in Fig. 5. This confirms that there are significant variations with respect to a constant period. Even when we remove the SuperWASP transits (which have large error bars) or the Mars transit (which is partial) from the regression, the remaining timing measurements exhibit significant TTVs. This means that the TTVs are significantly detected here.

We consider the mean ephemeris as our final ephemeris for WASP-148b. It agrees with the ephemeris obtained by the fit above considering Keplerian orbits and a time shift δT_c for each light curve.

Fig. 4 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.

Fig. 5 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 into the region of the last four transits.

5.4.2 System parameters

The results from the sampling of the posterior distribution are summarized in Table 6 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 higher than any point outside. The reported orbital period and time of transit for WASP-148b are the averaged periods and times obtained from the precedent section. The times for inferior conjunction or transit, T_c, and pericenter passage, T_p, reported in Table 6 are those measured near our most accurate transit (the time for pericenter passage is observed with Rise). For WASP-148c, these parameters are slightly more accurate when they are measured at an epoch located in the middle of the SOPHIE observations: T_c = 2 456 968.77 ± 0.90 BJD and T_p = 2 456 965.0 ± 1.3 BJD.

The results for the nuisance model parameters related to the data sets (flux jitter amplitudes) are reported in a separate table (Table 3). The smallest nuisance model parameter 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 LD 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 these band passes (e.g., Claret et al. 2013). The residual dispersion of the RVs is 13.9 m s⁻¹, which is slightly higher than the typical estimated RV error bars (Sect. 2.2). The RV jitter we fit to account for this 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 semimajor axis of the orbits, we sampled values of the required parameter from normal distributions with a mean and width corresponding to the values reported in Sect. 3.2. The transit of WASP-148b allows stellar radius and density to be directly measured as R_* = 0.92 ± 0.07R_⊙ and ρ_* = 1.89 ± 0.48 g cm⁻³. These values are more accurate but agree with those obtained from stellar analysis in Sect. 3.2 (R_* = 1.03 ± 0.20R_⊙ and g cm⁻³). Similarly, the a∕R_* ratio computed from WASP-148b transits translates into 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.055R_Jup, translating into a bulk density of 0.95 g cm⁻³. WASP-148chas a sky-projected mass of M_c sini_p = 0.397 ± 0.044 M_Jup. The orbits of both planets are significantly eccentric, with the MAP estimate of planet b and c being e_b = 0.22 ± 0.06 and e_c = 0.36 ± 0.09. Based on the determination of the ω angles, the pericenters might be aligned, that is, ω_b − ω_c = 0 is within the 95%-HDI, but the uncertainties remain too large for a significant determination. Figure 6 presents the histograms of the marginal distributions of the eccentricities and of ω_b − ω_c. We also report in Table 6 the equilibrium temperature T_eq of both planets. It was computed while the planet was at the semimajor axis of the star, assuming blackbody, a Bond albedo of 0.1, and uniform heat redistribution to the planetary nightside. We find T_eq= 940 ± 80 K and 590 ± 50 K for WASP-148band WASP-148c, respectively.

The 13.9-m s⁻¹ dispersion of the RV residuals might include the signature of additional planets, such 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). We therefore attempted to fit our dataset with three-planet 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 this fit is similar to the two-planet fit. Except for its eccentricity, the resulting properties of the 34.5 d planet are similar to those of WASP-148c presented above, in particular its mass. The fit circular 17.3 d planet produces a semiamplitude of 9 ± 2 m s⁻¹, which corresponds to a Neptune mass. This signal is at the detection limit according to the accuracy of our RV data, and the periodograms presented in Fig. 3 show no indication for a signal at about 17 days. The planet WASP-148c on a 34.5 d orbit is not placed in doubt, therefore, 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 cannot be entirely excluded from the available RVs. However, we favor here the solution with only two planets, both on eccentric orbits. Section 6 below and the TTVs give additional arguments in favor of eccentric orbits.

Fig. 6 Histogram of the marginal posterior samples of the orbital eccentricities of WASP-148b and WASP-148c (top) and of the difference in the argument of pericenter (bottom).

6 Dynamical analyses

The orbital solution given in Table 6 shows a compact system (a_b = 0.084 and a_c = 0.210 AU) with eccentric orbits (e_b = 0.22 and e_c = 0.36). The stability of the system is not straightforward because the planetary masses are on the same order as the mass of Saturn (M_b = 0.29 M_Jup and M_c sini_p = 0.40 M_Jup). 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 these dynamical aspects here, with particular focus on the instabilities that may arise, and their effects on the TTVs.

6.1 N-body characterization of the planetary system

The solution given in Table 6 was obtained assuming noninteracting Keplerian orbits. We first performed an n-body Newtonian fit to the RV data taking the mutual interactions of the planets into account. We performed this fit using the approach presented by Correia et al. (2010) and fixing the epoch of WASP-148b transits from Table 6.

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), which means that both fits are equivalent from this point of view. The obtained orbital parameters show no significant differences with those obtained for the Keplerian parameters. This justifies the hypotheses made in Sect. 5 and the reliability of the parameters reported in Table 6 over timescales of a few years, which is that of our datasets. We study the dynamical effects on longer timescales below.

6.2 Stability analysis

We performed a global frequency analysis (Laskar 1990, 1993) in the vicinity of the best fit (Table 6), in the same way ashas been achieved for other planetary systems (e.g., Correia et al. 2005, 2010). This allows us to analyse and estimate thestability of the orbital solution.

The system was integrated on a regular 2D mesh of initial conditions, with varying semimajor axes and eccentricity of WASP-148c, while the other parameters were retained at their nominal values. We used the symplectic integrator SABA1064 of Farrès et al. (2013) with a step size of 5 × 10⁻³ yr and general relativity corrections. Each initial condition was integrated over 50 kyr, and a stability indicator was derived with the frequency analysis of the mean longitude, 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 7 shows the wide vicinity of the nominal solution, together with the reduced-χ² level curves (whose minimum gives the best-fit 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 show the extremely stable trajectories. We observe the large 4:1 MMR and its chaotic separatrix. This 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 6 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. 7) because the outer planet is highly eccentric. This suggests the eccentricity of WASP-148c might be slightly lower to bring the system to an even more stable region. Overall, this analysis allows further constraints to be placed on the planetary parameters by reducing the region of parameter space in which 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₁ and g₂ (Table 7). Because of the proximity of the two orbits, there is a strong coupling within the secular system (see Laskar 1990). Planets WASP-148b and c both precess with the same frequency, g₁. The two pericenters are thus locked and Δω = ω_b − ω_c oscillates around 0° (aligned ellipses), with a maximum amplitude of about 45° (Fig. 8). The secular period for the eccentricity and Δω oscillations is 2π∕(g₂ − g₁) = 881 yr. On this timescale, the eccentricities of WASP-148b and c therefore roughly oscillate between 0.1 and 0.45, and 0.25 and 0.4, respectively. Whereas most hot Jupiters have circularized orbits as a result of tidal dissipation (e.g., Dawson & Johnson 2018), the eccentricity of WASP-148b might at least partially be explained by its interactions with WASP-148c.

Fig. 7 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 semimajor axis ac and eccentricity ec 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 color scale corresponds to values between − 9 (black) and − 3 (red) for the decimal logarithm of the stability index D used in Correia et al. (2010). The red zone corresponds to highly unstable orbits, while the dark blue region can be assumed to be stable on a billion-year timescale. The reduced-χ2 level curves of the Newtonian fit are also plotted.

6.3 Additional constraints

Because 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 6). 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, therefore for simplicity, we fixed Ω_b = 0° and ΔΩ = Ω_c. We thus built a 2D stability map for the two unknown parameters to determine how dynamics can constrain their possible values.

Figure 9 explores the stability in the (i_c, Ω_c) domain by keeping the remaining parameters fixed at the values shown in Table 6. We also show the reduced-χ² 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 place strong constraints on the determination of this minimum with the available data. We are therefore unable to determine these parameters from our Newtonian fit at present (see, e.g., Correia et al. 2010). However, there is only one subset of (i_c, Ω_c) values for which the system can be stable: 55° ≲ i_c ≤ 125° and |Ω_c|≲ 35° (Fig. 9, lower panel). We note that a stable zone for retrograde orbits also exists (around Ω_c = 180°; Fig. 9, upper panel), but it is excluded because the 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 sini_c < 0.49 M_Jup (Table 6), 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 alone.

Another constraint can be derived for the mutual inclination between orbital planes, I. Assuming i_b = 90°, we have (2)

In Fig. 9 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 that the orbital plane of the two planets cannot have a mutual inclination higher than this value.

Fig. 8 Secular evolution of the WASP-148 system, assuming coplanar orbits (Table 6). We show the eccentricity of WASP-148b (top) and WASP-148c (middle), and the angle Δω = ωc − ωb (bottom). These three parameters oscillate with an 881 yr period.

Fig. 9 Global stability analysis of the WASP-148 planetary system. We fixed all orbital parameters of the solution shown in Table 6 and varied the only two unconstrained parameters, both from the outer planet WASP148c: the longitude of its ascending node, Ωc, and its inclination, ic. The upper panel shows the whole (ic, Ωc) domain, whereas the lower panel zooms ino the most stable regions. The step size was 0.2° in the node and 0.5° in the inclination. For each initial condition, the system was integrated over 50 kyr and a stability criterion was derived with the frequency analysis of the mean longitude. Because of the dynamical invariant, the figure is symmetric with respect to the (ic = 90°, Ωc = 0) center. White dashed curves give the isolines of constant mutual inclination I = 10°, 20°, 30°, and 40°. The color scale corresponds to values between − 9 (black) and − 3 (red) for the decimal logarithm of the stability index D used in Correia et al. (2010). The red zones correspond to highly unstable orbits, while the dark blue region can be assumed to be stable on a billion-year timescale. The reduced-χ2 level curves of the Newtonian fit are also plotted in white.

6.4 Transit-timing variations

As shown inSects. 2 and 5.4.1, whereas transits of WASP-148c have not been detected nor ruled out with the existing data, WASP-148b does transit in front of the host star and shows significant TTVs. These TTVs are likely to be produced by gravitational interactions with WASP-148c and can be used to constrain itsorbit, in particular, the eccentricity or the true mass (e.g., Lissauer et al. 2011). We model here those TTVs.

Using the MAP solution from Table 6 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. 10 we show the variations corresponding to each solution. Over our ten-year 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 a configuration of two planets like this near MMR, and it agrees with the assumption made in Sect. 5.4.1.

Figure 10 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 might partially be explained as follows: we did not fit the TTVs simultaneously with the other data, but instead measured the TTVs and then fit them. It might 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. 10 presents the predicted TTVs assuming lower eccentricities, within the 95% HDI. Here their amplitudes agree better 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 that TESS will soon observe it (Sect. 7), do not justify running 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 when we assume that both orbits are circular, the amplitude of the WASP-148b TTVs would be negligible, on the order of a few seconds or below. The detection of TTVs therefore supports the assumption that the orbits are eccentric (Sect. 5.4.2). If WASP-148c were eventually be discovered to be transiting its host star, it should also present TTVs. They would be anticorrelated with those of WASP-148b, but have a larger amplitude because its orbital period is longer, which has also been observed for Kepler-9, for example.

Fig. 10 Transit-timing variations for the WASP-148b planet. In the upper panel, the curves correspond to numerical simulations for the solution shown in Table 6 with ic = 90°, Ωc = 0° (coplanar orbits, in red), or with ic = 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: eb = 0.18, ec = 0.29 (in orange) or eb = 0.11, ec = 0.21 (in green).

7 Conclusions

We have presented the discovery and characterization of the WASP-148 exoplanetary system. This is based on ten years of photometricand 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 every 8.80 days, and an outer planet that is 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 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 MMR. This orbital configuration is stable, but shows significant deviations from purely Keplerian orbits. In particular, the inner planet exhibits TTVs of about 15 minutes.

This configuration makes WASP-148 a unique case. As systems with the greatest similarities with WASP-148, but still with significant differences, we 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 that are about three times the Earth radius and orbital periods of 17.3 and 33.0 days (Wu & Lithwick 2013; Xie 2014). This configuration causes TTVs of a few minutes, which implies masses similar to those of the 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 evaluated, respectively. Here the eccentricities and the planetary densities are low. These three systems show TTVs and are near the 2:1 MMR. Their planets are located in the period valley, which is a domain of orbital periods between 10 and 100 days that is 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 shares 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 to other hot Jupiters, most of whose orbits are circular.

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 in its sectors 24–26 from 2020 April 16 to July 4. It will thus observe nine consecutive transits of WASP-148b, which means that it will significantly increase the number of available transits and dramatically improve their accuracy. This will improve the transit parameters, and more importantly, it is expected to 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 these new transit light curves will allow the dynamical model of the WASP-148 system to be refined, and in particular, a full TTV analysis can be achieved. In addition, TESS will cover two inferior conjunctions of WASP-148c (on 2020 May 1 and June 4). This might reveal that this outer planet does transit its host star (and also presents TTVs), or it might 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°, this 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 that so far could not be significantly detected in our RV data but may be found in an 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 planetary 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 may 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 might be about 10 m s⁻¹ or larger, which means that it can be reached by several spectrographs.

Acknowledgements

We thank the Observatoire de Haute-Provence (CNRS) staff for its support. This work was supported by the “Programme National de Planétologie” (PNP) of CNRS/INSU, and the CNRS-PICS program between France and Argentina (PICS07826). A.C.M.C. acknowledges support by CFisUC strategic project (UID/FIS/04564/2019), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal. O.D.S.D. acknowledges support from FCT and FEDER/COMPETE2020 in the form of the work contract DL 57/2016/CP1364/CT0004 and the projects EPIC (PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953) and GEANES (PTDC/FIS-AST/32113/2017 - POCI-01-0145-FEDER-032113). This work was granted access to the HPC resources of MesoPSL financed by the Région Île-de-France and the project Equip@Meso of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche (ANR-10-EQPX-29-01). This work is partly financed by the Spanish Ministry of Economics and Competitiveness through project ESP2016-80435-C2-2-R and PGC2018-098153-B-C31. This article is based on observations made with different facilities, including the SOPHIE spectrograph at the 1.93-m telescope of Observatoire Haute-Provence, France, the SuperWASP-North instrument located on La Palma in the Canary Islands, Spain, the Near Infra-red Transiting ExoplanetS (NITES) Telescope located at La Palma in the Canary Islands, Spain, the FastCam and Wide-FastCam instruments, at Telescopio Carlos Sánchez operated on the island of Tenerife by the IAC in the Spanish Observatorio del Teide, and the Observatoire Hubert Reeves in Mars, France.

References

All Tables

All Figures

Fig. 1 Photometric and RV observations of WASP-148, together with their two-planet Keplerian fit described in Sect. 5 (the parameters are reported in Table 6). The two upper panels show photometric transit dataof WASP-148 acquired by a series of observatories (see Table 3). The title of each panel shows the observatory name and observation date. The relative flux is provided for each transit. Below each panel, the residuals of the MAP model are shown. For the transit observed by the Sánchez telescope, we present the full dataset as blue empty circles, and a binned version as purple points. The lower set of panels shows the SOPHIE RVs and their 1σ error bars (Table 2). At left they are plotted as a function of time; the residuals to the MAP model are also plotted. The two panels at the right are the phase-folded RV curves for WASP-148b (P = 8.80 d) and WASP-148c (P = 34.5 d) after the effect of the other planet is removed. In 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).

In the text

	Fig. 2 Fitted spectral energy distribution of the star WASP-148 (gray line) overplotted with measured magnitudes (black circles, Table 1).
In the text

Fig. 3 Lomb-Scargle periodograms of the SOPHIE RVs of WASP-148. Upper panel: periodogram computed based on the initial RVs without any fit removed. Second and third panels: periodograms computed based on the residuals of the fits including WASP-148b or WASP-148c, respectively. Bottom panel: periodogram after subtraction of the two-planet fit. The horizontal dashed lines correspond to the false-alarm probability of 1 × 10−3.

In the text

	Fig. 4 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.
In the text

	Fig. 5 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 into the region of the last four transits.
In the text

	Fig. 6 Histogram of the marginal posterior samples of the orbital eccentricities of WASP-148b and WASP-148c (top) and of the difference in the argument of pericenter (bottom).
In the text

Fig. 7 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 semimajor axis ac and eccentricity ec 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 color scale corresponds to values between − 9 (black) and − 3 (red) for the decimal logarithm of the stability index D used in Correia et al. (2010). The red zone corresponds to highly unstable orbits, while the dark blue region can be assumed to be stable on a billion-year timescale. The reduced-χ2 level curves of the Newtonian fit are also plotted.

In the text

	Fig. 8 Secular evolution of the WASP-148 system, assuming coplanar orbits (Table 6). We show the eccentricity of WASP-148b (top) and WASP-148c (middle), and the angle Δω = ωc − ωb (bottom). These three parameters oscillate with an 881 yr period.
In the text

Fig. 9 Global stability analysis of the WASP-148 planetary system. We fixed all orbital parameters of the solution shown in Table 6 and varied the only two unconstrained parameters, both from the outer planet WASP148c: the longitude of its ascending node, Ωc, and its inclination, ic. The upper panel shows the whole (ic, Ωc) domain, whereas the lower panel zooms ino the most stable regions. The step size was 0.2° in the node and 0.5° in the inclination. For each initial condition, the system was integrated over 50 kyr and a stability criterion was derived with the frequency analysis of the mean longitude. Because of the dynamical invariant, the figure is symmetric with respect to the (ic = 90°, Ωc = 0) center. White dashed curves give the isolines of constant mutual inclination I = 10°, 20°, 30°, and 40°. The color scale corresponds to values between − 9 (black) and − 3 (red) for the decimal logarithm of the stability index D used in Correia et al. (2010). The red zones correspond to highly unstable orbits, while the dark blue region can be assumed to be stable on a billion-year timescale. The reduced-χ2 level curves of the Newtonian fit are also plotted in white.

In the text

Fig. 10 Transit-timing variations for the WASP-148b planet. In the upper panel, the curves correspond to numerical simulations for the solution shown in Table 6 with ic = 90°, Ωc = 0° (coplanar orbits, in red), or with ic = 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: eb = 0.18, ec = 0.29 (in orange) or eb = 0.11, ec = 0.21 (in green).

In the text