HD 191939 revisited: New and refined planet mass determinations, and a new planet in the habitable zone

HD 191939 (TOI-1339) is a nearby (d=54pc), bright (V=9mag), and inactive Sun-like star (G9 V) known to host a multi-planet transiting system. Ground-based spectroscopic observations confirmed the planetary nature of the three transiting sub-Neptunes (HD 191939 b, c, and d) originally detected by TESS and were used to measure the masses for planets b and c with 3$\sigma$ precision. These previous observations also reported the discovery of an additional Saturn-mass planet (HD 191939 e) and evidence for a further, very long-period companion (HD 191939 f). Here, we report the discovery of a new non-transiting planet in the system and a refined mass determination of HD 191939 d. The new planet, HD 191939 g, has a minimum mass of 13.5$\pm$2.0 M$_\oplus$ and a period of about 280 d. This period places the planet within the conservative habitable zone of the host star, and near a 1:3 resonance with HD 191939 e. The compilation of 362 radial velocity measurements with a baseline of 677 days from four different high-resolution spectrographs also allowed us to refine the properties of the previously known planets, including a 4.6$\sigma$ mass determination for planet d, for which only a 2$\sigma$ upper limit had been set until now. We confirm the previously suspected low density of HD 191939 d, which makes it an attractive target for attempting atmospheric characterisation. Overall, the planetary system consists of three sub-Neptunes interior to a Saturn-mass and a Uranus-mass planet plus a high-mass long-period companion. This particular configuration has no counterpart in the literature and makes HD 191939 an exceptional multi-planet transiting system with an unusual planet demographic worthy of future observation.


Introduction
The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) is a NASA-sponsored space telescope whose original mission was a 2-yr full-sky survey searching for transiting planets. One of TESS main scientific goals is to look for small planets (R P < 4 R ⊕ ) around bright stars suitable for radial velocity (RV) follow up and atmospheric characterisation. Since the beginning of operations in 2018, TESS has discovered several multi-planetary transiting systems around bright host stars (e.g. Dragomir et al. 2019;Günther et al. 2019;Quinn et al. 2019;Luque et al. 2021). This type of system is an excellent laboratory for planetary astrophysics. Multi-planet systems share the same initial conditions (e.g. protoplanetary disc), allowing for comparison between sibling planets, and also the testing of planet formation and evolution processes and theories.
Here, we focus on the HD 191939 system, one of those multi-planetary transiting systems with small planets discovered by TESS. HD 191939 is a bright (V = 9 mag), nearby (d = 54 pc), inactive solar-like star (G9 V). Badenas-Agusti et al. (2020, hereafter BA20) confirmed the presence of three transiting sub-Neptune-sized planets, HD 191939 b, c, and d, with periods of 8.88, 28.58, and 38.35 days, respectively. Their radii are very similar, ranging 3.16-3.42 R ⊕ . The brightness of the host star made long-term RV monitoring campaigns with different high-resolution spectroscopy facilities feasible, allowing the determination of the planetary masses. The physical properties of the planetary system were studied by Lubin et al. (2022, hereafter L22), finding masses of 10.4±0.9 M ⊕ and 7.2±1.4 M ⊕ for planets b and c, respectively. L22 only presented an upper limit of 5.8 M ⊕ (2σ confidence) for planet d. Furthermore, L22 found evidence for two extra planets via RV measurements: planets e and f. HD 191939 e has a period of 101 days and a minimum mass of 108±3 M ⊕ . The long-term trend showed by the RVs was related to the presence of HD 191939 f, a high-mass planet with an unconstrained period of between 1700 and 7200 days. In Table 1, we compile a comprehensive list of HD 191939 stellar properties from the literature.

