Photodynamical analysis of the nearly resonant planetary system WASP-148: Accurate transit-timing variations and mutual orbital inclination

WASP-148 is a recently announced extra-solar system harbouring at least two giant planets. The inner planet transits its host star. The planets travel on eccentric orbits and are near the 4:1 mean-motion resonance, which implies significant mutual gravitational interactions. In particular, this causes transit-timing variations of a few minutes, which were detected based on ground-based photometry. This made WASP-148 one of the few cases where such a phenomenon was detected without space-based photometry. Here, we present a self-consistent model of WASP-148 that takes into account the gravitational interactions between all known bodies in the system. Our analysis simultaneously fits the available radial velocities and transit light curves. In particular, we used the photometry secured by the TESS space telescope and made public after the WASP-148 discovery announcement. The TESS data confirm the transit-timing variations, but only in combination with previously measured transit times. The system parameters we derived agree with those previously reported and have a significantly improved precision, including the mass of the non-transiting planet. We found a significant mutual inclination between the orbital planes of the two planets: I=41.0 +6.2 -7.6 deg based on the modelling of the observations, although we found I=20.8 +/- 4.6 deg when we imposed a constraint on the model enforcing long-term dynamical stability. When a third planet was added to the model - based on a candidate signal in the radial velocity - the mutual inclination between planets b and c changed significantly allowing solutions closer to coplanar. We conclude that more data are needed to establish the true architecture of the system. If the significant mutual inclination is confirmed, WASP-148 would become one of the only few candidate non-coplanar planetary systems. We discuss possible origins for this misalignment.


Introduction
While the orbit of a single planet around its host star is well reproduced by a Keplerian model with a constant orbital period, multi-planetary systems have more complex orbits. Indeed, the mutual gravitational interactions between planets imply small deviations from the Keplerian orbits and in particular slight variations of the orbital periods. Such effects are negligible for most systems, but they are amplified when orbital periods are exactly or nearly commensurable. When such a resonant or nearly resonant system includes transiting planets, their orbital periods can be accurately measured, allowing their variations to be detected. This was predicted in particular by Holman & Murray (2005) and Agol et al. (2005), who showed how interactions in multiplanetary systems might cause transit-timing variations (TTVs).
When they are measured, TTVs are a powerful tool to characterise planetary systems. In particular, they allow constraints to be put on masses as well as on the presence of additional planets. This was done on the first system in which TTVs were detected (Holman et al. 2010) as well as on dozens of other detections reported thereafter (e.g. Lissauer et al. 2011;Nesvorný et al. 2013;Wang et al. 2014;Hadden & Lithwick 2014;Gillon et al. 2017;Freudenthal et al. 2019;Jontof-Hutter et al. 2021). Among different methods employed to analyse TTV observations, the pho-todynamical modelling 1 (Carter et al. 2011) of light curves can be used together with radial velocities to constrain planetary system parameters without using external inputs on the masses or radii of the host stars often derived from stellar evolution models (Agol et al. 2005). One of the strengths of this technique is that it fits the full transit light curves (including transit durations) instead of only using the transit timings. Photodynamical analyses have been used to characterise several TTV systems (e.g. Almenara et al. 2018a,b). Here, we apply it to the system WASP-148 recently reported by Hébrard et al. (2020).
Most of the confirmed TTV detections to date have been discovered using light curves obtained from space telescopes (namely Kepler or Spitzer). On the other hand, the TTVs in the WASP-148 system were detected using ground-based telescopes only. The WASP-148 system includes (at least) two giant planets. The inner one, WASP-148 b, is a hot Saturn of 0.72±0.06 R J and 0.29 ± 0.03 M J that transits its host with an orbital period of 8.80 days. The outer planet, WASP-148 c, has an orbital period of 34.5 days and a minimum mass of 0.40 ± 0.05 M J and a true mass < 0.60 M J . It was discovered and characterised using radial velocities obtained with the SOPHIE spectrograph. No tran-sits of this planet were detected. The orbits of both planets have significant eccentricities (e b = 0.22 ± 0.06 and e c = 0.36 ± 0.09) and their orbital periods fall near the 1:4 mean-motion resonance. This particular configuration induces amplified dynamical effects, and TTVs were detected with an amplitude of about ±15 minutes and a period of roughly 460 days (Hébrard et al. 2020).
Several analyses of the available data sets have been presented by Hébrard et al. (2020). The first one simultaneously fitted transit light curves and radial velocities, but it did not include any mutual interactions. The model allowed for small, artificial, ad hoc shifts in the transit times in order to reproduce the TTVs. The second model was fitted to the radial velocities taking mutual interactions into account, constrained by the average period and phase derived from the transit light curves. Both analyses gave similar results and the derived system properties imply large TTVs for both planets, whereas only those of WASP-148 b could actually be observed as no transits of WASP-148 c were detected. However, Hébrard et al. (2020) did not present a complete model fitted to both radial velocity and transit light curve data sets that takes the gravitational interactions between the planets into account. Finally, Hébrard et al. (2020) indicate that there is a hint of a possible third planet in the system with a period near 150 days and a minimum mass around 0.25 M J , without confirming that however.
Following the detection and characterisation of the WASP-148 system by Hébrard et al. (2020), observations of that star secured with the Transiting Exoplanet Survey Satellite (TESS) were released. TESS provides continuous, high-quality photometry over 28 days (Ricker et al. 2015). In the case of WASP-148, it observed seven new transits of WASP-148 b, and provided a unique opportunity to search for possible transits of WASP-148 c or potential additional planets in the system. Maciejewski et al. (2020) published an analysis of the TTVs and radial velocities of WASP-148 presented by Hébrard et al. (2020) in addition to TESS data, as well as two new ground-based transit observations. They took into account the gravitational interactions between the planets, but they used only the transit timings instead of the whole transit light curves. Recently, Wang et al. (2022) have shown the orbit of WASP-148 b is aligned and prograde.
Here, we present new analyses of the WASP-148 system, applying the full photodynamical approach on the available transit light curves and radial velocities. The article is organised as follows: In Section 2 we describe the data used; in Section 3 we determine the stellar parameters; in Sections 4 and 5 we analyse the radial velocity and transits without accounting for the gravitational interactions between the planets, respectively; in Section 6 we detail the photodynamical modelling; and in Section 7 we present the results. Finally, we discuss the results of our work in Section 8.

