The CARMENES search for exoplanets around M dwarfs. Dynamical characterization of the multiple planet system GJ 1148 and prospects of habitable exomoons around GJ 1148 b

Context. GJ 1148 is an M-dwarf star hosting a planetary system composed of two Saturn-mass planets in eccentric orbits with periods of 41.38 and 532.02 days. Aims. We reanalyze the orbital configuration and dynamics of the GJ 1148 multi-planetary system based on new precise radial velocity (RV) measurements taken with CARMENES. Methods. We combined new and archival precise Doppler measurements from CARMENES with those available from HIRES for GJ 1148 and modeled these data with a self-consistent dynamical model. We studied the orbital dynamics of the system using the secular theory and direct N-body integrations. The prospects of potentially habitable moons around GJ 1148 b were examined. Results. The refined dynamical analyses show that the GJ 1148 system is long-term stable in a large phase-space of orbital parameters with an orbital configuration suggesting apsidal alignment, but not in any particular high-order mean-motion resonant commensurability. GJ 1148 b orbits inside the optimistic habitable zone (HZ). We find only a narrow stability region around the planet where exomoons can exist. However, in this stable region exomoons exhibit quick orbital decay due to tidal interaction with the planet. Conclusions. The GJ 1148 planetary system is a very rare M-dwarf planetary system consisting of a pair of gas giants, the inner of which resides in the HZ. We conclude that habitable exomoons around GJ 1148 b are very unlikely to exist.