Observations
To remove some additional variability present in the PDC-SAP and to save computational time fitting the 11 TESS sectors, we detrended the PDC-SAP light curves, with the planetary transits masked, performing a Gaussian process (GP) regression model using a Matern kernel 3/2 from celerite (Foreman-Mackey et al. 2017). The model only considered a relative flux offset and a jitter term, and two GP hyperparameters, which are shared between the different TESS sectors. The priors and posteriors from the detrended process are shown in Table A Quirrenbach et al. 2014Quirrenbach et al. , 2020 located at the Calar Alto Observatory, Almería, Spain. CARMENES has two spectral channels: the optical channel (VIS), which covers the wavelength range from 0.52 to 0.96 µm with a resolving power of R = 94 600, and the nearinfrared (NIR) channel, which goes from 0.96 to 1.71 µm with a resolving power of R = 80 400. The star was monitored from 6 November 2019 to 4 September 2021. During this time, we obtained 138 high-resolution spectra.
The observations were carried out as part of observing programs F19-3.5-014, S20-3.5-011 (PI: Nowak), and F20-3.5-013 (PI: Luque). The exposure times were set to 900 s, leading to a signal-to-noise ratio (S/N) per pixel of 41-188 at 7370 Å. The observations were reduced using the CARMENES pipeline caracal (Caballero et al. 2016) and we processed the VIS and NIR spectra with serval 2 (Zechmeister et al. 2018), which is the standard CARMENES pipeline to derive relative RVs and several activity indicators using template matching: chromatic RV index (CRX), differential line width (dLW), and Hα, Na D1 and Na D2, and Ca II IRT line indices. We also used the RACCOON code 3 (Lafarga et al. 2020) to measure the CCF_FWHM, CCF_CTR, and CCF_BIS spectral activity indicators via cross-correlation. In the analysis presented here, we used the activity indicators from VIS and NIR extracted with serval and RACOON, and the RVs measured from CARMENES VIS spectra with serval. Because the precision in the RVs obtained from the VIS is higher than that obtained with the NIR, we only used the CARMENES VIS RVs, which have smaller error bars. CARMENES VIS RVs are corrected using measured nightly zero-point corrections as discussed in Trifonov et al. (2020;shown in

High-resolution spectroscopy with HARPS-N
HD 191939 was observed with the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N; Cosentino et al. 2012) mounted on the 3.6 m Telescopio Nazionale Galileo (TNG) in Roque de los Muchachos Observatory, La Palma. The star was monitored from 30 May 2020 to 13 September 2021. During this time, we obtained 42 highresolution (R = 115 000) spectra.
The observations were carried out as part of observing programs CAT19A_162 program (PI: Nowak), ITP19_1 (PI: Pallé) and CAT21A_119 (PI: Nowak). The exposure times varied from 284 to 1800 s, depending on weather conditions and scheduling constraints, leading to a S/N per pixel of 27-134 at 5500 Å. The spectra were extracted using the off-line version 3.7 of the HARPS-N DRS pipeline (Cosentino et al. 2014). Doppler measurements (absolute RVs) and spectral activity indicators (CCF_FWHM, CCF_CTR, CCF_BIS, and Mount-Wilson S-index) were measured using an online version of the DRS, the YABI tool 4 , by cross-correlating the extracted spectra with a G2 mask (Baranne et al. 1996). We also used the serval code to measure relative RVs by template matching, CRX, dLW, and Hα and sodium Na D1 and Na D2 indexes. The uncertainties of the RVs measured with serval are in the range of 0.5-3.1 m s −1 , with a mean value of 1.07 m s −1 . In the analysis presented here, we used the relative RVs measured from HARPS-N spectra with serval (shown in Fig. 1).

High-resolution spectroscopy with APF and HIRES
L22 also performed a ground-based follow-up campaign with two different high-resolution spectrographs. They obtained 73 RV measurements with HIRES and 107 RV measurements with the Automated Planet Finder (APF, Vogt et al. 2014) telescope. In total, these observations covered a baseline of 415 days, and their details are explained in Sect. 2.2 of L22. The time series of APF and HIRES are displayed in Fig. 1 along with the CARMENES and HARPS-N RVs.

Analysis and results
3.1. Stellar rotation and activity indicators BA20 and L22 reported that HD 191939 is a slow rotator star with low or null stellar and chromospherical activity. BA20 4 Available at http://ia2-harps.oats.inaf.it:8000 derived a P rot / sin i = 79±66 days, where the large uncertainties come from the large error on the v sin i.
We searched for modulations in the different activity indices derived from CARMENES and HARPS-N spectra using the generalised Lomb-Scargle (GLS; Zechmeister & Kürster 2009) periodogram 5 . We computed the theoretical false-alarm probability (FAP) as it is described in Zechmeister & Kürster (2009). The GLS periodograms of the activity indices are shown in Fig. A.1. None of the indices present a significant peak at either the periodicities of the known planets or at the ∼300 days signal, attributed to HD 191939 g. As is expected for a low-activity star, the GLS periodograms remain below the 10% level of significance. Only some CARMENES indices show a broad peak near 400 days that is not detected in the HARPS-N indices, and CARMENES Ca II IRT 1 and HARPS-N Hα display significant periodicities near 100 days and 200 days, respectively. However, none of those peaks have a counterpart in the RVs periodogram analyses.