Observations
We used the ground-based photometry and SOPHIE radial velocities presented in Hébrard et al. (2020). We converted the time of the photometry from BJD UTC to BJD TDB 2 (Eastman et al. 2010), except for the Telescopio Carlos Sánchez transit, and NITES observation on June 13, 2016, which were already in BJD TDB . In addition to the data in Hébrard et al. (2020), we added the photometry from TESS and four transits observed with the 1.5 m Ritchey-Chrétien Telescope at the Sierra Nevada Ob-servatory (OSN150), the first two of which have been presented in Maciejewski et al. (2020).

TESS
TESS observed WASP-148 in sectors 24, 25, and 26, with a total time span of 79.3 days (TIC 115524421, TOI-2064). Two transits were lost during the interruption of the observations that occur every TESS orbit (around 14 days) at the middle of sectors 25 and 26. Each TESS sector lasts two orbits of the satellite. Around the perigee of the TESS orbit, data collection is paused. In total, seven new transits of WASP-148 b were observed. The photometry is available in the full-frame images (FFIs) at 30-minute cadence. The FFIs were calibrated by the Science Processing Operations Center (SPOC) at NASA Ames Research Center (Jenkins et al. 2016). We used eleanor (Feinstein et al. 2019) to extract the light curve of WASP-148 from the TESS FFIs. We chose the point spread function photometry that has the lower dispersion for this object. eleanor corrects the times for the object coordinates, which otherwise are set to the centre of the CCD in FFIs. There is a star 4.5 mag fainter in the Gaia G-band at 26 of WASP-148 that partially contaminates the TESS photometry 3 . This is taken into account in our modelling (Sections 5 and 6). The TESS data show no evidence for any transit of WASP-148 c, which is in agreement with Maciejewski et al. (2020).

OSN150
Two new precise photometric time series for transits of WASP-148 b were acquired in March and June 2021 using the 1.5 m Ritchey-Chrétien telescope (OSN150) at the Sierra Nevada Observatory (OSN, Spain). The instrument was equipped with a Roper Scientific VersArray 2048B CCD camera with a 2048 × 2048 × 13.5 µm back-illuminated matrix. The field of view was 7 .92 × 7 .92 with the pixel scale of 0 . 232 per pixel. The instrument was mildly defocussed to allow for longer exposure times and a lower fraction of time lost for CCD readout. The observations were gathered without any filter to increase the signalto-noise ratio for transit-timing purposes. The observing runs were scheduled following an ephemeris from Maciejewski et al. (2020). The telescope was auto-guided to keep the star at the same position in the CCD matrix. The details on the individual runs are given in Table 1. The observations started about 90 minutes before a transit ingress and lasted about 90 minutes after an egress. This out-of-transit monitoring was secured for detrending purposes. On March 22, 2021 the observations were stopped about 12 minutes after the transit due to dawn.
Photometric data reduction was performed with the AstroIm-ageJ software (Collins et al. 2017) following a standard calibration procedure. The science frames were de-biased and flat-field calibrated using sky flat frames. The light curves were generated with the differential aperture photometry method. The aperture size and a collection of comparison stars were set after a series of test runs to minimise the data point scatter. The fluxes were de-trended against airmass and time using out-of-transit data only and then they were normalised to unity outside the transits. Timestamps were converted into barycentric Julian dates in barycentric dynamical time BJD TDB . For the homogeneity of our   Notes. Date UT is given for the beginning of an observing run. X shows the changes of the target's airmass during a run. N obs is the number of useful scientific exposures. t exp is the exposure time. Γ is the median number of exposures per minute. pnr is the photometric noise rate (Fulton et al. 2011) in parts per thousand (ppth) of the normalised flux per minute of the observation.
analysis, we also reprocessed the two transit light curves from Maciejewski et al. (2020) following the same de-trending procedure.