Introduction
In the past twenty years, M-dwarf stars have been the primary targets for a number of planet search surveys via the precise Doppler method (Marcy et al. 1998(Marcy et al. , 2001Delfosse et al. 1998;Endl et al. 2003;Kürster et al. 2003;Bonfils et al. 2005;Butler et al. 2006;Reiners et al. 2018a). These efforts are motivated by the fact that M dwarfs are the A&A 638, A16 (2020) 0.354 ± 0.015 Schw19 P rot [d] 71.5 ± 5.1 DA18 References. 2MASS: Skrutskie et al. (2006) predominant population of stars in the solar neighborhood. Their significantly lower mass and luminosity when compared to Sunlike stars make them suitable for the detection of lower mass planets orbiting in the so-called habitable zone (HZ), where the surface temperature would be favorable for liquid water to exist. Planets in the HZ have already been discovered around low-mass M dwarf neighbors such as Proxima Centauri (Anglada-Escudé et al. 2016), LHS 1140 (Dittmann et al. 2017), HD 147379 (Reiners et al. 2018b), GJ 752A (Kaminski et al. 2018), and GJ 357 (Luque et al. 2019), among others. Of special interest for Doppler surveys are also potentially habitable M-dwarf multiple planet systems consisting of Earth-size planets such as GJ 1069 (Dreizler et al. 2020), or the ultra-cool M dwarfs TRAPPIST-1 (Gillon et al. 2017) and Teegarden's star .
As of October 2019, more than 100 M-dwarf planetary systems have been discovered via the radial velocity (RV) method (Martínez-Rodríguez et al. 2019) 1 . Many more are expected to be detected with the transit technique thanks to the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), which already delivered the first planet detections around M-dwarf stars such as L 95-98 (Kostov et al. 2019;Cloutier et al. 2019), TOI 270 (LEHPM 3808; Günther et al. 2019), and GJ 357 (Luque et al. 2019), which are in fact all multi-planet systems.
The recent bulk of multi-planetary discoveries around M dwarfs opens a great opportunity to study their orbital architecture in detail. Understanding the dynamics of multi-planet systems around M dwarfs provides an important clue to planet formation around these low-mass stars (Lissauer 2007;Raymond et al. 2007;Zhu et al. 2012;Anglada-Escudé et al. 2013;Coleman et al. 2017;Morales et al. 2019). Furthermore, the planetary dynamics affect the evolution of any natural satellites; these exomoons can provide critical information on the formation and evolution of a planetary system (see Heller et al. 1 See also http://exoplanet.eu 2014), and may even be habitable themselves (Williams et al. 1997;Heller & Barnes 2013;Forgan & Dobos 2016). While no exomoon detection has been confirmed (see, e.g., Teachey & Kipping 2018;Heller et al. 2019), theoretical studies have shown that the planets' orbital oscillations can drive the orbital evolution of satellites (Zollinger et al. 2017) and that the tidal evolution of the satellite's orbit about the host planet can also be significant (e.g., Martínez-Rodríguez et al. 2019).
In this paper, we provide the first detailed analysis of the dynamics of the rare GJ 1148 system, which was reported to host two Saturn-mass planets orbiting around the M dwarf on moderately eccentric orbits (Trifonov et al. 2018a). We provide new precise Doppler measurements from the CARMENES 2 instrument (Quirrenbach et al. 2016;Reiners et al. 2018a) and focus on the long-term dynamical properties of the system, which could provide important information on the primordial protoplanetary conditions needed to form two Saturn-mass planets around a low-mass star.
The paper is organized as follows. In Sect. 2 we introduce the updated physical properties of the M-dwarf star GJ 1148 and provide a short summary of the literature on the known planetary system. In Sect. 3 we present the observational data, and in Sect. 4 we present our data analysis. In Sect. 5 we detail our dynamical results of the GJ 1148 system. In Sect. 6 we provide our conclusions and a brief discussion.

The GJ 1148 planetary system
The planet GJ 1148 b was discovered by Haghighipour et al. (2010) based on 37 velocities taken with the Keck/HIRES spectrograph (Vogt et al. 1994). It is an eccentric Saturn with a period of ∼41.4 d and, interestingly, mostly resides in the optimistic habitable zone (HZ). After the public release of the HIRES velocity archive by Butler et al. (2017), it became clear that the 2 Calar Alto high-Resolution search for M dwarfs with Exo-earths with Near-infrared and optical Echelle Spectrographs, http://carmenes. caha.es   The photometric measurements from SuperWASP, NSVS, and HATNet are densely sampled, and thus most of the GLS power appears significant, but it cannot be associated with the RV signals. The most likely stellar rotation frequency is shown with a red dashed line, while the red shaded area denotes its uncertainties, as given by Díez Alonso et al. (2019). extended HIRES data time series yields an additional RV signal with a period of ∼530 d, likely due to a second massive planet in the system. In Trifonov et al. (2018a), we added 52 CARMENES RVs and performed a combined Doppler analysis, which revealed the two-planet configuration of the GJ 1148 system. We found the following orbital configuration for GJ 1148 b: m b sin i = 0.304 M Jup , P b = 41.380 d, e b = 0.379, a b = 0.166 au, and for GJ 1148 c: m c sin i = 0.214 M Jup , P c = 532.6 d, e c = 0.341, a c = 0.913 au. Figure 1 shows the architectures of multi-planet systems known to orbit M dwarfs, including GJ 1148. To date the observational results indicate that the planetary population around M dwarfs consists predominately of low-mass planets in the range from a few Earth masses to Neptune mass. Together with the GJ 876 system (Rivera et al. 2010;Nelson et al. 2016;Millholland et al. 2018) the GJ 1148 system belongs to the rare population of M dwarfs known to host a pair of Jovian planets. In contrast, observations of solar-mass main-sequence stars and more massive giant and sub-giant stars show a higher frequency of Jupiter-mass planets, which seems to be correlated with stellar host mass (Fischer & Valenti 2005;Reffert et al. 2015). These findings suggest that more massive stars tend to have more massive planets, likely due to a more massive primordial protoplanetary disk from which planets form. In this context M dwarfs are not expected to have a large population of Jupiter-mass planets, which makes systems like GJ 1148 an important discovery that may provide insights into the primordial disk properties around low-mass stars needed to accumulate gas giant planets.
Apart of the GJ 1148 and the GJ 876 M-dwarf planetary systems, there are some other peculiar exceptions; for example, GJ 1046 hosts a brown dwarf companion on 169-day period orbit (Kürster et al. 2008;Trifonov et al. 2018b), and there is an unexpected giant planet around the very low-mass M-dwarf star GJ 1132 (Morales et al. 2019). The GJ 317 system also has a Jupiter-mass planet in orbit (Johnson et al. 2007 Trifonov et al. (2018a), there are 24 additional CARMENES RVs, and the HIRES data were corrected for small systematics. The data uncertainties include the estimated RV jitter.
The observational evidence suggest that the GJ 1132 and GJ 317 M dwarfs are very likely accompanied by another long-period giant planet whose orbit has not yet been precisely determined. In our first paper on GJ 1148 (Trifonov et al. 2018a), we performed dynamical modeling of the GJ 1148 RV data and we studied the long-term stability of the best fit to the system, but at that time we did not perform a detailed long-term dynamical analysis of the GJ 1148 system. In this second paper, we aim to analyze the possible orbital configurations and planetary mass regimes based on an extensive long-term dynamical analysis that includes additional CARMENES Doppler data from recent observations.
Based on analyses of CARMENES spectra taken between January 2016 and January 2019, Schweitzer et al. (2019) and Passegger et al. (2018) updated the physical parameters of the GJ 1148 host star. The results from these two studies report that GJ 1148 is a typical M4.0V star with a stellar mass of 0.354 ± 0.015 M , effective temperature T eff = 3358 ± 51 K, and luminosity L = 0.0143 ± 0.0003 L .

HIRES data
HIRES is a high-resolution spectrograph mounted on the 10 m Keck-I telescope, Keck Observatory, Hawaii, USA, (Vogt et al. 1994); via an iodine (I 2 ) cell (Marcy & Butler 1992) it can achieve RV precision down to ∼3 ms −1 (Butler et al. 1996). Butler et al. (2017) have published 125 precise HIRES RV measurements of GJ 1148 as part of a large archive of precise Doppler measurements and stellar-line activity-index measurements (∼65 000 HIRES spectra of ∼1700 stars), obtained between 1996 and 2014. Recently, the large sample size of the HIRES RV archive has allowed Tal-Or et al. (2019) to identify small but significant systematic nightly zero-point variations in the data on the order of ∼1 ms −1 , and to correct for them, making the HIRES data even more precise.  Fig. 3. Due to its lower precision the near-IR data was not used in the data analysis, but is only overplotted with an optimized RV offset to the best two-planet dynamical model constructed from the optical data shown in Fig. 3. The CARMENES near-IR RVs are in an excellent agreement with the model from the available optical data, fully recovering the two planetary signals.
We note that for GJ 1148 in particular, the small systematic corrections found in Tal (Butler et al. 2006), and thus, in this paper we decided to perform our analysis with the data provided by Tal-Or et al. (2019).

CARMENES data
A detailed description of the CARMENES instrument can be found in Quirrenbach et al. (2016), while the CARMENES survey and goals are described in more detail in Reiners et al. (2018a). The standard CARMENES data reduction flow is described in Zechmeister et al. (2018).
Briefly, the CARMENES instrument consists of two highresolution spectrographs, which are designed to perform stateof-the-art RV measurements in the optical (0.52-0.96 µm) and in the near-IR (CARMENES-NIR; 0.96-1.71 µm). We obtained 76 optical and 68 near-IR CARMENES RV measurements of GJ 1148 between January 2016 and May 2019. The uneven number of optical and near-IR CARMENES data is due to NIR channel calibration issues (typically during the first six-months of the CARMENES operation). For eight of the near-IR spectra we were unable to obtain meaningful RVs measurements, and thus these data were discarded. The first 52 optical RV measurements were published in Trifonov et al. (2018a). We computed all RVs with the latest pipeline and the SpEctrum Radial Velocity AnaLyser (SERVAL; Zechmeister et al. 2018) version, thus creating a revised and expanded data set. Therefore, we present 24 new and 52 reprocessed optical RVs and a new set of near-IR data of GJ 1148.
Additionally, using SERVAL we recomputed stellar activity index time series derived from the available CARMENES spectra: the chromatic index (CRX), the differential line width (dLW), calcium infra-red triplet (cairt), and the Hα measurements. For a detailed description of the SERVAL activity indicators we refer to Zechmeister et al. (2018). All CARMENES Doppler measurements, activity index data, and their individual formal uncertainties used for our analysis are available in Tables A.1-A.2. As in the HIRES data, the new CARMENES data are consistent with the two planetary signals with periods of ∼41.3 d and 532.6 d, respectively. There are no significant GLS peaks in the residuals of the simultaneous two-planet Keplerian fit applied to the HIRES and CARMENES optical data. In Fig. 2  The HIRES S-and H-index time series of Butler et al. (2017) lack any significant periodicity that could be associated with the two strong planetary signals seen in the RVs (Trifonov et al. 2018a). From Fig. 2, it is evident that the spectral activity indicator time series from CARMENES also do not yield significant peaks at the known planetary periods (blue dashed lines).

Periodogram analysis
Photometric observations of GJ 1148 were performed by Haghighipour et al. (2010), who found a periodic signal at P rot = 98.1 d, suggesting that this is the likely rotational period of the star (see Fig. 5 and Sect. 5 in their paper for details). However, from the Hungarian-made Automated Telescope Network (HAT-Net; Bakos et al. 2004)  In the three bottom right panels of Fig. 2, we show the GLS periodograms of the SuperWASP, the NSVS light-curve data of GJ 1148 compiled by Díez Alonso et al. (2019), and the HATNet photometry. These photometric measurements were taken with a higher cadence when compared to the RV data, and thus many frequencies in the GLS power spectrum appear significant. A well-defined significant peak is only detected in the HATNet time series, whose maximum power is within the P rot = 71.5 ± 5.1 d uncertainties (red shaded frequency range in Fig. 5. Orbital evolution of the GJ 1148 system for a Myr long N-body integration using the Wisdom-Holman scheme. From top to bottom: evolution of the planetary semi-major axes, eccentricities, and the apsidal alignment angle ∆ω = b − c . The gray area is composed of 1000 randomly chosen stable configurations from the MCMC run, which represent the possible orbital outcome consistent with the RV data. The solid color curves show the evolution of the best-fit configuration.  However, there is no definitive agreement between the SuperWASP, NSVS, and the HATNet photometry on the one hand, and the spectral Doppler measurements and activity indicators on the other. We concluded that there is no photometric periodic signal that could be firmly associated with the GJ 1148 b and c planets.

Orbital update of the GJ 1148 system
For the analysis of the precise Doppler data of GJ 1148 we use The Exo-Striker fitting toolbox 4 (Trifonov 2019). The Exo-Striker provides a large variety of fitting and sampling schemes coupled either with a standard Keplerian model or with a more complex N-body model. The latter takes into account the gravitational perturbations that occur in multi-planet systems while it models the radial component of the stellar velocity caused by the companions. For a widely separated system such as GJ 1148 there is little or no practical incentive in using the more expensive dynamical model over the simple Keplerian model. The Keplerian and N-body models are indistinguishable within the temporal baseline of the combined HIRES and CARMENES observations (Trifonov et al. 2018a). The dynamical model, however, has the advantage of being able to fit for mutually inclined orbits which distribution and long-period stability we aim to study. Therefore, our natural choice is to use the N-body RV model for the orbital analysis of the GJ 1148 system.
To derive the best fit we adopt a maximum likelihood estimator (MLE) scheme, which optimizes the N-body model's likelihood value (− ln L) via the Nelder-Mead algorithm (also known as Simplex, Nelder & Mead 1965). As in Trifonov et al. (2018a) we optimize the osculating semi-amplitude K, orbital period P, eccentricity e, argument of periastron ω, and mean anomaly M for the first HIRES observational epoch and the HIRES and CARMENES data velocity offsets. In these analyses we also include the HIRES and CARMENES velocity jitter variance as additional fitting parameters (Baluev 2009). Additionally, when we study mutually inclined orbits we also fit for the orbital inclinations i and the difference of the orbital lines of node ∆Ω. To estimate the parameter uncertainties of our best fits and to perform a parameter distribution analysis, we rely on a Markov chain Monte Carlo (MCMC) sampling using the emcee sampler (Foreman-Mackey et al. 2013), which is integrated within The Exo-Striker.
Our new best-fit analysis is largely consistent with our previous results for the GJ 1148 system presented in Trifonov et al. (2018a). Assuming an edge-on coplanar configuration, our new N-body model fitted to the CARMENES and HIRES data yields osculating orbital parameters for the inner planet GJ 1148 b:  Table 2 summarizes the best-fit parameters and uncertainties of our new coplanar edge-on N-body model of the GJ 1148 system. The best fit to the new CARMENES and HIRES RV time series is shown in Fig. 3, together with a phasefolded representation of the GJ 1148 b and c planetary signals. From Fig. 3 it is clear that the CARMENES and HIRES data are fully consistent with the presence of the two planetary signals.

CARMENES near-IR data
We tested the 68 CARMENES near-IR RVs for consistency with the two-planet signals. The near-IR data was not used in our analysis together with the available optical data, for two reasons. First, the precision of the CARMENES near-IR RV data of GJ 1148 is significantly lower than that achieved in the optical channel of CARMENES. This is expected for GJ 1148: for stars of spectral class M 4.0 V the spectroscopic information needed for precise RV measurements is still more abundant in the CARMENES-VIS spectra (see Fig. 6 in Reiners et al. 2018b). Thus, the weight of the near-IR data is much lower in the modeling. Second, the CARMENES optical and NIR data are taken practically simultaneously (usually with a very small deviation in their time stamps due to the slightly different photon mid-point center). This causes difficulties for the N-body modeling, which must automatically adapt the integration time steps to very small values in order to go between two neighboring CARMENES data points. This makes the modeling of the near-IR data together with the higher quality optical data impractical. Instead, we decided to test the near-IR data against the best-fit model constructed from the optical, following the methodology adopted by Trifonov et al. (2015) for the VLT-CRIRES (Kaeufl et al. 2004) near-IR data. So we applied this model to the near-IR data by optimizing only the RV offset and we studied its statistical properties. Red curves represent the evolution obtained from the Wisdom-Holman N-body integration evolution of the best fit (see Fig. 5), while red curves are the Wisdom-Holman N-body integration with GR precession included. Gray dashed and dotted curves represent the evolution of the best fit but derived from the octupole secular perturbation theory and the octupole evolution with a GR precession correction, respectively. Figure 4 shows the CARMENES near-IR RV measurements and a phase-folded representation of the GJ 1148 b and c planetary signals from the best-fit model from Table 2. A visual inspection shows that the signal of GJ 1148 b is clearly present in the CARMENES-NIR data. The GJ 1148 c planetary signal is somewhat sparsely sampled around periastron by the CARMENES-NIR data, yet there is good agreement in phase and amplitude between the data and the model.
The statistical properties of the near-IR data with respect to the planetary signals of GJ 1148 are summarized in Table 3. Assuming no planets (i.e., the Null hypothesis), the near-IR data have a large rms = 25.36 m s −1 , while subtracting only the innermost planetary signal from the data is strongly preferred yielding an rms = 9.34 m s −1 , and a p = 1.3 −38 . Finally, applying the bestfit two-planet model results in an rms = 8.06 m s −1 and p = 5.3 −40 with respect to the flat model. We note that in the two-planet model the CARMENES-NIR jitter is only 0.01 m s −1 , meaning that the internal data uncertainties (σ NIR = 7.7 m s −1 ) are adequately estimated. In terms of relative probability, R = e ∆ ln L , between a model that assumes only the presence of GJ 1148 b and the best-fit two-planet model, requires ∆ ln L > 5 to claim significance (see Baluev 2009). The best two-planet model applied to the near-IR data yields ∆ ln L ∼ 11 (see Table 2) over the one-planet model, making the two-planet model the much better choice. We therefore conclude that the CARMENES near-IR RVs fully agree with the best-fit model from the CARMENES optical data and the HIRES data. Therefore, the near-IR data strengthens the multi-planet hypothesis around GJ 1148.

Coplanar edge-on configurations
In Trifonov et al. (2018a) we show that the best-fit orbital configuration of the GJ 1148 system is stable for at least 10 Myr, but A16, page 7 of 18 A&A 638, A16 (2020) exhibits strong long-term secular dynamical interactions leading to large oscillations of the planetary eccentricities on a timescale of ∼80 000 yr. In our new dynamical analysis, we inspect the dynamics of the new best fit, and also the overall stability of the MCMC posterior parameter distribution consistent with the data.
Dynamical analyses are performed using a custom version of the Wisdom-Holman algorithm (Wisdom & Holman 1991), which works in the Jacobi coordinate system (e.g., Lee & Peale 2003). Our long-term stability integrations are performed for a maximum of 10 Myr with an integration time step equal to 0.4 d. We find this time step and maximum integration time to be sufficient for an accurate dynamical test of the system. For each MCMC integration, we automatically monitor the evolution of the planetary semi-major axes and eccentricities as a function of time to ensure that the system remains regular and well separated at any given time of the simulation. Any deviation of the planetary semi-major axes by more than 20% from their starting values, or eccentricities leading to crossing orbits, were considered unstable.
The coplanar edge-on MCMC posterior distributions around the best fit from Table 2, and their respective correlations are shown in Fig. A.1. We found that all the coplanar edge-on MCMC samples are stable for 10 Myr. This is perhaps not a surprise given the relatively large orbital separation between the two massive planets and taking into account that the best-fit model is very well constrained by the RV data. From the coplanar edge-on perspective, we cannot draw conclusions based upon long-term stability, but the overall dynamics of the MCMC posterior parameter distribution allow us to to shed more light on the possible dynamical architecture of the system. Figure 5 shows the 1 Myr orbital evolution of the best fit and 1000 randomly chosen configurations from our MCMC test. The planetary semi-major axes remain nearly constant with no major deviations from the starting values over the integration. The eccentricities, however, exhibit large variations in the range of e b = 0.01 to 0.60 and e c = 0.10 to 0.50. Interestingly, most of the MCMC configurations exhibit apsidal alignment with the angle ∆ω = b − c librating around 0 • with a typical semi-amplitude of 60 • ; very few configurations circulate between 0 • and 360 • . This dynamical behavior can be nicely explained using the octupole secular theory. For comparison, we also integrate the time evolution of the orbital elements for the planar hierarchical three-body problem, using the octupolelevel secular perturbation equations derived by Ford et al. (2000).
In Fig. 6 we show the evolution of the planetary eccentricities e b and e c and the apsidal alignment angle ∆ = b − c . We integrate the best-fit solution both with N-body and secular codes. For the secular code, we use SecuLab 5 up to octupole order, with corrected double averaging terms (Luo et al. 2016;Grishin et al. 2018b). Another octopule code without the corrected terms (Lee & Peale 2003) gives almost identical results. The linear secular theory (up to second order in eccentricities) predicts a secular period of 73 400 yr (see Eq. (7.9)-(7.28) in Murray & Dermott 1999). The secular integration yields a timescale of 75 800 yr. The differences are due to the large eccentricities of both planets.
General relativistic precession can also play a role in the long-term evolution of the system. The apsidal precession rate 5 https://github.com/eugeneg88/SecuLab Fig. 7. Trajectories in the phase-space diagrams of the planetary eccentricities e b and e c vs. the apsidal alignment angle ∆ = bc for the GJ 1148 system constructed from the octupole-level secular perturbation theory, with GR correction included, assuming the same total angular momentum. Blue crosses indicate the N-body best-fit parameters of e b , e c and bc and their uncertainties, while the red paths represent the trajectories starting at the best fit. The blue dotted lines represent the maximum planetary eccentricities, which can be attained for the total angular momentum of the system. of planet b is (Blaes et al. 2002;Liu et al. 2015) and the typical timescale is T GR = 2π/ω GR ≈ 1.54 Myr. Including general relativity (GR) precession in the secular code yields a slightly shorter secular period of 72 500 yr. The simplified expression for the secular timescale with GR is T sec,GR = T sec /(1 + T sec /T GR ) ≈ 72 200 yr, which slightly underestimates the timescale obtained from numerical integration. This is due to the large variations in the eccentricity, which affects the GR precession rate as ∝ (1 − e 2 b ) −1 . Overall, the Wisdom-Holman N-body code and the secular code are consistent with each other. Apart from the slightly different timescales, the phase and amplitude of the osculating parameters of the octupole prediction are clearly consistent with those from the N-body dynamical evolution of the GJ 1148 system, including GR corrections. Moreover, the Wisdom-Holman code with GR corrections applied on the GJ 1148 system is the first public code available which has been compared to secular theory. Only very recently have other public N-body codes started incorporating GR corrections (Tamayo et al. 2019). Figure 7 shows the dynamical behavior of the planetary eccentricities and apsidal alignment for the co-planar configuration in the phase-space diagrams of e b and e b on the y-axis and the apsidal alignment angle ∆ = bc on the x-axis. We show the trajectories calculated from the octupole-level secular perturbation theory with GR precession, assuming the same total angular momentum as for the GJ 1148 system, but with different initial e b , e c , and ∆ . The blue crosses in Fig. 7 indicate the position of the best fit and its uncertainties, while the red trajectory is the evolution starting from the best fit. Fig. 7 shows that ∆ of

Stability of the mutually inclined configurations
We now study the mutual inclination limits of the GJ 1148 system by allowing the planetary inclinations i b and i c and the difference of the line of nodes ∆Ω = Ω b -Ω c to vary together with the remaining parameters in the dynamical modeling. Initially we tried a simplex optimization with i b , i c , and ∆Ω being free, but this did not yield a significantly different best fit when compared to the co-planar edge-on model. The likely reason is that there is no sharp − ln L minimum in the i b , i c , and ∆Ω space to which the simplex algorithm could to converge. This is expected since constraining the mutual inclinations from RV data requires strong interactions between the planets over the temporal baseline of the observations. The best example is the other massive pair around an M dwarf; for the GJ 876 system, Rivera et al. (2010), Nelson et al. (2016), and Millholland et al. (2018) successfully measured the orbital inclinations from RV data. For GJ 876, the timescales of the perturbations are shorter than the observational data time span, and thus it is possible to constrain the mutual inclination. As we showed in Sect. 5.1, the timescales of the perturbations in the GJ 1148 system are much longer, and thus it is hardly possible to constrain the inclinations by modeling the RV data.
Therefore, instead of optimizing the − ln L via the simplex algorithm, we simply ran an MCMC test starting from the best coplanar edge-on fit and we allowed emcee (Foreman-Mackey et al. 2013) to explore the parameter space consistent with the combined HIRES and CARMENES RVs. For all parameters in this test we adopted flat priors with rather relaxed boundaries; the Notes. The posterior confidence intervals for the other planetary parameters are consistent with those obtained for the coplanar edge-on case (see Table 2), and are thus not shown here.
where ∆Ω = Ω b -Ω c is the difference between the orbital line of nodes. Since only ∆Ω is important for deriving ∆i, Ω b can remain fixed at 0 • . Figure 8 shows the results from our dynamical MCMC allowing for mutual inclinations. The samples that survived 10 Myr without leading to a strong chaotic behavior are found only in prograde configurations in the range 0 • < ∆i 60 • or, less likely but still possible, on a retrograde mutually inclined A16, page 9 of 18 A&A 638, A16 (2020) Fig. 9. Face-on representation of the orbital evolution of the best fit of the GJ 1148 system with a time step of 9000 yr. covering a full secular cycle of the planetary eccentricities, which has a period of ≈ 72 500 yr. Ellipses illustrate the planetary orbits with the line connecting the central star (red dot) and the planetary argument of periastron . The timescale for orbital precession is ≈ 52 000 yr, thus shorter than the secular period of the system. The gray shaded disk represents the optimistic habitable zone around the GJ 1148 M dwarf. Most of the time GJ 1148 b orbits mostly within the HZ. configuration in the range 130 • ∆i < 180 • . With these results, we cannot firmly constrain the mutual inclinations between the planets, but we can derive limits on their dynamical masses. Most of the configurations in the MCMC posterior distribution can be found in the range between 57.36 • < i b < 122.97 • , and between 58.18 • < i c < 119.64 • . The sin i factor for these inclinations is 1 ≥ sin i > ∼ 0.85, which leads to maximum dynamical planetary masses of m b = 0.432 M Jup and m c = 0.313 M Jup . Table 4 summarizes the stable posterior 1σ confidence intervals of the dynamical masses and semi-major axes constrained by the possible orbital inclinations of the GJ 1148 system.

Prospects of habitable exomoons around GJ 1148 b
The GJ 1148 system is very interesting, particularly because the inner giant GJ 1148 b resides inside the optimistic HZ. Figure 9 shows a face-on representation of the secular cycle of the GJ 1148 system. Each panel of Fig. 9 shows a time step of 9000 yr of the evolution, covering a full cycle of the planetary eccentricities which has a period of ≈72 500 yr and a full cycle of the orbital precession of the system, which has a period of ≈52 000 yr. The two gray overlapping annuli in the panels shows the optimistic and conservative HZ range around the GJ 1148 M dwarf. For the HZ estimate we made use of the Habitable Zone Calculator 6 , which relies on the work by Kopparapu et al. (2013), who provide HZ estimates around MS stars with effective temperatures in the range of 2600 K < T eff < 7200 K. From Fig. 9 it is evident that GJ 1148 b is orbiting within the optimistic HZ, but it leaves the conservative HZ temporarily when its orbit is most eccentric.
While GJ 1148 b is unlikely to be hospitable to life as we know it, the question arises of whether it could host habitable moons. We refer the reader to Zollinger et al. (2017) for a more detailed description of exomoons around Saturn-mass planets orbiting M-dwarf stars and to Heller & Barnes (2013, It is reasonable to assume that GJ 1148 b was formed farther out, beyond the protoplanetary ice-line around the GJ 1148 M dwarf where it can grow massive enough to become a Saturn-mass planet, and migrated inwards to more temperate orbits. Given our best local example, Saturn, it is also reasonable to assume that GJ 1148 b has icy exomoons, or had them in the past. In the HZ, such exomoon bodies may not be icy anymore, but potentially habitable, ocean-like worlds. The formation of moons around planets is similar to the formation of planets in protoplanetary disks. Due to type I migration torques, the regular moons, which are close to the planet, are thought to have formed in situ during the late stages in a starved circumplanetary disk before its dispersal (Canup & Ward 2002). Conversely, the irregular satellites, which usually orbit farther out, are smaller in size, and have eccentric and inclined orbits, are thought to have been dynamically captured (see Jewitt & Haghighipour 2007 for a review).
Additionally, exomoons will have some benefits with respect to planetary bodies orbiting around an M dwarf inside the HZ. For example, an M-dwarf HZ planet will be tidally locked to the star, with only one side facing the star. It is still debatable whether such tidally locked planets could host life due to the large temperature contrast between the day and night sides of the planet, but most recent research suggests tidally locked planets can maintain surface liquid water (Joshi et al. 1997;Pierrehumbert 2011;Yang et al. 2013;Del Genio et al. 2019). On the other hand, exomoons will be tidally locked to the giant planet, not to the star. This may allow more uniform irradiation of the exomoon while it orbits around its planetary host, thus overcoming the problems the HZ planets may have. In addition, HZ planets around M dwarfs may experience intensive stellar flares typical for M-dwarf stars. The massive Saturn-mass planet may, however, have a strong magnetic field, which could act as an effective shield protecting the habitable moon from stellar activity (Heller & Zuluaga 2013). Of course, these assumptions are highly speculative, while many physical considerations must be taken into account to test how feasible these scenarios are. From the dynamical perspective, long-term stability and tidally induced satellite orbit decay can limit the existence of habitable exomoons around GJ 1148 b (Barnes & O'Brien 2002;Sasaki & Barnes 2014;Kane 2017), and we look further into these possible limitations.

Stability border
The natural scale for the Hill stability of an exomoon is the Hill radius, which for GJ 1148 b is where q b = a b (1 − e b ) is the periastron distance, m b is the planetary mass of GJ 1148 b, and M is the stellar mass of GJ 1148. Grishin et al. (2017) have found an analytic fit for the stability border as a function of arbitrary mutual inclination. The stability border is around 0.5 R Hill for prograde coplanar configurations; decreases to 0.4 R Hill for highly inclined configurations with i ∼ 90 • ; and increases to 0.7-1 R Hill for retrograde coplanar configurations, whose orbits at > ∼ 0.6 R Hill are chaotic and have high eccentricities. Retrograde configurations have long been known to be more stable than prograde ones (Henon 1970;Hamilton & Burns 1991), but are also less likely to occur; thus, we do not A16, page 11 of 18 A&A 638, A16 (2020) study retrograde exomoon orbits and focus only on the co-planar prograde case. For the GJ 1148 system, the outer planet GJ 1148 c could further reduce the stability border, because GJ 1148 c secularly drives the eccentricity and changes the Hill radius of GJ 1148 b. The Hill stability of exomoons was studied in a different context by Hamers et al. (2018) around circumbinary planets, which also further reduced the fraction of stable orbits, but for different reasons involving resonances between the inner binary star and the orbit of the exomoon.
In order to explore the stability of exomoons around GJ 1148 b, we performed N-body integrations of test particles around the temperate planet GJ 1148 b. In this test, we randomly placed 20 000 test particles with semi-major axes in the range 0.001-0.007 au on initially circular orbits. The semi-major axis range was chosen carefully to avoid the planetary Roche limit, inside of which the exomoon cannot exist, and to be inside the planetary Hill radius. The integrations were done using the same Wisdom-Holman algorithm (Wisdom & Holman 1991) that we used in our long-term stability tests based on the MCMC analysis, but modified to handle the evolution of additional test particles in the Jacobi coordinate system (Lee & Peale 2003). We inverted the hierarchy of the system, making the GJ 1148 b planet the central body, orbited by the non-mutually-interacting massless test particles, the GJ 1148 stellar body, and the outer planetary body GJ 1148 c, in that order. The initial parameters of the massive bodies were adopted from the best coplanar edgeon dynamical model. Since the integrations were done in Jacobi frame, the evolution of the three massive bodies will be the same as shown in Figs. 5, 6, and 9 for the best fit, where the central body is the star. The time step chosen for this test is very small, only 0.01 d, which is needed for the accurate N-body integration of the closest test particles, which have periods of only ∼0.8 d. Figure 10 shows the results from the test particle simulations. Each panel of Fig. 10 tracks the evolution of the test particles with a step of 8600 yr to match the snap-shot evolution of the massive bodies in the GJ 1148 system shown in Fig. 9. The blue dot indicates the position of GJ 1148 b assuming that it is the central body, the gray dashed line marks the planetary Roche limit, while the blue dashed line marks the planetary 0.5 R Hill . The 0.5 R Hill limit scales with (1-e b ) due to the dynamical perturbations of GJ 1148 c, i.e., it is farther out when e b ≈ 0, and closer in when e b > 0. Test particles which fall outside the 0.5 R Hill limit are quickly ejected, while those just inside are excited to large eccentricities, and also become unstable at some point during the integration. To summarize, the stability border of the coplanar prograde case is around ∼0.0023 au which is around 0.5 min {R Hill }, where the minimum is obtained for the maximum value of e b during the secular period.
Most likely the larger exomoons will be well within the Hill sphere. Highly inclined exomoons may still be stable if they are within the Laplace radius (Tremaine et al. 2009) where the inner quadrupole of the planet induced by its oblateness overtakes the outer quadrupole of the star. Here J 2 is a dimensionless coefficient that is related to the planet's oblateness, and R b is the radius of the planet. For Jupiter-mass planets, the Laplace radius is a few dozen times the planet's radius, r L ∼ 30 R b . On the other hand, even moderate inclinations slightly outside the Laplace radius can result in unstable orbits either by the inner quadrupole (Tremaine et al. 2009) or by chaotic evolution with an additional fourth body (Grishin et al. 2018a). In A16, page 12 of 18 T. Trifonov et al.: Dynamical characterization of GJ 1148 Fig. 12. Spin evolution of GJ 1148 b with adopted R = 9.44 R ⊕ , k 2 = 0.341, and Q = 18 000 due to star-planet tides (left) and tidal evolution of Mars-like (m = 0.107 M ⊕ , R = 0.53 R ⊕ , k 2 = 0.14, Q = 140) and Titan-like (m = 0.0225 M ⊕ , R = 0.404 R ⊕ , k 2 = 0.589, Q = 58.9) exomoons due to exomoon-planet tides (right). The initial rotational period of GJ 1148 is P rot = 11 h (consistent with Saturn), while of exomoons it is P rot = 24 h (consistent with Mars). Asymptotic rotation at 2/3 of the orbital period is reached for GJ 1148 b after ∼850 Myr. The Saturn-like planet therefore reaches P rot = 27.5 d, which leads to strong orbital decay of the N-body stable exomoon orbits due to tides with the planet. The initial exomoon orbits are adopted from those that survived the N-body stability test in Sect. 5.3.1 (see also Fig. 10).
either case the result can be the ultimate ejection or collision of the exomoon.

Tidal heating
Tidal torques exerted on an exomoon by both the host planet and the host star deform the exomoon. This dissipates heat in the interior of the exomoon, which increases its surface temperature above the value resulting from stellar illumination alone. The amount of heat dissipated depends on the tidal efficiency factor k 2 /Q: this value was determined to be ≈0.015 for Io, and between 0.0026 and 0.0127 for Enceladus (Nimmo et al. 2018).
We wrote an algorithm that computes the surface temperature of an exomoon of given physical and orbital characteristics around any particular exoplanet, then investigated how the surface temperature changes with semi-major axis, eccentricity, and tidal efficiency factor of the moon. We adopted two fiducial moons: one with the physical characteristics of Titan (Titan-like) and one with the physical characteristics of Mars (Mars-like). We varied the semi-major axis within the stability region (i.e., 0.001 AU ≤ a m ≤ 0.002 AU) with the base value at 0.0015 AU; we varied the eccentricity in the range 0.0001 ≤ e m ≤ 0.01, with the base value at 0.001 (one-fifth of the Enceladus eccentricity). We computed surface temperature with and without tidal heating; in the former case k 2 /Q = 0.01 for the Titan-like moon and k 2 /Q = 0.001 for the Mars-like moon (literature values for Mars and other solar system bodies are listed in Lainey 2016), while in the latter case, k 2 /Q was set to zero for both fiducial moons. The results are displayed in Fig. 11. White dots denote the fiducial values; we note that the color bars are not identical for all figures. Light gray regions indicate surface temperatures in excess of the boiling point of water at 1 atmospheric pressure.
In all cases, regardless of tidal heating, a moon at 0.0015 AU with eccentricity around 0.001 has a surface equilibrium temperature of around 255 K. Since this is above the snow line temperature, surface ice may evaporate to form an atmosphere. Hence, in terms of ambient temperature habitable exomoons around GJ 1148 b are possible: if a sufficiently strong greenhouse effect takes place, surface conditions on a hypothetical exomoon GJ 1148 b-I may allow liquid water.

Tidal torque drift due to planet-moon tides
It is important to note that in our stability test, we neglected tidal interactions between the bodies; however, they are important in such a compact system and pose further obstacles for the existence of potentially habitable moons.
There are two critical tidal timescales, which are relevant for exomoons around GJ 1148 b: (i) the spin-orbit synchronization timescale of GJ 1148 b and (ii) the tidal torque drift timescale for the "stable" exomoon due to planet-moon tides. To derive these timescales we use the EQTIDE 7 code (Barnes 2017), which calculates the tidal evolution of two bodies based on models by Ferraz-Mello et al. (2008). To calculate (i), for the starplanet pair we adopted the mean osculating semi-major axis a b = 0.166 au and mean eccentricity e b = 0.361 for GJ 1148 b. Stellar physical parameters were adopted from Table 1, while for the planet the initial spin of P rot = 11h was chosen, and the planetary radius set to R b = 58 230 km (consistent with Saturn). From Zollinger et al. (2017) we adopted Saturn-like physical parameters, such as love number k 2 = 0.341, and dissipation factor Q = 18 000 (see their Table 2). To calculate the tidal evolution in (ii), we adopted the stable test particles (i.e., exomoons) from Sect. 5.3.1, and we adopted Mars-like and Titan-like physical parameters (i.e., for Mars m = 0.107 M ⊕ , R = 0.53 R ⊕ , and for Titan m = 0.0225 M ⊕ , R = 0.404 R ⊕ ), and for both an initial rotation period of 1 d, which is consistent with the rotational period of Mars. For Mars-like moons we adopt k 2 = 0.14 from Murray & Dermott (1999), while for Titan-like moons we adopt k 2 = 0.589 from Iess et al. (2012). For consistency with the adopted k 2 /Q values in Sect. 5.3.2, we adopt Q = 58.9 for Titan, and Q = 140 for Mars. Figure 12 shows the results of our tidal evolution calculations. The left panel of Fig. 12 shows the planetary rotational evolution of GJ 1148 b due to star-planet tides. After ∼850 Myr, GJ 1148 b reaches a rotation period that is 2/3 of the orbital period, and remains there with P rot = 27.5 d. During the integration the planetary semi-major axis and eccentricity are mostly unaffected. An asymptotic rotation period that is shorter than synchronous and 2/3 of the orbital period is expected for e b > ∼ 0.24 in the constant Q tidal model (Goldreich & Peale 1966;Cheng et al. 2014). The time for GJ 1148 b to reach asymptotic rotation is inversely proportional to the initial P rot , as long as the initial P rot is much less than 27.5 d, and it depends on the other parameters of GJ 1148 b according to Eq. (3) of Barnes (2017) and Eq. (15) of Cheng et al. (2014). The rotational period of GJ 1148 b is thus very likely much longer than the orbital periods of the hypothetical exomoons, which could be dynamically stable only with orbital periods between 0.7 and 2 d. The right panel of Fig. 12 shows that the longer rotational period of GJ 1148 b (P rot = 27.5 d) leads to strong orbital decay of the stable exomoon orbits due to tidal interactions with the planet. An exomoon eventually reaches the Roche limit where it is tidally disrupted by the gas giant. Not even one hypothetical "stable" exomoon in the context of Sect. 5.3.1 had survived this test. The maximum time a Mars-like exomoon could survive is ∼55 M yr, while for Titanlike moons the maximum survival time is longer, ∼255 M yr. The latter is longer by roughly the mass ratio of Mars to Titan, which can be understood from Eq. (2) of Barnes (2017) and Eq. (16) of Cheng et al. (2014). These timescales are optimistic since the orbital decay would start before the planet reaches the asymptotic spin state. In both cases the survival times are much shorter than the age of the system. Therefore, given the relatively fast orbital decay in the small stable region around the planet, we conclude that exomoons around GJ 1148 b are unlikely to exist.
It should be noted that Martínez-Rodríguez et al. (2019) also considered the tidal evolution of a hypothetical exomoon of GJ 1148 b and reported a timescale for orbital migration larger than 14 Gyr, a dramatically different result. The source of the discrepancy is that they assumed a much longer orbital period for the moon than we do (8 days versus < 2 days), which changes the tidal evolution by orders of magnitude. In addition, they did not calculate the stability limit for the moon's orbit.

Conclusions and discussion
We present an updated orbital solution and dynamical analysis of the GJ 1148 M-dwarf planetary system based on new CARMENES optical RV data. We also present CARMENES near-IR RVs for GJ 1148, which are in excellent agreement with the available optical data from HIRES and CARMENES. The GJ 1148 system is in a peculiar configuration consistent with a pair of Saturn-mass gas giants, which exhibit large eccentricity oscillations over a secular timescale of ∼67 000 yr. Even so, the system is well separated and remains stable for a large set of coplanar and mutually inclined configurations. Characteristic for the GJ 1148 system is the evident apsidal alignment due to the large eccentricities, but overall we did not find any high-order mean-motion resonant commensurabilities in which the system could be trapped. It is possible that the observed high eccentricities are the result of planet-planet scattering events, which could possibly have occurred in the early phases of the history of the system.
We find that the existence of potentially habitable exomoons around GJ 1148 b is unlikely due to the narrow stability region around the close and eccentric planet. Additionally, the tidally driven orbit decay timescales are much shorter than the age of the GJ 1148 system. When considering both disturbances by an outer planet in a system with significant eccentricities and tidal effects, there may be no space left for moons. Even if exomoons had resided in the stable region around the inner planet in the early phase of the history of the system, they would eventually spiral towards the planet.
TESS will observe GJ 1148 in Sector 22 from 18 February 2020 to 18 March 2020. Given the stellar radius estimated by Schweitzer et al. (2019) and the planetary orbital eccentricities and semi-major axes we obtain in this work, we find a geometric transit probability of ∼1% for GJ 1148 b and ∼0.1% for GJ 1148 c (see, e.g., Barnes 2007). Given the short duration (∼27 days) of the TESS photometric observations, we can expect to detect at most a single transit event for each planet. Based on the posterior distribution of the orbital parameters from this work, we predict that a potential transit of GJ 1148 b will occur on 13 February 2020 around 11:30:00.0 UT ± 4h, and the next one will be on 25 March 2020 ∼20:40:00.0 UT ± 4h. So both are unfortunately outside the expected TESS window. For GJ 1148 c, the situation is even worse with the closest potential transit calculated to occur 84 ± 12 d before TESS starts the observations of Sector 22. Therefore, there is no chance to observe potential transits of GJ 1148 b & c with TESS. Nevertheless, TESS could reveal shorter period, low-mass transiting planets whose RV signature could be below the RV jitter we record with HIRES and CARMENES, and thus remain undetected by the Doppler method.