Radial velocity analysis
The analyses presented in this work combine the APF and HIRES RVs from L22 with those obtained with CARMENES and HARPS-N. With a total of 362 RVs covering a baseline of 677 days (Fig. 1), we were able to improve the planetary mass determination but also look for the presence of planets with periods between those of planets e (101 days) and f (>1700 days).
We analysed the planetary signals in the RVs computing the GLS periodogram and modelling the detected signals with juliet 6 (Espinoza et al. 2019). This python library is based on other public packages for transit light curve (batman, Kreidberg 2015), RV (radvel, Fulton et al. 2018), and GP (george, Ambikasaran et al. 2015;celerite, Foreman-Mackey et al. 2017) modelling. Juliet uses nested sampling algorithms (dynesty, Speagle 2020; MultiNest, Feroz et al. 2009;Buchner et al. 2014) to explore all the parameter space and compute the Bayesian model log-evidence (ln Z), which allows us to compare models with different numbers of free parameters. If the difference between two models, for example M1 and M2, is ∆ ln Z = ln Z M2 − ln Z M1 > 5, then the M2 model is strongly preferred statistically over the M1 model (Trotta 2008). If ∆ ln Z ≤ 1, the two models are statistically indistinguishable  (e) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 days, 28.6 days, 38 days (vertical cyan line), and 101 days signals. ( f ) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 days, 28.6 days, 38 days, 101 days, and 300 days (vertical orange line) signals. In all panels, the 10%, 1%, and 0.1% FAP levels are indicated by dotted, dash-dotted, and dashed grey horizontal lines, respectively. The vertical black dotted line indicates the dataset baseline. We highlight the different scale of the y axis in each panel. and the preferred one is the simplest model with the least free parameters. We consider that M2 has moderate evidence over M1 for intermediate cases (∆ ln Z ∼ 2.5).
Because the main purpose of this preliminary study is to look for additional signals, we considered circular orbits for simplicity. Eccentric models are explored during the joint fit (see Sect. 3.3) after exploring the signals present in the data. For this RV-only analysis, we fixed the period (P) and the central time of transit (t 0 ) for the three transiting planets based on a photometriconly analysis. The precision derived for P and t 0 from the light curves is significantly higher than from the RVs alone. Fixing these two parameters saves computational time without an impact on the RV modeling. For the signals with no counterpart in the photometry, we set normal priors for P and uniform priors for t 0 and we set uninformative priors for the semi-amplitude (K) of the fitted signals. For each spectrograph, we also included an instrumental jitter term (σ) and a systemic velocity (γ) term. The procedure described below to fit the periodicities detected in the periodograms is illustrated in Fig. 2. The following points also refer to the panels of Fig. 2.
(a) As in L22 and to save computational time, we modelled and subtracted the long-term RV trend detected by the four instruments with a linear (γ) term and a quadratic (γ) term. The linear and quadratic term model is statistically preferred (∆ ln Z > 5) over an only linear or only quadratic term model. The conspicuous signal at 101 days due to planet e dominates the RV periodogram.
(b) After fitting planet e model, the planet b periodicity is clearly seen in the periodogram of the residuals.
(c) When planet b is removed, the planet c signal is the most significant peak in the periodogram.
(d) After fitting planet c, all periodicities in the periodogram of the residuals remain well below the 10% level of significance, except for a peak at ∼300 days (FAP < 1%). As in Fig. 1 from L22, the signature of planet d at 38.35 days is not detected or significant in the RVs.
(e) When all the previously known planets are removed, the periodogram of the residuals is still dominated by the signal at ∼300 days which slightly increased its significance until FAP ∼ 0.1%. The ∼300 days peak is already visible and significant after removing planet e.
(f) The periodogram after fitting the ∼300 days signal is flat without significant peaks. We refer to the ∼300 days signal as HD 191939 g. We analyse this signal in detail below.
To crosscheck our results, we also analysed the RV dataset using Exo-Striker 7 (Trifonov 2019), obtaining similar results. After fitting the transiting planets b, c, and d as well as planets e and f, the GLS periodogram only showed a peak around ∼300 days, which is the signal we named HD 191939 g. Furthermore, the signal at 17.7 days found by L22 (Sect. 8 therein) is observed in the residuals in Fig. 2 but at a very low significance level (FAP≫10%). Our RV dataset does not support the scenario of a non-transiting planet with a period between transiting planets c and d.
We repeated the analysis presented above only with the APF and HIRES datasets, obtaining similar results to L22; that is, there is no statistically significant evidence (FAP>10%) for an additional signal at ∼300 days (see Fig. A.3). The non-detection of that signal in L22 may be due to the lower number of RVs used and the shorter baseline of those observations. The baseline of the 180 APF and HIRES RVs is 415 days, which is less than 1.5 periods of the ∼300 days signal. Here, we used 362 RVs with a baseline of 677 days, which cover two complete periods.
To unveil the nature of the ∼300 days signal, we computed a set of models that simultaneously fit the known planets and also account for the additional signal with two different strategies: fitting the signal by adding a GP term or with a Keplerian orbit. We considered three different GP kernels: exponential (GP exponential ), Matern 3/2 (GP Matern ), and quasi-periodic (GP qp ). When fitting the signal with a Keplerian orbit or with a GP qp kernel, we tested two different priors for the period parameter or rotational hyperparameter: an uninformative prior and a normal prior centred at 300 days. The priors and posteriors of the hyperparameters used in the GP models are shown in Table A.2. Table 2 presents the Bayesian log-evidence for the explored models and the measured K for each of the fitted signals. The derived Ks for the Table 2. Comparative between Bayesian log-evidences (∆ ln Z) and planet semi-amplitudes (K p ) for the different explored models.