Stellar parameters
To determine the stellar parameters of WASP-148, we used the precise parallax determination by Gaia (Gaia Collaboration et al. 2016, which was not exploited by Hébrard et al. (2020). Stellar atmosphere models and stellar evolution models are also required to model the observed spectral energy distribution (SED). We constructed the SED of WASP-148 using the magnitudes from Gaia Early Data Release 3 (Gaia EDR3, Riello et al. 2021), the 2-Micron All-Sky Survey (2MASS, Skrutskie et al. 2006;Cutri et al. 2003), and the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010;Cutri & et al. 2013).
The measurements are listed in Table A.1. We modelled these magnitude measurements using the procedure described by Díaz et al. (2014), with informative priors for the effective temperature (T eff ), surface gravity (log g), and metallicity ([Fe/H]) from Hébrard et al. (2020), and for the distance from Gaia EDR3 (Lindegren et al. 2021). We used non-informative priors for the rest of parameters. The priors are listed in Table 2. We decided to use an additive jitter for each set of photometric bands (Gaia, 2MASS, and WISE), which had the effect of slightly broadening the posteriors of the stellar parameters. We used the two stellar atmosphere models, PHOENIX/BT-Settl (Allard et al. 2012) and ATLAS/Castelli & Kurucz (Castelli & Kurucz 2003), and two stellar evolution models, Dartmouth (Dotter et al. 2008) and PARSEC (Chen et al. 2014). We obtained posterior samples for the four combinations of stellar atmosphere models and stellar evolution models using the Markov Chain Monte Carlo (MCMC) algorithm from Díaz et al. (2014). The posteriors of the stellar parameters of the individual combinations agree within 1-σ (see Figure 1). We merged the results assuming an equal probability for each model combination (labelled as 'merged' in Figure 1). The posteriors' median and 68.3% credible intervals (CI) for model parameters and for derived physical quantities of interest are listed in Table 2. Those results agree with the ones reported by Hébrard et al. (2020), but they are more accurate, mainly due to the use of the parallax from Gaia.
The data with the maximum a posteriori (MAP) stellar atmosphere model is shown in Figure 2. Before Gaia EDR3, we carried out the same analysis with Gaia DR2 (Gaia Collaboration et al. 2018), obtaining similar results. In particular, we obtained the same stellar radius error, despite a factor ∼3 increase in the precision of the parallax. Thus, at least for WASP-148, the precision in the parallax is not the limiting factor to improve the stellar radius determination with the SED technique.
We used the stellar rotation period, P rot = 26.2±1.3 days, derived in Hébrard et al. (2020) and the stellar mass derived in this section (Table 2) to estimate a gyrochronological age, neglecting the influence of the planets, of 4.0 +0.9 −0.7 Gyr (Barnes 2010; Barnes & Kim 2010, using a P 0 between 0.12 and 3.4 days), where we added a systematic 10% error to the statistical one (Meibom et al. 2015). The isochronal age (Table 2) agrees with the gyrochronological age within the uncertainties, but it is less precise.

SOPHIE-only analysis
When radial velocity and transit observations are modelled simultaneously (as was done by Hébrard et al. 2020), both data sets constrain the eccentricity. In this section, we aim for a planetary eccentricity determination using only radial velocity data. We analysed the SOPHIE radial velocities with a two-planet Keplerian model, without taking mutual gravitational interactions into account and using juliet (Espinoza et al. 2019), radvel (Fulton et al. 2018), and dynesty (Speagle 2020). We used normal priors from Hébrard et al. (2020) for the period and the time of inferior conjunction of planet b, and non-informative  priors for the rest of the parameters. We included a jitter parameter that was added quadratically to the velocity uncertainties. We performed a second analysis using a Gaussian process (GP) regression model with a quasi-periodic (QP) kernel to model the error terms. We employed the QP kernel included in celerite (Foreman-Mackey et al. 2017), with noninformative priors for the hyperparameters. We found an odds ratio of ln Z jitter /Z jitter+QP = 0.5 ± 0.4. Neither of the models is strongly favoured over the other (Kass & Raftery 1995), and we therefore adopted the simpler model of uncorrelated error terms with an additive jitter. The posterior median and 68.3% credible interval are shown in Table 3. They agree with the results presented by Hébrard et al. (2020).
We investigated the significance of the main peak in the residuals of the two-planet model ( Fig. 3 in Hébrard et al. 2020), thus we repeated the analysis adding a third Keplerian with a period of ∼150 days. We found a minimum mass of 0.24 ± 0.05 M J , similar to the 0.25 M J value reported by Hébrard et al. (2020). The model comparison gives an odds ratio of ln Z 3 planets /Z 2 planets = 1.6 ± 0.4, which is positive evidence in favour of the three-planet model but below the strong evidence cutoff to be preferred over the two-planet model, so we continued with the adopted two-planet model. Notes. The parameters are: Orbital period, time of conjunction, eccentricity, argument of pericentre, radial velocity semi-amplitude, systemic velocity, and jitter.

Transit-only analysis
We derived the transit times of planet b without considering planet c or the radial velocities. This modelling is independent of the two-planet system hypothesis adopted in the photodynamical modelling (Section 6), and it can be used to verify if it is an appropriate hypothesis given the data.
In Hébrard et al. (2020), the error term of the transit photometry data were modelled using a simple additive jitter term. However, if systematics are not correctly accounted for, the posteriors of the transit modelling can be biased (Barros et al. 2013). This should be particularly severe for incomplete transits. Here, we reanalyse the transits presented in Hébrard et al. (2020) and the TESS observations using a more sophisticated error model.
To model the transits of planet b simultaneously with the systematics, we used juliet (Espinoza et al. 2019), using batman (Kreidberg 2015) for the transit model, and we chose the approximate Matern kernel GP included in celerite (Foreman-Mackey et al. 2017). We used different GPs for each groundbased transit, and for each TESS sector. The timing of each individual transit is a free parameter. We used a prior for the stellar density from Section 3, and for √ e cos ω and √ e sin ω from Hébrard et al. (2020). Otherwise, we adopted non-informative or large priors for the rest of the parameters. We oversampled (Kipping 2010a) the model of the MARS transit by a factor 3, which has an 88 second cadence, and the TESS data by a factor 30, which have about a 30 minute cadence. The remaining transit observations have cadences between 0.9 and 37 seconds, and we deemed them unnecessary to oversample their model. We set a dilution factor for the wide-field observations: WASP (one dilution for the three WASP transits), and TESS (one dilution per sector). To reduce the number of free parameters, we made the choice to use one set of limb darkening parameters for all transit observations without a filter ('clear'), but they certainly do not correspond to exactly the same instrument bandpass. In total, there are 76 free parameters. To sample from the posterior, we used the nested sampling code dynesty (Speagle 2020). The complete list of parameters, priors, and posteriors are shown in Table A.2. Figure 3 shows the data and the model posterior. We compare the timing of the ground-based observations with the results in Hébrard et al. (2020), and found differences within 1σ, except for the partial transit observed with MARS where it is 12.0 ± 2.5 min later (see Figure 5). The transit normalisation can affect the transit-timing determination (Barros et al. 2013). The MARS transit can be more affected because it is a partial transit and the baseline is not well-defined.
We repeat this analysis for the four OSN150 transits, which were included at a later stage in this work. The results are presented in Table A.3 and Figure 4.
This analysis could be used to de-trend the transits for the photodynamical modelling. However, in this process, one set of parameters of the transit model needs to be fixed, including a transit time. Therefore, if the de-trended transits are then modelled, the results can be biased. Instead, we decided to include the GP in the photodynamical modelling.

Photodynamical modelling
While Sections 4 and 5 do not take mutual interactions into account, here we report our fits of the observed photometry and radial velocity measurements accounting for the gravitational interactions between the three bodies known in the system using a photodynamical model. Its positions and velocities in time were obtained through an n-body integration. The sky-projected po-    Figure 3, but for the OSN150 transits.
sitions were used to compute the light curve (Mandel & Agol 2002) using a quadratic limb-darkening law (Manduca et al. 1977), which we parametrised following Kipping (2013). To account for the integration time, the model was oversampled by a factor of 30 and 3 for the TESS and MARS data, respectively, and then binned back to match the cadence of the data points (Kipping 2010a). The line-of-sight projected velocity of the star issued from the n-body integration was used to model the radial velocity measurements. We used the n-body code REBOUND (Rein & Liu 2012) with the WHFast integrator (Rein & Tamayo 2015) and an integration step of 0.01 days, which results in a maximum error of ∼20 ppm for the photometric model. The light-time effect (Irwin 1952) is included, although with an amplitude of ∼0.05 s in the transit timing; this is a negligible effect for this system. The model was parametrised using the stellar mass and radius, planet-to-star mass ratios, planet b-to-star radius ratio, and Jacobi orbital elements (Table 4)  180 • , and we limited the inclination of the outer one i c > 90 • . We used a GP regression model, with an approximate Matern kernel (celerite, Foreman-Mackey et al. 2017) for the model of the error terms of the transit light curves. We used different kernel hyperparameters for each transit, except for the TESS data for which we used different kernel hyperparameters for each sector. We added one dilution factor 4 for each TESS sector, and another one for the WASP transits. For each photometric data set, we added a transit normalisation factor and an additive jitter parameter. For the radial velocity, we added a systemic radial velocity and an additive jitter parameter. In total, the model has 90 free parameters. We used normal priors for the stellar mass and radius from Section 3, a non-informative sinusoidal prior for the orbital inclinations, and non-informative uniform prior distributions for the rest of the parameters. The joint posterior distribution was sampled using the emcee algorithm (Goodman & Weare 2010;Foreman-Mackey et al. 2013) with 1 000 walkers with starting points based on the results of Hébrard et al. (2020) and Section 5.

Results
In Table 4 we list the prior, the median, and the 68% credible interval of the inferred system parameters' marginal distributions. The one-and two-dimensional projections of the posterior sample are shown in Figure Figure 6 shows the posterior of the planet orbits. With an inferred impact parameter of 21.2 ± 4.5, the transits of planet c are highly disfavoured, which is in agreement with the null result of the transit search in Maciejewski et al. (2020). 4 The definition of the dilution factor and transit normalisation factor is different from the one in juliet. For the photodynamical model, we  and planet c (orange). The origin is the system barycentre, and the orbits are projected in the sky plane seen by the observer (X-Y, left) and X-Z plane (right, system top view, the movement is clockwise, the positive Z-axis points towards the observer). A thousand random orbits were drawn from the posterior samples, and the MAP is shown as a black orbit. The black points mark the position of the star (size to scale) and the planets at the central time of the RISE transit for the MAP (the size of planet b is enlarged by a factor of 10, and the size of planet c is not known). Bottom: Idem. as for 1000 stable solutions (Section 7.5).
The derived parameters agree with the ones reported by Hébrard et al. (2020), but they have a significantly improved precision. A difference, however, is the slightly larger radius ratio R p /R . Instead of M p sin i, the photodynamical modelling allows the true mass of planet c to be measured at M c = 0.424 ± 0.046 M J , which is in agreement with the dynamical upper limit of 0.60 M J reported by Hébrard et al. (2020). Most of the posterior samples have the apses of the orbits that librate around alignment. From tests we determined that the improved precision in model posteriors over the earlier study of Hébrard et al. (2020) is due to both the additional data and the photodynamical modelling that account for the interactions between the planets. The posterior transit times, which rely on the three-body system hypothesis, agree with the individually derived transit times (Section 5). In addition, the periodicity of the TTVs agrees with the ≈450 days super-period (Lithwick et al. 2012) for two planets near the 4:1 mean-motion resonance 5 . The periodicity of the TTVs also agrees with the one reported by Hébrard et al. (2020).