Model
∆ ln Z K b (P = 8.8 days) K c (P = 28.6 days) K d (P = 38 days) K e (P = 101 days) K g (P ≃ 300 days) Notes. We used the ln Z from the four-planet model as a reference. The adopted model used in the joint fit is marked in boldface (see Sect. 3.2 for details about the selection of the final model). The last row shows the K p from the joint fit in Sect. 3.3 for illustrative purposes. (a) Using an uninformative prior for the periodic parameter. (b) Using a normal prior for the periodic parameter centred at 300 days. different models are consistent within errors, ensuring that the planet mass is not model-dependent or affected by the ∼300 days signal.
The five-planet model and the GP qp model are statistically preferred over the four-planet model, the GP exponential model, and the GP Matern model. The ∆ ln Z between the five-planet model and its analogous GP qp model is less than 2. Therefore, the GP qp model is moderately preferred over the five-planet model. However, the GP qp model has one free parameter more than the five-planet model. Moreover, closer inspection of the posterior distributions of the GP qp hyperparameters shows that the P rot is not constrained (see Fig. A.2). In the uninformative GP qp model, the P rot posterior distribution is mainly flat. When we used a normal prior for P rot , the posterior distribution is equal to the input prior, which is in contrast to the period of the additional signal, which is well determined in both five-planet models.
Therefore, we chose the five-planet model over the GP qp model due to its more physical plausibility. The complete RV model includes quadratic and linear terms, planets HD 191939 b, c, d, and e, and the new planet HD 191939 g.

Joint fit
We simultaneously modeled the detrended two-minute cadence TESS photometry and the APF, HIRES, CARMENES, and HARPS-N RVs using juliet to refine the parameters for the HD 191939 system. For the joint fit, we considered transits for the planets b, c, and d, and we adopted a five-planet model with a quadratic and linear trend for the RVs.
We adopted a quadratic limb-darkening law for the TESS light curve. The limb-darkening coefficients were parameterised with the uniform sampling prior (q 1 ,q 2 ) introduced by Kipping (2013). Additionally, rather than directly fitting the impact parameter of the orbit (b) and the planet-to-star radius ratio (p = R p /R ⋆ ), we considered the uninformative sample (r 1 , r 2 ) parameterization introduced in Espinoza (2018). The parameters r 1 and r 2 ensure a full exploration of the physically plausible values of p and b, with uniform priors sampling. We used the value and uncertainties of the stellar density (ρ ⋆ ) from L22 in Table 1 to set a normal prior for ρ ⋆ . We fixed the dilution factor to 1 based on BA20 and we added a relative flux offset (µ) and a jitter term (σ) to TESS data.
Systems with multiple transiting planets, such as HD 191939, normally present low eccentricities, but not necessarily zero (Van Eylen & Albrecht 2015; Xie et al. 2016;Hadden & Lithwick 2017). We therefore considered Keplerian orbits for the five planets with a beta prior distribution with α = 1.52 and β = 29 for the orbital eccentricity ecc (Van Eylen & Albrecht 2015; Van Eylen et al. 2019). We also computed a joint fit with circular orbits. However, the eccentric joint-fit model is statistically preferred over the non-eccentric one (∆ ln Z = ln Z ecc − ln Z no ecc > 8) and the results from both models are coincident within their uncertainties.
The priors used in the joint fit are listed in Table A.3. The posterior distributions and the derived parameters for the planetary system are reported in Table 3. As is shown in Table 2, the semi-amplitudes obtained from the joint fit are consistent with the RV-only models explored in Sect. 3.2. The stellar density is consistent with those derived by L22 and by BA20 within 1σ. The adopted stellar properties used to derive the planetary parameters are the ones from L22 (see Table 1) for a better comparison with the results presented there. The best-fit models for phase-folded light curves and phase-folded RVs are shown in Figs. 3 and 4, respectively. Photometric and RV time series along with the best-fit models are shown in Figs. A.4-A.6, respectively.
We checked for transit-like events in the TESS PDC-SAP and SAP for the non-transiting planets with no results. We estimated the flux decrease (∆F ≃ (R p /R ⋆ ) 2 ) and transit duration t T of planets e and g with their predicted radii (see Sect. 4.1). HD 191939 e would produce a flux decrease of ∼1.6% over ∼8 h. According to the derived ephemeris, a transit of planet e was expected during Sector 16 (see Fig. A.4) and is clearly not detected. Because the predicted radius for HD 191939 g is similar to that of the inner planets, the flux decrease for planet g would be similar (∼1300 ppm). However, the transit would span ∼12 h due to its large period. Planet g was expected to transit at some point in Sectors 18 and 19 (see Fig. A.4). Although there is no clear evidence of a complete transit, ingress, or egress during these sectors, the 1σ uncertainty for the transit midpoint (∼16 days) comprises observing gaps where the transit could have happened. These observing gaps represent 18% of the ± 1σ expected transit time region. Upcoming TESS sectors will allow us to confirm whether we were unlucky or HD 191939 g does not transit, as in HD 191939 e.    below 13 M J . L22 performed a joint RV and astrometric analysis to impose some limits on the properties of planet-candidate f. From that analysis, the period should be between 1700 and 7200 days, which is still much longer than the baseline of our new combined RV dataset.
We also tried to analyse the properties of planet-candidate f by fitting the RV long-term trend with a Keplerian signal instead of a quadratic term plus a linear term. Although the RV time series from this work have a longer base line, they do not show a peak to peak of HD 191939 f signal that could help to constrain its period and mass. Therefore, the values reported for planet-candidate f (see last column in Table 3) should be considered as updated upper and lower limits. The best-fit model for the RVs (Fig. A.7) and the phase-folded RVs for planetcandidate f (Fig. A.8)  Earth-like composition 100%H 2 O 5%H 2 + 95% Earth-like 5%H 2 + 95% Water-rich with the 1700-7200 days period range, and the semi-amplitude of ≳36 m s −1 is consistent with the previous lower limit of >23 m s −1 from APF and HIRES RVs. The properties of the planetary candidate HD 191939 f are still not well determined, and therefore further observations of HD 191939 are needed to sample its long period and better constrain its properties.