Transit-timing variations
The posterior timing from the photodynamical modelling is quite wide at some epochs (lower panel of Figure 7 and A.4). The uncertainty in the posterior transit times is related with the  We note that P and T 0 should not be confused with the period or the time of conjunction, respectively, and they were only used to reduce the correlations between jump parameters, replacing the semi-major axis and the mean anomaly at t ref .  (Hébrard et al. 2020 in red, and Section 5 in black). The thick grey line represents a linear ephemeris computed using only the transits observed by TESS (whose residuals are shown in the small panel in the upper right). In the lower panel, the posterior median transit-timing value was subtracted to visualise the uncertainty of the distribution. The posterior median transit time was also subtracted from each observed epoch for the individual transit-time determinations to allow for better comparison with the photodynamical modelling. The orange curve in the upper panel represents the variation in the times of inferior conjunction for planet c. knowledge of the system parameters. Future transit observations should favour the epochs where the posterior transit uncertainty is large to further improve the characterisation of the system. There is room for improvement with ∼1 minute transit-timing precision observations. Predictions of transit times up to 2026 are listed in Table A.5.
If only the TESS observations are considered, TTVs would not have been detected (Figure 7). TESS nearly continuous photometry observations of WASP-148, with a time span of approximately nine planet b periods, are insufficient to detect the TTVs of planet b. This is due to the particular configuration of the TESS observations, which accidentally only cover a part of increasing TTVs (see upper panel Figure 7). Thus, it is possible that a similar situation is occurring in other systems on which TTVs are also missed. The period derived using only the TESS data lasts 8.80604 ± 0.00014 days, and from all the observed transit times it is 8.8038083 ± 0.0000026 days (Section 7.2). This is a difference of 193 ± 12 seconds. Using the TESS ephemeris to predict future transit timing induces an offset of ∼2.2 hours/year.
The posterior TTVs' evolution in the near future is shown in Figure A

Model-independent linear ephemeris
The predictions issued from the photodynamical model and presented in Table A.5 have the disadvantage of having been obtained under the assumption that the two-planet model is correct. The results from the transit-only analysis presented in Section 5 can be used to provide a model independent ephemeris which should provide valid predictions, albeit less precise ones, even if the system is discovered to contain additional planets in the future.
A value of a mean period and time of transit can be straightforwardly produced by fitting a slope to the transit times in Table A.2. However, in the presence of transit-timing variations, a slope is not a flexible enough model to describe the transit times. Therefore, the parameters inferred from such a model are likely biased (e.g. Bishop 2007). Hébrard et al. (2020) dealt with this by inflating the error bars of the individual transit times to reach χ 2 ∼ 1. Here we decided to use a non-parametric model to describe the variation of the transit times over the linear ephemeris model.
More precisely, we chose a GP regression model whose mean function was specified to be a linear function with parameters -slope and intercept-to be inferred from the data (see Rasmussen & Williams 2005, Section 2.7). If normal priors are chosen for these parameters, the computation of the marginal likelihood can be performed analytically. We chose diffuse priors with widths of 100 days, centred at 8.804 days for the slope (period) and at the observed time of the second TESS transit observed in sector 25 for the intercept.
We tried several kernel functions to define the covariance function of the GP. We used a modified version of the implementation in scikit-learn (Pedregosa et al. 2011) and optimised the hyperparameters using the L-BFGS-B (Byrd et al. 1995;Zhu et al. 1997) and the Sequential Least SQuares Programming (SLSQP) algorithms (Kraft et al. 1988), implemented in the scipy package (Virtanen et al. 2020). All tested kernels produced almost identical results for the parameters of the linear ephemeris. We discuss below the results issued from the kernel choice that provided the largest value of optimised marginal likelihood.
The kernel function, k(x, x ), producing the largest marginal likelihood value was a exp-sine-squared (ESS) kernel, without exponential decay: where P and are the kernel hyperparameters that correspond to the period and length-scale, respectively, and A ESS is the amplitude of the covariance function. The posterior mean and 68.3% intervals for the TTVs are shown in Fig. 8. A kernel function with an additional decay term produced similar results, with a decay timescale far exceeding the time span of the observations. Under the assumption of normal priors, the linear model parameter posterior is also normal. We find the posterior ephemeris: from this model have a covariance which includes additional terms coming from the interaction between the mean function and the non-parametric part of the model. In fact, the predictions of the transit times remain precise to better than 10 minutes for over 2640 transits of planet b, that is to say over 60 years. Comparison with the transit times reported in Table A.5 shows that this method predicts transit times that are in agreement with the fully photodynamical one to better than 5 minutes up to the end of 2026.
In comparison, an ordinary least squares (OLS) fit to the transit times, with uncertainties scaled to have a reduced χ 2 of one as presented in Hébrard et al. (2020), produces a period of P = 8.803824 ± 0.000029 days, which is in agreement with the result from above and with the values presented in Hébrard et al. (2020). The OLS intercept estimator has a standard deviation six times smaller than the one in our model, which is probably unrealistic. However, the predictions from the OLS model remain precise to better than 10 minutes for less than four years.
The optimised covariance amplitude is A ESS = 25.3 minutes, and the period P is 49.7 orbits, corresponding to 437 days. The length-scale was fixed to 2.0. The fact that the model with the largest marginal likelihood does not include a long-term evolution of the transit times means that with the current data, such a trend is not detected.

Mutual inclination
We inferred a mutual inclination between the planets 7 (Figure 9). With these values, a coplanar system, as assumed in Maciejewski et al. (2020), was discarded. Maciejewski et al. (2020) tried a noncoplanar model and found a best fit solution (with Ω c − Ω b ≈ −17 • and i c = 47 • , or Ω c − Ω b ≈ +17 • and i c = 133 • ) which corresponds to a mutual inclination of ∼46 • , although they also found that the Bayesian information criterion (BIC) disfavours the non-coplanar solution. However, for high-dimensional models such as these ones, the BIC is known to provide unreliable results (Díaz et al. 2016;Nelson et al. 2020). The coplanar solution in Maciejewski et al. (2020) has a lower eccentricity for planet b (Figure 10). From a stability analysis of the orbital solution they derived, Hébrard et al. (2020) found an upper limit of 35 • for the mutual inclination. We note that the reduced-χ 2 level 7 cos I = cos i b cos i c + sin i b sin i c cos (Ω b − Ω c ) curves of the Newtonian fit plotted in Fig. 9 of Hébrard et al. (2020) favoured mutual inclinations around 30 • .
We tried to investigate which observable favours the significant mutual inclination. For this, we repeated the photodynamical analysis (Section 6) assuming coplanar orbits; we fixed the longitude of the ascending node of both planets to the same value, and we matched the orbital inclination of planet c to the one of planet b, with the latter still being a free parameter. The results are presented in Table A .6,and Figures 10,A.4,and A.5. The TTVs' posteriors of coplanar and inclined orbits are roughly the same for the observed transits ( Figure A.4). On the other hand, the transit duration posteriors are different, the transit duration of the observed transits remains almost constant for the coplanar model, whereas it decreased for the model with inclined orbits (Figure A.5). However, we found that the precision of individually determined transit durations 8 is not enough to confirm the results of the modelling with inclined orbits independently. The photodynamical modelling is in principle more sensible than the analysis on individual transits (Almenara et al. 2015), but the heterogeneity of the transit observations analysed in this work call for caution. Small variations in transit shape or timing intrinsic to the different observations could be wrongly interpreted by photodynamical modelling as an evolution in the orbital parameters 9 . If not correctly taken into account, the limb darkening dependence with the observation bandpass could be misinterpreted as changes in the impact parameter, and therefore high mutual inclinations between the planets. However, no correlation is seen between the mutual inclination and the limbdarkening parameters. In addition, the posterior distributions of the limb-darkening parameters agree with those expected from theoretical computation (Claret & Bloemen 2011;Claret 2017), although they are much wider.
In addition, models starting with a co-planar configuration were run. We found that the co-planar region of parameter space is left quickly by the MCMC walkers. The detection of a nonnegative mutual inclination does not seem to be produced by an inadequate exploration of a space parameter, or by bad mixing or lack of convergence of the MCMCs. Future observations focussed on distinguishing between the increasingly and differently predicted transit duration could conclude about the mutual inclination.

Photodynamical modelling without stellar priors
The photodynamical modelling of photometry and radial velocity data allows one to measure the absolute radius and mass in multi-planetary systems (Agol et al. 2005;Almenara et al. 2015). To test the precision of this determination in this system, we ran the same analysis as in Section 6, but with non-informative priors for the stellar mass and radius. Masses, radii, densities, and orbital parameters are listed in Table A 32% relative uncertainty for the star and planet b. The bulk densities of the star and planet b were determined with a precision of 20% and 39%, respectively. The masses of the star, planet b, and planet c were determined with a precision of 89%, 58%, and 82%, respectively. The precision on planet masses outperform the one on the star (Almenara et al. 2018b). The posterior eccentricities are also more precise for the photodynamical modelling with the stellar priors ( Figure 10). This means that the photoeccentric effect (Dawson & Johnson 2012) puts constraints on the eccentricities in addition to the ones coming from the TTVs. The same is true for the mutual inclination.