Discussion
By compiling several RV observations, we detect the signal of a new, likely non-transiting planet, namely HD 191939 g, with a period of 284 +10 −8 days. We derived a minimum mass of 13.5±2.0 M ⊕ with a precision of 15%. Moreover, we refined the planet properties of the previously known planets, HD 191939 b, c, d, and e. In particular, we were able to determine the semiamplitude and planetary mass of HD 191939 d with a 4.6σ level of significance and confirm the low density of this sub-Neptune planet. Figure 5 puts in context the different planets of the HD 191939 system as compared to the known exoplanets with masses and radii measured with a precision of better than 30%, from 1 M ⊕ up to 1.5 M J , and 1 R ⊕ up to 1.5 R J , along with theoretical composition models of Zeng et al. (2019)

HD 191939 g: a new planet
With an orbital period of ∼280 days and a minimum mass of 13.5 M ⊕ , HD 191939 g joins the selected group of exoplanets that could only be detected thanks to a large number of RV measurements spanning a relatively wide time interval. Because the RV method is more sensitive to shorter period planets and also to the more massive ones, there are only a few long-period planets with intermediate ( Because HD 191939 g and e are only detected in the RV measurements and TESS photometry does not show evidence of their transits, we were not able to determine radii. We therefore marked their positions with vertical bands in the mass-radius diagram (Fig. 5 bottom). However, their radii can be forecasted using empirical mass-radius relations for planets. We used the probabilistic planet mass-radius relation given in Chen & Kipping (2017) via its python implementation 9 . The code predicted planetary radii for planets g and e of 3.7 +1.5 −1.1 R ⊕ and 12.9 +4.7 −3.9 R ⊕ , respectively. From their minimum masses and forecasted radius, the expected mean bulk densities for planets g and e are 1.5 +2.7 −0.8 g cm −3 and 0.30 +0.70 −0.16 g cm −3 , respectively. To crosscheck these forecasted results for planet g, we also estimated its radius using the mass-radius relation for sub-Neptune-sized planets of Wolfgang et al. (2016). This method predicts a planetary radius of ∼3.4 R ⊕ for planet g, which is consistent with the above estimation. If we extrapolate this mass-radius relation to planet e, the predicted radius is about ∼17 R ⊕ , also falling within the uncertainties. However, we stress that these estimated values are merely illustrative, and should not be considered as the actual planetary radii and densities.
HD 191939 g is the only planet in the system in the conservative habitable zone (HZ) of the star; that is, by definition its T eq is compatible with the presence of liquid water (T ∈ [273, 373] K). With a semi-major axis of ∼0.82 AU, planet g is in the outer edge of the HZ, which we set at ∼0.44-0.84 AU. Figure 6 displays a face-on view of the HD 191939 system, where the HZ of the star is marked. However, we stress that despite being in the HZ, HD 191939 g, being a gaseous planet, cannot be considered as a habitable planet.
Long-period, intermediate-mass planets are located in a lonely region of the mass-period diagram, with HD 191939 g at the centre of this group (see Fig. 5). These objects have the commonality that they are outer planets in their respective planetary systems and their orbits are in or near the HZ of the host star. In the mass-radius diagram, planets with masses of ∼13.5 M ⊕ are above the Earth-like composition line, supporting the idea that HD 191939 g is likely a gaseous planet. Moreover, its 1σ mass uncertainty overlaps with Kepler-90 g and HIP 41378 f mass determinations (see Fig. 5). These planets are two of the lowest density (puffy) planets known, with ρ p = 0.15±0.05 g cm −3 (Liang et al. 2021) and ρ p = 0.09±0.02 g cm −3 (Santerne et al. 2019), respectively.

HD 191939 d: a puffy planet
L22 already noted that HD 191939 d was probably a lowdensity planet. Here, we derived a bulk planetary density of ρ d = 0.57±0.13 g cm −3 , confirming this previous result with better uncertainty. With such low density, planet d is at the edge of the planetary mass-radius distribution (Fig. 5 bottom). Thermal expansion of the atmosphere is a possible mechanism to explain planet inflation leading to puffy atmospheres, such as in the case of ultra-hot Jupiters. However, given the relatively cold equilibrium temperature of HD 191939 d (T eq ≃ 540 K), this explanation is unlikely. The nearest planet to HD 191939 d in the mass-radius diagram is Kepler-79 e (KOI-152 e; Jontof-Hutter et al. 2014), which has very similar properties (3.49±0.14 R ⊕ , ρ p = 0.53±0.15 g cm −3 , T eq ≃ 480 K). The most relevant characteristic of the Kepler-79 system is the low density of its planets, whose masses were calculated from transit-time variations. Their densities range between ρ p = 0.09 and 1.43 g cm −3 and the densest planet is the innermost one. This similarity with the HD 191939 transiting planets reinforces the hypothesis that the non-transiting planets of the system are also of a gaseous-like composition.
The brightness (J = 7.6 mag) and low level of stellar activity of the host star offer an excellent opportunity to inspect and study the atmosphere of a puffy planet. To quantify the viability of these observations, we computed the transmission spectroscopy metric (TSM) proposed by Kempton et al. (2018). The estimated TSM value for HD 191939 d is 227, which is well above the threshold of 90 indicated by Kempton et al. (2018) and planets b and c (TSM b = 153; TSM c = 107). Moreover, HD 191939 d has a much better TSM value than other known puffy planets such as HIP 41378 d (TSM = 71), HIP 41378 e (TSM = 57), Kepler-79 planets (TSM = 7-60), or Kepler-90 g (TSM = 27). We note that the TSM is simply proportional to the expected transmission spectroscopy S/N, assuming standardised planetary atmosphere models (e.g. clear atmosphere with solar composition). Observational surveys do not support a strong correlation between expected transmission spectroscopy S/N and actual atmospheric detectability (Tsiaras et al. 2018).
We searched for planets with a radius of ∼3 R ⊕ and/or a T eq ∼ 500-600 K in the database of exoplanet atmospheric observations of ExoAtmospheres 10 . Only the warm sub-Neptune GJ 1214 b (2.74±0.05 R ⊕ , T eq ≃ 596 K; Cloutier et al. 2021) fitted our conditions, although it is a denser planet (ρ p = 2.20±0.17 g cm −3 ). For GJ 1214 b (TSM = 440), only a tentative detection of He I could be set recently (Orell-Miquel et al. 2022) after many non-detection results (Bean et al. 2010;Kreidberg et al. 2014;Petit dit de la Roche et al. 2020;Kasper et al. 2020). When we looked for puffy planet observations, we found that the atmosphere of HIP 41378 f was analysed via transmission spectroscopy at low resolution with the Hubble Space Telescope (HST). However, HIP 41378 f, with a higher TSM (=342) than HD 191939 d, showed a featureless NIR spectrum with a median precision of 84 ppm (Alam et al. 2022).
We explored the potential of HD 191939 d for transmission spectroscopy with the James Webb Space Telescope (JWST) through spectral simulations for a range of atmospheric scenarios. We adopted TauREx 3 (Al-Refaie et al. 2021) to compute our set of model atmospheres using the atmospheric chemical equilibrium (ACE) module (Agúndez et al. 2012), including collisionally induced absorption by H 2 --H 2 and H 2 --He (Abel et al. 2011(Abel et al. , 2012Fletcher et al. 2018), and Rayleigh scattering. The benchmark model assumes a clear atmosphere with solar composition, which displays the largest spectral features. The other models include the dampening effects on the transmission spectrum due to enhanced metallicity or haze in the HD 191939 d atmosphere. The haze was modelled with TauREx 3 using a Mie scattering contribution with the formalism of Lee et al. (2013). A super-solar metallicity is indeed predicted for low-mass, lowdensity planets such as HD 191939 d based on the core accretion theory of planet formation (Fortney et al. 2013;Thorngren et al. 2016). The equilibrium temperature of 540 K also favours the formation of high-altitude photochemical haze in the HD 191939 d atmosphere (Gao & Zhang 2020;Ohno & Tanaka 2021;Yu et al. 2021).
We used ExoTETHyS (Morello et al. 2021) to simulate the corresponding JWST spectra, as observed with the NIRISS-SOSS (0.6-2.8 µm), NIRSpec-G395M (2.88-5.20 µm), and MIRI-LRS (5-12 µm) instrumental modes. The procedure to select the spectral bins and estimate the error bars was identical to that of previous papers (e.g. Espinoza et al. 2022;Luque et al. 2022). In particular, we obtained error bars of 10-12 ppm per spectral point for the NIRISS-SOSS and NIRSpec-G395M  Notes. These include clear or hazy atmospheres with scaled solar metallicities. We adopted the formalisms of Lee et al. (2013) for the hazy models, where α denotes the particle size in µm, the mixing ratio is χ c = 10 −12 , and the extinction coefficient is Q 0 = 40. Lines in bold correspond to the synthetic spectra shown in Fig. 7. Nominal wavelength ranges of (a) HST WFC3-G141 scanning mode, (b) JWST NIRISS-SOSS mode, (c) JWST NIRSpec-G395M mode, and (d) JWST MIRI-LRS mode.
modes at median resolving power of R ∼ 50, and of 27 ppm for the MIRI-LRS bins with sizes of 0.1-0.2 µm. We note that the predicted error bars in the NIR spectrum of HD 191939 d are seven times smaller than those reported for HIP 41378 f by Alam et al. (2022). Figure 7 shows the synthetic transmission spectra for three selected atmospheric configurations, one of which with simulated JWST observations overplotted. These spectra exhibit strong absorption features due to H 2 O and CH 4 , which are an order of magnitude larger than the predicted error bars. Table 4 reports the amplitudes of spectral modulations at low resolution (R ∼ 170) for the full set of synthetic spectra within the nominal wavelength ranges of HST and JWST instrumental modes. Higher metallicities lead to smaller absorption features over the entire spectral range, as they increase the mean molecular weight, thereby reducing the atmospheric scale height. The haze mostly affects the visible and NIR portion of the spectrum, flattening the absorption features and introducing a possible slope. The mid-IR spectrum is less severely affected by haze, but also depends on its physical properties. Based on the predicted spectroscopic amplitudes, even a clear atmosphere with 1000× solar metallicity or a hazy one with 100× solar metallicity surrounding HD 191939 d would be detectable with a single JWST visit. Multiple instruments can break the degeneracy between haze or clouds and metallicity effects. Similar considerations could also apply to other puffy planets with a flat near-infrared spectrum observed with HST WFC3-G141 (Alam et al. 2022;Chachan et al. 2020;Libby-Roberts et al. 2020), albeit with quantitative differences.

Architecture of the planetary system
BA20 and L22 already noted that the transiting planets are close to a near mean motion resonance of 1:3:4 (P b = 8.88 days, P c = 28.58 days, P d = 38.35 days). HD 191939 e, with a much longer period, seems disconnected from that resonance chain. However, the discovery of HD 191939 g reveals that the nontransiting planets of the system appear to be in a period ratio of 1:3 (P e = 101 days, P g ≃ 280 days). There are other cases of multi-planetary systems where inner planets are gathered in a different resonance chain from the outer ones: Kepler-90 planets (Cabrera et al. 2014) are in 2:3:4 and 4:5 periods, and HIP 41378 planets (Vanderburg et al. 2016;Santerne et al. 2019) are in 1:2:4 and 3:4:6 periods. Furthermore, L22 noted that planetary systems hosting puffy planets tend to have their planets in resonance (e.g. Kepler-79, Jontof-Hutter et al. 2014;Kepler-51, Masuda 2014;Kepler-87, Ofir et al. 2014), which seems to also be the case for HD 191939. L22 searched in the literature for planetary systems similar to that described in their work; their best match was Kepler-68. Mills et al. (2019) described this system as two sub-Neptunes interior to a 634d period Jovian planet, and with strong evidence for an object with >0.6 M J in a very long-period orbit (≫3000 days). Another similar system analysed in Mills et al. (2019) is Kepler-65: three sub-Neptunes near orbital resonance of 1:3:4 interior to a non-transiting planet with a mass of 212 M ⊕ and a period of 259 days. However, the detection of a Uranus-mass planet between the warm Saturn and the massive long-period planet makes the HD 191939 system more exceptional.
We performed a new search in the NASA Exoplanet Archive 11 database looking for systems with 2-4 intermediate planets (2-25 M ⊕ or 2-8 R ⊕ ) interior to a gas-giant planet (>50 M ⊕ or >8 R ⊕ ) plus a planet with comparable properties to those of the inner ones. Although a system with three intermediate planets interior to a gas giant plus another intermediate planet was not found, our search returned one system that suited our initial conditions: KOI-94 (Weiss et al. 2013).
KOI-94 has two inner planets with masses of 10.5±4.6 M ⊕ and <21.3 M ⊕ , and then a warm Saturn-like planet (106±11 M ⊕ , 11.2±1.1 R ⊕ ) and an outer planet slightly more massive than the inner ones (35 +18 −28 M ⊕ , 6.6±0.6 R ⊕ ), but without evidence for a long-period companion like HD 191939 f. KOI-94 planets are gaseous with low densities ranging between ∼0.35 and 1 g cm −3 (except for KOI-94 b; 10.1±5.5 g cm −3 ), and the mass of KOI-94 d is consistent with the minimum mass derived for HD 191939 e. KOI-94 b, c, and d are close to a mean motion resonance of 1:3:6 (P b = 3.7 days, P c = 10.4 days, P d = 22.3 days), with the giant planet also close to 5:2 with the outermost planet KOI-94 e (P e = 54 days). In both systems, a more massive planet divides the planets with similar masses (and likely similar characteristics). This mass distribution across the planets is illustrated in Fig. 8. Although the KOI-94 system is more compact than HD 191939, the semi-major axis (a) scaled to the semimajor axis of the most massive known planet (KOI-94 d and HD 191939 e, respectively) shows that the a of the outer planet (KOI-94 e and HD 191939 g, respectively) is approximately twice that of the massive one and the a of the inner planets are ∼0.5 that of the massive one. The planets in these systems present resonances between periods with the massive planet linked to the outer one.
Moreover, in addition to the spectral type of the host star, the HD 191939 system presents some similarities (excluding 11 https://exoplanetarchive.ipac.caltech.edu/ planet f) to our own Solar System. The smaller planets are interior to the massive planet and the intermediate-mass planets are further out relative to the massive one. Still, the planetary system HD 191939 is more compact, with all the constrained planets in the system confined within ∼0.82 au (HD 191939 g semimajor axis), which is comparable to the Venus orbital distance of 0.72 au.

Conclusions
The multi-planetary system around HD 191939 was previously known to host three transiting sub-Neptunes with very similar radii (HD 191939 b, c, and d) and a non-transiting Saturn-mass planet (HD 191939 e), and also showed evidence for an external long-period planet (HD 191939 f). In this paper, we revisited the system using new RV data from CARMENES and HARPS-N spectrographs in addition to archival data from APF and HIRES. The combined dataset, containing 362 RV measurements spanning over ∼2 yr, allowed the detection of a new non-transiting planetary signal (HD 191939 g) with a period of ∼280 days and a minimum mass of M p sin i = 13.5±2.0 M ⊕ . The planet-to-star distance of HD 191939 g places this new planet in the conservative HZ around the host star. However, our measurements suggest HD 191939 g is likely a gaseous planet.
We also present refined mass and bulk properties for planets HD 191939 b, c, and e. Additionally, we improve the mass determination of HD 191939 d at the 4.6σ level of significance, for which only an upper limit was known.
We determine a mass for HD 191939 d of 2.80±0.60 M ⊕ , leading to a mean bulk density of ρ d = 0.57 g cm −3 . Due to its low density and host-star brightness, HD 191939 d is one of the best puffy targets for atmospheric exploration via transmission spectroscopy. Although the detection of spectral features in puffy atmospheres seems to be challenging, JWST may be capable of detecting the atmosphere of HD 191939 d based on the predicted spectroscopic amplitudes. In particular, our simulations suggest that JWST instruments may break the degeneracy between hazes, clouds, and metallicity effects with a single visit.
With a period of 101 days, HD 191939 e was disconnected from the near resonance chain of the three inner transiting planets (1:3:4). However, the detection of HD 191939 g in a 280-day orbit indicates that these two outer non-transiting planets (e and g) are in a separate relation, close to a 1:3 period resonance. HD 191939 does not seem to be unique in this respect, as there are other multi-planetary systems in the literature where the inner and outer planets are in different resonance chains. Moreover, puffy planets tend to be in resonant orbits, which reinforces the hypothesis of a low mean density for planets e and g.
The singular system architecture of three sub-Neptunes interior to a Saturn-mass planet and a Uranus-mass planet, together with the existence of a very long-period massive companion, makes the HD 191929 system unique. The diversity of planets around this star makes this system a prime target for more follow-up observations.    , and J represent fixed, uniform, normal, beta and Jeffrey's distributions, respectively. The parametrization for (p, b) using (r 1 , r 2 ) (Espinoza 2018) and (q 1 , q 2 ) quadratic limb darkening (Kipping 2013) Fig. 6 but including planet f. Planet f parameters are from Table 3. A40, page 20 of 20