Stability analysis
The dynamical analysis of the WASP-148 orbital solution reported in Hébrard et al. (2020) has shown that the system is stable, despite significant mutual gravitational interactions between the planets. For this study, we repeated a similar stability study, but instead we performed a global frequency analysis (Laskar 1990(Laskar , 1993  (black) and 0.9 (red) for the decimal logarithm of the stability index D used in . The red zones correspond to highly unstable orbits, while the dark blue region can be assumed to be stable on a billion-year timescale.
D was derived with the frequency analysis of the mean longitude, that is 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 significantly vary for chaotic trajectories. In Figure 11 we show the distribution of the 50,000 samples projected in a (Ω c , i c ) diagram, which corresponds to the two less constrained parameters. The colour index gives the value of D for each solution. The values of log 10 D < −6 for both WASP-148 b and c correspond to stable systems on scales of billions of years . Only a small region from this diagram is stable for Ω c < 205 • and i c < 120 • , corresponding to 1239 solutions (2.5% of the total). In Figure 9 we show the probability density function (PDF) of this stable subset of solutions. The solutions cluster around a mutual inclination of 20.8 ± 4.6 • . We hence conclude from that three-body analysis that WASP-148 can only be stable for mutual inclinations below about 30 • . In the last two columns in Table 4, we provide the system parameters corresponding to the MAP, the median, and the 68% credible interval of the stable solutions. As in Hébrard et al. (2020), we also explored the stability around the stable MAP solution (Table 4), which we refer to as the nominal solution henceforward. As expected, in the (a c , e c ) domain, we confirm that it lies in a stable area, with the orbits close to the 4:1 mean-motion resonance. In the (i c , Ω c ) domain (Figure 12), we confirm that stable orbits must have a mutual inclination of I 30 • . In comparison to the stability map shown in Fig. 9 in Hébrard et al. (2020) with i c < 90 • , here we only zoom into the most stable regions, and imposed i c > 90 • : this is equivalent from a dynamical point of view, but now the stable solutions are centred around Ω c = 180 • rather than Ω c = 0 • . The main difference in the new analysis is the fact that the nodes and the inclination are now constrained by the observations. The two stability analyses presented above use different approaches: Fig. 11 presents a stability analysis carried out on different solutions allowing all parameters to vary, whereas ig. 12. Global stability analysis of the WASP-148 planetary system. We fixed all orbital parameters of the stable MAP solution (Table 4), and we varied the inclination, i c , and the longitude of the ascending node, Ω c , of planet-c. The step size is 0.25 • in both axes. 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. White dashed curves give the isolines of constant mutual inclination I = 10 • , 20 • , 30 • , and 40 • . The white dot marks the position of the stable MAP solution from Table 4. The colour bar corresponds to the one in Figure 11. the other parameters are fixed at their MAP values. We show here that they provide similar results in terms of the derived parameters for stability, which supports the reliability of both approaches and their results.

Orbital evolution
To explore the dynamics of the system, we analyse 1000 stable samples from the joint parameter posterior distribution of the photodynamical model and performed numerical integrations 10 for 1 kyr after t ref .
The results for the selected parameters are 10 With the same n-body integrator and time step used in Section 6. plotted in Figure A.6 and A.7. The posterior median and 68.3% credible interval of the mutual inclination over 1 kyr is 20.3 ± 4.5 • , and for the eccentricities it is e b = 0.126 ± 0.076 and e c = 0.217 +0.018 −0.027 . The mutual inclination remains above ∼10 • at the 95.4% credible interval. Over the 1 kyr integration, the orbital inclination of planet c is too low for transits to occur for most of the samples. Interestingly, planet b only transits for a small fraction of that time (Figure A.6); in particular, it will not transit anymore after about 200 years, then it will transit again in about 600 years.

Tidal evolution
The semi-major axis of the innermost planet is only 0.082 au, which means that the planet is close enough to the star to undergo some tidal evolution. As a result, the eccentricity can be damped and the inner orbit circularised (e.g. Hut 1981). Adopting a value identical to Jupiter's value (Lainey et al. 2009) for the tidal quality factor Q = 10 5 , we get a characteristic timescale ∼10 Gyr (e.g.  for the circularisation, which is comparable to the lifetime of the system. Therefore, it is not surprising that at present the innermost planet still shows a significant eccentricity. Moreover, the orbits of the two planets strongly interact, and secular or resonant effects can also excite the eccentricity of the innermost planet (e.g. Correia et al. 2012Correia et al. , 2013. To check this scenario, using a direct three-body model with linear tides (Correia 2018), and general relativity corrections, we ran a simulation over 600 Myr, starting with the initial conditions of the nominal solution from Table 4. We adopted a Love number of k 2 = 0.5 and a time lag of ∆t = 1 s (equivalent to Q ∼ 10 5 ). We observed that the eccentricity of the innermost planet undergoes large oscillations owing to secular interactions. Moreover, after 130 Myr, the system crosses a resonance which pumps the inner planet eccentricity and damps the mutual inclination ( Figure 13). As a result, at the end of the simulation, the average inner planet eccentricity is higher than the initial one. We conclude that the presently observed non-zero value is compatible with the tidal evolution of the WASP-148 system. As expected, the semi-major axis of the innermost planet also slightly decreases, and we have a b ≈ 0.073 au after 600 Myr. Since the age of the star at present is already 4 Gyr, we can assume that the initial semi-major axis was slightly larger than the present value, and so the two planets could even be trapped in the 4:1 mean motion resonance. This is an interesting formation scenario that deserves more attention in future work on the system.

Three-planet model
All these results rely on the three-body system hypothesis. If the dynamics of planet b is affected by additional planets other than planet c, increasing the mutual inclination could provide the additional variability in the model required to fit the data. The photodynamical model is flexible and could overfit the data. The argument of an inaccurately specified model is discussed by Petit et al. (2020) in the context of the determination of the eccentricity of the K2-19 planets.
To test the influence in the system parameters of an unaccounted for additional planet in the system, we repeated the photodynamical analysis (Section 6) with a third planet, which we started at the period of 150 days, which is the peak in the residuals of a 2 Keplerian fit of the radial velocities (Sect. 4.1 and Fig. 3 in Hébrard et al. 2020). The results are presented in Table A.8 and Figures 9, 10 Table 4, we show the evolution of the mutual inclination (top), the eccentricities (middle), and semi-major axes (bottom) over 600 Myr. We used a direct three-body model with linear tides (Correia 2018) with general relativity corrections, and adopted k 2 = 0.5 and ∆t = 1 s for both planets.
a period of 151.2 +2.1 −1.7 days and a mass of 0.262 +0.10 −0.076 M J , which is in agreement with the estimations reported in Hébrard et al. (2020). Its mutual inclination relative to the planets b and c is not well-constrained. The mutual inclination between planet b and c is reduced from 41.0 +6.2 transit timing and duration variations do not allow one to distinguish between the two-and three-planet models with the current data (Figures A.4 and A.5). More data are needed to assert the presence of additional planets in the system.

Discussion
The WASP-148 system is composed of a G5V star orbited by a hot Saturn (0.287 +0.022 −0.076 . Also, assuming only two planets in the system, their orbits have a mutual inclination of 20.8 ± 4.6 • . Both planets are located in a scarcely populated part of the mass-separation diagram (Figure 14), which cannot be fully explained by a detection bias 12 , thus indicating a low probability outcome for the planetary formation process. The scarcity of planets such as WASP-148 b can be partially explained by the low efficiency of the Lidov-Kozai mechanism (Lidov 1962;Kozai 1962) to form hot Saturns because most migrating planets are tidally disrupted (Anderson et al. 2016). WASP-148 c was placed in the 'period valley' (Udry et al. 2003) between the hot and warm Jupiters clumps.
The WASP-148 giant planets are expected to form in a protoplanetary disk beyond the snow line, which is located at a few au (Lecar et al. 2006), and then migrate inwards (Goldreich & Tremaine 1980). Convergent migration leads to planets being 11 The period ratio at t ref is 3.92411 +0.00014 −0.00017 and computed as (a c /a b ) 3/2 , with a c /a b = 2.487868 +0.000061 −0.000073 . For the stable MAP solution, the period ratio (computed over a 1000 orbits of the planet c) is 3.922. 12 While the detection of systems similar to WASP-148 is accessible for high-precision stable radial velocities' instruments, they are difficult to detect with transit searches. Due to the mutually inclined orbits and the secular oscillation of the orbital inclination, the probability that both planets transit for a given observer is low (Section 7.6).
Article number, page 12 of 26 J.M.  captured in mean-motion resonances, but usually of first order. High order resonances such as the 4:1 one require high initial orbital eccentricities (Rein & Papaloizou 2010). However, disk migration does not favour significant orbital eccentricity excitation because planet-disk interactions tend to damp eccentricities (Bitsch et al. 2013;Dunhill et al. 2013). After the disk disappears, planet-planet scattering can cause high eccentricity migration for the inner planet (Rasio & Ford 1996), or the eccentricities can be excited by chaotic secular interactions (Wu & Lithwick 2011). Both mechanisms are compatible with the noncoplanarity of the observed orbits.
Another possibility that does not require planet-planet scattering or secular chaos is explored in Lee & Peale (2002): the planets are captured in the mean-motion resonance and are massive enough to open gaps in the disk. With the mechanism for damping the eccentricity being reduced, the orbits can become elliptical during the inward migration within the resonance (Artymowicz 1992). Capture in the 4:1 mean-motion resonance is then able to produce some inclination excitation given that the inner planet is not too massive (Thommes & Lissauer 2003;Libert & Tsiganis 2009).
Finally, at a later stage, tidal interactions with the star should shrink the inner planet orbit afterwards and pull the system out of resonance. Tidal interactions usually also damp the eccentricity, but the presence of the outer planet delays this process due to secular and resonant interactions.
There are only a few other planetary systems for which a large mutual inclination has been reported: 102.3 +7.4 , and some systems with ultra short period planets with mutual inclinations larger than 10 • reported in Dai et al. (2018). HD 3167 is a system with three small-sized planets, including two transiting that allowed Rossiter-McLaughlin measurements. π Mensae is a hierarchical system with a period ratio of ∼330, and a large difference in planet masses. υ Andromedae is composed of three planets, and the mutual inclination is measured for the outer pair of super Jovian planets with a period ratio of 5.3. Kepler-108, composed of two Saturn-mass planets with a period ratio of 3.88, is the one that resembles WASP-148 the most, but with longer period planets in lower eccentricity orbits. Interestingly, the period ratio is also close to the 4:1 mean-motion resonance.
As shown in Section 7.8, additional planets not accounted for in the modelling can affect the determination of the system parameters, in particular the mutual inclination. Thus more data are needed to conclude on the true architecture of the system.
Acknowledgements. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https: //www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. Part of these simulations have been run on the Lesta cluster kindly provided by the Observatoire de Genève. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. AC acknowledges support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal. GM acknowledges the financial support from the National Science Centre, Poland through grant no. 2016/23/B/ST9/00579. MF acknowledges financial support from grant PID2019-109522GB-C5X/AEI/10.13039/501100011033 of the Spanish Ministry of Science and Innovation (MICINN). MF, VC and JS acknowledge financial support from the State Agency for Research of the Spanish MCIU through the Center of Excellence Severo Ochoa award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709     Notes. The table lists  Notes. N(µ, σ): Normal distribution with mean µ and standard deviation σ. T N(µ, σ, a, b): Normal distribution with mean µ and standard deviation σ, truncated between a lower a and upper b limit. U(a, b): A uniform distribution defined between a lower a and upper b limit. J(a, b): Jeffreys (or log-uniform) distribution defined between a lower a and upper b limit.
Article number, page 17 of 26 A&A proofs: manuscript no. output