Open Access
Issue
A&A
Volume 669, January 2023
Article Number A40
Number of page(s) 20
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202244120
Published online 04 January 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 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 (RP < 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.

In this work, we combined the previously published RV observations with new time series obtained with the CARMENES and HARPS-N spectrographs, and with the existing TESS photometric data for this target. In particular, the combination of four RV datasets allowed us to refine the physical properties of the system and all the planetary masses. Furthermore, the significant increase in the number of observations and its baseline permitted a check for extra planetary signals beyond the period of planet e. This is also the case for the newly discovered HD 191939 g: a Uranus-mass planet in the habitable zone with an orbital period of ~300 days.

Table 1

Stellar parameters of HD 191939.

2 Observations

2.1 TESS photometry

Listed as TIC 269701147 in the TESS Input Catalog (TIC; Stassun et al. 2018), HD 191939 was observed by TESS in two-minute short-cadence integrations in Sectors 15–19 from 15 August 2019 to 24 December 2019, Sectors 21–22 from 21 January 2020 to 18 March 2020, Sectors 24–25 from 16 April 2020 to 8 June 2020, Sector 41 from 23 July 2021 to 20 August 2021, and Sector 48 from 28 January 2022 to 26 February 2022.

In this work, we made use of the Presearch Data Conditioning-corrected simple aperture photometry (PDC-SAP; Smith et al. 2012; Stumpe et al. 2012, Stumpe et al. 2014; Morris et al. 2017) reduced by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) at the NASA Ames Research Center and publicly available at the Mikulski Archive for Space Telescopes (MAST1).

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.1. We obtained the detrending model by evaluating the GP component at each time point, which includes the transit times. Finally, we divided the PDC-SAP light curve by the detrending model.

2.2 High-resolution spectroscopy follow up

2.2.1 High-resolution spectroscopy with CARMENES

HD 191939 was observed with the Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs (CARMENES; Quirrenbach et al. 2014, 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 ℛ = 94 600, and the near-infrared (NIR) channel, which goes from 0.96 to 1.71 µm with a resolving power of ℛ = 80400. 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 serval2 (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α, NaD1 and NaD2, and CaIIIRT line indices. We also used the RACCOON code3 (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 Fig. 1).

thumbnail Fig. 1

Time series of RV measurements taken by APF (blue circles), HIRES (green up triangles), CARMENES (red squares), and HARPS-N (orange down triangles).

2.2.2 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.6m 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 high-resolution (ℛ =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 tool4, 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 NaD1 and NaD2 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).

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

3 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 derived a Prot/ 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) periodogram5. 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 IRT1 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.

3.2 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 juliet6 (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 etal. 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 ZM2 − ln ZM1 > 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 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 (t0) for the three transiting planets based on a photometric-only analysis. The precision derived for P and t0 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 t0 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 (γ˙$\dot \gamma $) term and a quadratic (γ¨$\ddot \gamma $) 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-Striker7 (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 (GPexponential), Matern 3/2 (GPMatern), and quasi-periodic (GPqp). When fitting the signal with a Keplerian orbit or with a GPqp 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 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 GPqp model are statistically preferred over the four-planet model, the GPexponential model, and the GPMatern model. The ∆ ln Z between the five-planet model and its analogous GPqp model is less than 2. Therefore, the GPqp model is moderately preferred over the five-planet model. However, the GPqp model has one free parameter more than the five-planet model. Moreover, closer inspection of the posterior distributions of the GPqp hyperparameters shows that the Prot is not constrained (see Fig. A.2). In the uninformative GPqp model, the Prot posterior distribution is mainly flat. When we used a normal prior for Prot, 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 GPqp 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.

thumbnail Fig. 2

GLS periodograms of the time series of APF, HIRES, CARMENES, and HARPS-N RV measurements and the residuals after subtraction of different models. All the models include quadratic and linear terms to account for the long-term trend detected in the RV time series. (a) GLS periodogram of RVs after removing the long-term trend. (b) GLS periodogram of the RV residuals after fitting the 101 days signal (vertical magenta line). (c) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 days (vertical red line) and 101 days signals. (d) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 days, 28.6 days (vertical green line), and 101 days signals. (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.

Table 2

Comparative between Bayesian log-evidences (∆ ln Z) and planet semi-amplitudes (Kp) for the different explored models.

3.3 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 parame-terised with the uniform sampling prior (q1,q2) introduced by Kipping (2013). Additionally, rather than directly fitting the impact parameter of the orbit (b) and the planet-to-star radius ratio (p = Rp/R*), we considered the uninformative sample (r1, r2) parameterization introduced in Espinoza (2018). The parameters r1 and r2 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 Zecc − ln Zno 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.4A.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 ≃ (Rp/R*)2) and transit duration tT 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.

Table 3

Parameters and 1σ uncertainties for the juliet joint fit model for HD 191939 planetary system.

3.4 Planet-candidate f constraints

Based on high-angular-resolution imaging, BA20 (Sect. 2.5 therein) did not find companions 5 mag fainter than HD 191939 at 200 mas or 8.4 mag fainter at 1″. However, the RV time series present a long-term trend that L22 cautiously referred to as a planetary object, HD 191939 f, because its mass is most likely below 13 MJ. 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 planet-candidate f (Fig. A.8) also confirm that HD 191939 f properties are not constrained. The period of ≳2200 days is in agreement 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.

thumbnail Fig. 3

Photometry data phase-folded to the period P and central time of transit t0 (shown above each panel, t0 units are BJD – 2 457 000) derived from the joint fit model. Two-minute cadence TESS phase-folded photometry for HD 191939 b (panel a), c (panel b), and d (panel c) along with the best-fit model. Orange points show binned photometry for visualisation. The error bars include the photometric jitter term added in quadrature.

thumbnail Fig. 4

RV data phase-folded to the period P and central time of transit t0 (shown above each panel, t0 units are BJD – 2 457 000) derived from the joint-fit model. Here, we show APF (blue circles), HIRES (green up triangles), CARMENES (red squares), and HARPS-N (orange down triangles) RVs phase-folded for HD 191939 b (panel a), c (panel b), d (panel c), e (panel d), and g (panel e) along with the best-fit model (black line) and the 3σ confidence interval (shaded grey area). The error bars include the instrumental jitter term added in quadrature.

thumbnail Fig. 5

Mass-period (top panel) and mass-radius (bottom panel) diagrams for well-characterised planets with masses and radii measured with a precision better than 30%, from 1 M up to 1.5 MJ and 1 R up to 1.5 RJ, from the TEPCat database (February 2022; Southworth 2011) and http://exoplanet.eu. HD 191939 planets are colour coded and marked with a filled circle with error bars (the radii of HD 191939 e and g are forecasted). Vertical colour bands indicate the ±1σ mass regions of HD 191939 g and e. Temperate planets with Teq = 250–395 K are marked by blue squares. The mass-radius panel also shows theoretical composition models at 300 K from Zeng et al. (2019).

4 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 2848+10$284_{ - 8}^{ + 10}$ 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 semi-amplitude 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 MJ, and 1 R up to 1.5 RJ, along with theoretical composition models of Zeng et al. (2019)8.

4.1 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 (~5–20 M) and well-determined masses, namely: HD 31527 d (P≃274 days, 16±3 M; Mayor et al. 2011), HD 10180 g (P≃602 days, 21±3 M; Lovis et al. 2011), GJ 3138 d (P≃257 days, 10±2 M; Astudillo-Defru et al. 2017), Barnard b (P≃232 days, 3.2±0.4 M; Ribas et al. 2018), GJ 273 days (P~414 days, 11±4 M) and e (P~542 days, 9±4 M; Tuomi et al. 2019), and GJ 687 c (P≃727 days, 16±4 M; Feng et al. 2020). The intermediate- and long-period planet group is also supplemented with some transiting planets discovered by the Kepler space mission (Borucki et al. 2010), as spectroscopic observations confirmed the planetary nature of HIP 41378 d (P = 278 days, <4.6 M), e (P≃369 days, 12±5 M), and f (P=542 days, 12±3 M; Vanderburg et al. 2016; Santerne et al. 2019). However, the faintness of some of the Kepler targets complicates RV follow-up campaigns. Finally, transit-time-variation analyses helped to compute the planetary masses for Kepler-87 c (P = 191 days, 6.4±0.8 M; Ofir et al. 2014), KOI-1783 c (P = 284 days, 15.03.6+4.3$15.0_{ - 3.6}^{ + 4.3}$ M; Vissapragada et al. 2020), and Kepler-90 g (P = 210 days, 15±1 M; Cabrera et al. 2014; Liang et al. 2021).

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 implementation9. The code predicted planetary radii for planets g and e of 3.71.1+1.5$3.7_{ - 1.1}^{ + 1.5}$ R and 12.93.9+4.7$12.9_{ - 3.9}^{ + 4.7}$ R, respectively. From their minimum masses and forecasted radius, the expected mean bulk densities for planets g and e are 1.50.8+2.7$1.5_{ - 0.8}^{ + 2.7}$ g cm−3 and 0.300.16+0.70$0.30_{ - 0.16}^{ + 0.70}$ 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 Teq 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-90g 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.

thumbnail Fig. 6

Diagram of the planetary system HD 191939. Planetary orbits are colour coded and planets (filled circles) are scaled according to their mass. The shaded area marks the region of Teq = 273–373 K. HD 191939 is marked at the centre with a black star. A diagram including planet f is shown in Fig. A.9.

4.2 HD 191939 d: a puffy planet

L22 already noted that HD 191939 d was probably a low-density 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 (Teq ≃ 540 K), this explanation is unlikely.

The nearest planet to HD 191939 d in the mass-radius diagram is Kepler-79 e (KOI-152e; Jontof-Hutter et al. 2014), which has very similar properties (3.49±0.14 R, ρp =0.53±0.15 g cm−3, Teq ≃ 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 (TSMb = 153; TSMc = 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-90g (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 Teq ~ 500–600 K in the database of exoplanet atmospheric observations of ExoAtmospheres10. Only the warm sub-Neptune GJ 1214 b (2.74±0.05 R, Teq ≃ 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 1214b (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 H2—H2 and H2—He (Abel et al. 2011, 2012; Fletcher 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, low-density 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 modes at median resolving power of ℛ ~ 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 191939d are seven times smallerthan those reported forHIP41378 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 H2O and CH4, which are an order of magnitude larger than the predicted error bars. Table 4 reports the amplitudes of spectral modulations at low resolution (ℛ ~ 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.

thumbnail Fig. 7

Synthetic transmission spectra for HD 191939 d. Models assuming a clear atmosphere with solar abundances (solid red line), a clear atmosphere with metallicity enhanced by a factor of 100 (blue dashed line), and a hazy atmosphere with solar abundances (green dotted line). Simulated measurements with error bars are shown for the observation of one transit with JWST NIRISS-SOSS, NIRSpec-G395M, and MIRI-LRS configurations.

Table 4

Transmission spectroscopy amplitudes at low resolution (ℛ ~ 170) for various models of the HD 191939 d atmosphere and nominal wavelength ranges of HST and JWST instrumental modes.

4.3 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 (Pb = 8.88 days, Pc = 28.58 days, Pd = 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 non-transiting planets of the system appear to be in a period ratio of 1:3 (Pe = 101 days, Pg ≃ 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 HIP41378 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 MJ 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 Archive11 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 (3528+18$35_{ - 28}^{ + 18}$ M, 6.6±0.6 R), but without evidence for a long-period companion like HD 191939f. 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-94d 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 (Pb = 3.7 days, Pc = 10.4 days, Pd = 22.3 days), with the giant planet also close to 5:2 with the outermost planet KOI-94 e (Pe = 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 semi-major 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 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.

thumbnail Fig. 8

Planet mass distribution across planets for HD 191939 (blue circles), KOI-94 (orange squares), Kepler-65 (green down triangles), Kepler-68 (red up triangles), and the Solar System (purple stars) systems. Masses are scaled to HD 191939 e minimum mass, and semi-major axes are scaled to that of the massive planet for each system.

5 Conclusions

The multi-planetary system around HD 191939 was previously known to host three transiting sub-Neptunes with very similar radii (HD 191939b, 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 Mp 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.

Acknowledgements

This work was supported by the KESPRINT collaboration, an international consortium devoted to the characterization and research of exoplanets discovered with space-based missions (www.kesprint.science). CARMENES is an instrument at the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the Deutsche Forschungsgemeinschaft (DFG) through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. This is University of Texas Center for Planetary Systems Habitability contribution number 0058. This research has made use of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia. G.M. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie grant agreement no. 895525. J.K. gratefully acknowledges the support of the Swedish National Space Agency (SNSA; DNR 2020-00104). P.K. acknowledges the support from the grant LTT-20015. K.W.F.L. was supported by Deutsche Forschungsgemeinschaft grants RA714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. C.M.P. gratefully acknowledges the support of the Swedish National Space Agency (DNR 65/19). The first author acknowledges the special support received by P. Conxa, P. Mercè, Jeroni, and Mercè. This work has made use of resources from AstroPiso collaboration. J.O.M. gratefully acknowledges the inspiring discussions with Yess, Alejandro, his colleagues, and friends.

Appendix A Additional figures and tables

Table A.1

Prior and posterior distributions for PDC-SAP detrending fit. Prior labels ℱ, 𝒩, and 𝒥 represent fixed, normal, and Jeffrey’s distributions, respectively.

Table A.2

Prior and posterior distributions for the RV GP models explored in Sect. 3.2. Prior labels 𝒰, 𝒩, and 𝒥 represent uniform, normal, and Jeffrey's distributions, respectively. Corner plots for the quasi-periodic model hyperparameters are shown in Fig. A.2.

thumbnail Fig. A.1

Generalised Lomb-Scargle periodograms of the activity indices from CARMENES (left) and HARPS-N (right), and S-index (bottom right) from APF, HIRES, and HARPS-N. In all panels, the broken vertical lines indicate the planetary signals at 8.9 (red), 28.6 (green), 38.4 (cyan), 101 (magenta), and 280 (orange) days. In all panels, the 10%, 1%, and 0.1% FAP levels are indicated by dotted, dash-dotted, and dashed grey lines, respectively. The vertical black dotted line marks the baseline for each dataset. We highlight the different scale in the y axis in each panel.

thumbnail Fig. A.2

Posterior distribution of hyperparameters from the uninformative (left) and constrained (right) GPqp models used in Sect. 3.2. Prior distributions are in Table A.2.

Table A.3

Priors for each parameter used in the juliet joint fit model for HD 191939 planetary system. Prior labels ℱ, 𝒰, 𝒩, ℬ, and 𝒥 represent fixed, uniform, normal, beta and Jeffrey's distributions, respectively. The parametrization for (p, b) using (r1, r2) (Espinoza 2018) and (q1, q2) quadratic limb darkening (Kipping 2013) are both explained in Section 3.3.

thumbnail Fig. A.3

Reproduction of the RV analyses presented in Fig. 1 from L22. All the models include quadratic and linear terms to account for the long-term trend. (a) GLS periodogram of APF and HIRES datasets. (b) GLS periodogram of the RV residuals after fitting the 101 d signal (vertical magenta line). (c) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d (vertical red line) and 101 d signals. (d) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d, 28.6 d (vertical green line), and 101 d signals. (e) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d, 28.6 d, 38 d (vertical cyan line), and 101 d 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 in the y axis in each panel.

thumbnail Fig. A.4

TESS photometry from Sectors 15–19, 21, 22, 24, and 25 along with the best-fit model (see Fig. A.5 for Sectors 41 and 48). Upward-pointing triangles mark the transits for HD 191939 b (red), c (cyan), and d (green). Downward-pointing triangles with error bars mark the expected t0 and ± 1σ uncertainty for the non-transiting planets HD 191939 e (magenta) and g (orange).

thumbnail Fig. A.5

As in Fig. A.4 but for Sectors 41 and 48.

thumbnail Fig. A.6

RV results from the joint fit model. Detrended RV time series (top panel) of APF (blue circles), HIRES (green up triangles), CARMENES (red squares), and HARPS-N (orange down triangles) along with the best-fit Keplerian model (black line) and the 3σ confidence interval (shaded grey area). The error bars include the instrumental jitter term added in quadrature.

thumbnail Fig. A.7

As in Fig. A.6 but considering a Keplerian signal for planet f (see Sect. 3.4).

thumbnail Fig. A.8

RVs phase-folded to the P and t0 (shown above the panel, P units are days and t0 units are BJD-2 457 000) for planet f along with the best-fit model (black line) and the 1σ, 2σ, and 3σ confidence intervals (shaded grey areas).

thumbnail Fig. A.9

As in Fig. 6 but including planet f. Planet f parameters are from Table 3.

References

  1. Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2011, J. Phys. Chem. A, 115, 6805 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2012, J. Chem. Phys., 136, 044319 [NASA ADS] [CrossRef] [Google Scholar]
  3. Agúndez, M., Venot, O., Iro, N., et al. 2012, A&A, 548, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 917, 37 [NASA ADS] [CrossRef] [Google Scholar]
  5. Alam, M. K., Kirk, J., Dressing, C. D., et al. 2022, ApJ, 927, L5 [NASA ADS] [CrossRef] [Google Scholar]
  6. Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252 [Google Scholar]
  7. Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017, A&A, 602, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Badenas-Agusti, M., Günther, M. N., Daylan, T., et al. 2020, AJ, 160, 113 [NASA ADS] [CrossRef] [Google Scholar]
  9. Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bean, J. L., Miller-Ricci Kempton, E., & Homeier, D. 2010, Nature, 468, 669 [NASA ADS] [CrossRef] [Google Scholar]
  11. Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [Google Scholar]
  12. Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Caballero, J. A., Guàrdia, J., López del Fresno, M., et al. 2016, SPIE Conf. Ser., 9910, 99100E [Google Scholar]
  14. Cabrera, J., Csizmadia, S., Lehmann, H., et al. 2014, ApJ, 781, 18 [Google Scholar]
  15. Cannon, A. J., & Pickering, E. C. 1993, VizieR Online Data Catalog: III/135A [Google Scholar]
  16. Chachan, Y., Jontof-Hutter, D., Knutson, H. A., et al. 2020, AJ, 160, 201 [NASA ADS] [CrossRef] [Google Scholar]
  17. Chen, J., & Kipping, D. 2017, ApJ, 834, 17 [Google Scholar]
  18. Cloutier, R., Charbonneau, D., Deming, D., Bonfils, X., & Astudillo-Defru, N. 2021, AJ, 162, 174 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cosentino, R., Lovis, C., Pepe, F., et al. 2012, SPIE Conf. Ser., 8446, 84461V [Google Scholar]
  20. Cosentino, R., Lovis, C., Pepe, F., et al. 2014, SPIE Conf. Ser., 9147, 91478C [Google Scholar]
  21. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
  22. Dragomir, D., Teske, J., Günther, M. N., et al. 2019, ApJ, 875, L7 [Google Scholar]
  23. Espinoza, N. 2018, RNAAS, 2, 209 [Google Scholar]
  24. Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262 [Google Scholar]
  25. Espinoza, N., Pallé, E., Kemmer, J., et al. 2022, AJ, 163, 133 [NASA ADS] [CrossRef] [Google Scholar]
  26. Feng, F., Shectman, S. A., Clement, M. S., et al. 2020, ApJS, 250, 29 [Google Scholar]
  27. Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  28. Fletcher, L. N., Gustafsson, M., & Orton, G. S. 2018, ApJS, 235, 24 [NASA ADS] [CrossRef] [Google Scholar]
  29. Foreman-Mackey, D., Agol, E., Angus, R., & Ambikasaran, S. 2017, AJ, 154, 220 [NASA ADS] [CrossRef] [Google Scholar]
  30. Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80 [Google Scholar]
  31. Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [Google Scholar]
  32. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Gao, P., & Zhang, X. 2020, ApJ, 890, 93 [Google Scholar]
  34. Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39 [NASA ADS] [CrossRef] [Google Scholar]
  35. Günther, M. N., Pozuelos, F. J., Dittmann, J. A., et al. 2019, Nat. Astron., 3, 1099 [Google Scholar]
  36. Hadden, S., & Lithwick, Y. 2017, AJ, 154, 5 [Google Scholar]
  37. Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27 [Google Scholar]
  38. Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, SPIE Conf. Ser., 9913, 99133E [Google Scholar]
  39. Jontof-Hutter, D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, ApJ, 785, 15 [NASA ADS] [CrossRef] [Google Scholar]
  40. Kasper, D., Bean, J. L., Oklopčić, A., et al. 2020, AJ, 160, 258 [Google Scholar]
  41. Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401 [Google Scholar]
  42. Kipping, D. M. 2013, MNRAS, 435, 2152 [Google Scholar]
  43. Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
  44. Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69 [Google Scholar]
  45. Lafarga, M., Ribas, I., Lovis, C., et al. 2020, A&A, 636, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Lee, J.-M., Heng, K., & Irwin, P. G. J. 2013, ApJ, 778, 97 [Google Scholar]
  47. Liang, Y., Robnik, J., & Seljak, U. 2021, AJ, 161, 202 [NASA ADS] [CrossRef] [Google Scholar]
  48. Libby-Roberts, J. E., Berta-Thompson, Z. K., Désert, J.-M., et al. 2020, AJ, 159, 57 [Google Scholar]
  49. Lovis, C., Ségransan, D., Mayor, M., et al. 2011, A&A, 528, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Lubin, J., Van Zandt, J., Holcomb, R., et al. 2022, AJ, 163, 101 [NASA ADS] [CrossRef] [Google Scholar]
  51. Luque, R., Serrano, L. M., Molaverdikhani, K., et al. 2021, A&A, 645, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Luque, R., Fulton, B. J., Kunimoto, M., et al. 2022, A&A, 664, A199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Masuda, K. 2014, ApJ, 783, 53 [Google Scholar]
  54. Mayor, M., Marmier, M., Lovis, C., et al. 2011, A&A, submitted, [arXiv: 1109.2497] [Google Scholar]
  55. Mills, S. M., Howard, A. W., Weiss, L. M., et al. 2019, AJ, 157, 145 [NASA ADS] [CrossRef] [Google Scholar]
  56. Morello, G., Zingales, T., Martin-Lagarde, M., Gastaud, R., & Lagage, P.-O. 2021, AJ, 161, 174 [NASA ADS] [CrossRef] [Google Scholar]
  57. Morris, R. L., Twicken, J. D., Smith, J. C., et al. 2017, Kepler Data Processing Handbook: Photometric Analysis, Kepler Science Document KSCI-19081-002 [Google Scholar]
  58. Ofir, A., Dreizler, S., Zechmeister, M., & Husser, T.-O. 2014, A&A, 561, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Ohno, K., & Tanaka, Y. A. 2021, ApJ, 920, 124 [NASA ADS] [CrossRef] [Google Scholar]
  60. Orell-Miquel, J., Murgas, F., Pallé, E., et al. 2022, A&A, 659, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Petit dit de la Roche, D. J. M., van den Ancker, M. E., & Miles-Paez, P. A. 2020, RNAAS, 4, 231 [NASA ADS] [Google Scholar]
  62. Quinn, S. N., Becker, J. C., Rodriguez, J. E., et al. 2019, AJ, 158, 177 [Google Scholar]
  63. Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, SPIE Conf. Ser., 9147, 91471F [Google Scholar]
  64. Quirrenbach, A., CARMENES Consortium, Amado, P. J., et al. 2020, SPIE Conf. Ser., 11447, 114473C [Google Scholar]
  65. Ribas, I., Tuomi, M., Reiners, A., et al. 2018, Nature, 563, 365 [Google Scholar]
  66. Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [Google Scholar]
  67. Santerne, A., Malavolta, L., Kosiarek, M. R., et al. 2019, arXiv e-prints [arXiv:1911.07355] [Google Scholar]
  68. Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
  69. Soubiran, C., Jasniewicz, G., Chemin, L., et al. 2018, A&A, 616, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Southworth, J. 2011, MNRAS, 417, 2166 [Google Scholar]
  71. Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
  72. Stassun, K. G., & Torres, G. 2018, ApJ, 862, 61 [Google Scholar]
  73. Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102 [Google Scholar]
  74. Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985 [Google Scholar]
  75. Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
  76. Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64 [NASA ADS] [CrossRef] [Google Scholar]
  77. Trifonov, T. 2019, Astrophysics Source Code Library [record ascl:1906.004] [Google Scholar]
  78. Trifonov, T., Tal-Or, L., Zechmeister, M., et al. 2020, A&A, 636, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Trotta, R. 2008, Contemp. Phys., 49, 71 [Google Scholar]
  80. Tsiaras, A., Waldmann, I. P., Zingales, T., et al. 2018, AJ, 155, 156 [NASA ADS] [CrossRef] [Google Scholar]
  81. Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, arXiv e-prints [arXiv: 1906.04644] [Google Scholar]
  82. Van Eylen, V., & Albrecht, S. 2015, ApJ, 808, 126 [Google Scholar]
  83. Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, AJ, 157, 61 [Google Scholar]
  84. van Leeuwen, F. 2007, A&A, 474, 653 [CrossRef] [EDP Sciences] [Google Scholar]
  85. Vanderburg, A., Becker, J. C., Kristiansen, M. H., et al. 2016, ApJ, 827, L10 [Google Scholar]
  86. Vissapragada, S., Jontof-Hutter, D., Shporer, A., et al. 2020, AJ, 159, 108 [NASA ADS] [CrossRef] [Google Scholar]
  87. Vogt, S. S., Radovan, M., Kibrick, R., et al. 2014, PASP, 126, 359 [NASA ADS] [CrossRef] [Google Scholar]
  88. Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [Google Scholar]
  89. Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19 [Google Scholar]
  90. Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proc. Natl. Acad. Sci. USA, 113, 11431 [NASA ADS] [CrossRef] [Google Scholar]
  91. Yu, X., He, C., Zhang, X., et al. 2021, Nat. Astron., 5, 822 [NASA ADS] [CrossRef] [Google Scholar]
  92. Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [CrossRef] [EDP Sciences] [Google Scholar]
  93. Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proc. Natl. Acad. Sci. USA, 116, 9723 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Stellar parameters of HD 191939.

Table 2

Comparative between Bayesian log-evidences (∆ ln Z) and planet semi-amplitudes (Kp) for the different explored models.

Table 3

Parameters and 1σ uncertainties for the juliet joint fit model for HD 191939 planetary system.

Table 4

Transmission spectroscopy amplitudes at low resolution (ℛ ~ 170) for various models of the HD 191939 d atmosphere and nominal wavelength ranges of HST and JWST instrumental modes.

Table A.1

Prior and posterior distributions for PDC-SAP detrending fit. Prior labels ℱ, 𝒩, and 𝒥 represent fixed, normal, and Jeffrey’s distributions, respectively.

Table A.2

Prior and posterior distributions for the RV GP models explored in Sect. 3.2. Prior labels 𝒰, 𝒩, and 𝒥 represent uniform, normal, and Jeffrey's distributions, respectively. Corner plots for the quasi-periodic model hyperparameters are shown in Fig. A.2.

Table A.3

Priors for each parameter used in the juliet joint fit model for HD 191939 planetary system. Prior labels ℱ, 𝒰, 𝒩, ℬ, and 𝒥 represent fixed, uniform, normal, beta and Jeffrey's distributions, respectively. The parametrization for (p, b) using (r1, r2) (Espinoza 2018) and (q1, q2) quadratic limb darkening (Kipping 2013) are both explained in Section 3.3.

All Figures

thumbnail Fig. 1

Time series of RV measurements taken by APF (blue circles), HIRES (green up triangles), CARMENES (red squares), and HARPS-N (orange down triangles).

In the text
thumbnail Fig. 2

GLS periodograms of the time series of APF, HIRES, CARMENES, and HARPS-N RV measurements and the residuals after subtraction of different models. All the models include quadratic and linear terms to account for the long-term trend detected in the RV time series. (a) GLS periodogram of RVs after removing the long-term trend. (b) GLS periodogram of the RV residuals after fitting the 101 days signal (vertical magenta line). (c) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 days (vertical red line) and 101 days signals. (d) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 days, 28.6 days (vertical green line), and 101 days signals. (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.

In the text
thumbnail Fig. 3

Photometry data phase-folded to the period P and central time of transit t0 (shown above each panel, t0 units are BJD – 2 457 000) derived from the joint fit model. Two-minute cadence TESS phase-folded photometry for HD 191939 b (panel a), c (panel b), and d (panel c) along with the best-fit model. Orange points show binned photometry for visualisation. The error bars include the photometric jitter term added in quadrature.

In the text
thumbnail Fig. 4

RV data phase-folded to the period P and central time of transit t0 (shown above each panel, t0 units are BJD – 2 457 000) derived from the joint-fit model. Here, we show APF (blue circles), HIRES (green up triangles), CARMENES (red squares), and HARPS-N (orange down triangles) RVs phase-folded for HD 191939 b (panel a), c (panel b), d (panel c), e (panel d), and g (panel e) along with the best-fit model (black line) and the 3σ confidence interval (shaded grey area). The error bars include the instrumental jitter term added in quadrature.

In the text
thumbnail Fig. 5

Mass-period (top panel) and mass-radius (bottom panel) diagrams for well-characterised planets with masses and radii measured with a precision better than 30%, from 1 M up to 1.5 MJ and 1 R up to 1.5 RJ, from the TEPCat database (February 2022; Southworth 2011) and http://exoplanet.eu. HD 191939 planets are colour coded and marked with a filled circle with error bars (the radii of HD 191939 e and g are forecasted). Vertical colour bands indicate the ±1σ mass regions of HD 191939 g and e. Temperate planets with Teq = 250–395 K are marked by blue squares. The mass-radius panel also shows theoretical composition models at 300 K from Zeng et al. (2019).

In the text
thumbnail Fig. 6

Diagram of the planetary system HD 191939. Planetary orbits are colour coded and planets (filled circles) are scaled according to their mass. The shaded area marks the region of Teq = 273–373 K. HD 191939 is marked at the centre with a black star. A diagram including planet f is shown in Fig. A.9.

In the text
thumbnail Fig. 7

Synthetic transmission spectra for HD 191939 d. Models assuming a clear atmosphere with solar abundances (solid red line), a clear atmosphere with metallicity enhanced by a factor of 100 (blue dashed line), and a hazy atmosphere with solar abundances (green dotted line). Simulated measurements with error bars are shown for the observation of one transit with JWST NIRISS-SOSS, NIRSpec-G395M, and MIRI-LRS configurations.

In the text
thumbnail Fig. 8

Planet mass distribution across planets for HD 191939 (blue circles), KOI-94 (orange squares), Kepler-65 (green down triangles), Kepler-68 (red up triangles), and the Solar System (purple stars) systems. Masses are scaled to HD 191939 e minimum mass, and semi-major axes are scaled to that of the massive planet for each system.

In the text
thumbnail Fig. A.1

Generalised Lomb-Scargle periodograms of the activity indices from CARMENES (left) and HARPS-N (right), and S-index (bottom right) from APF, HIRES, and HARPS-N. In all panels, the broken vertical lines indicate the planetary signals at 8.9 (red), 28.6 (green), 38.4 (cyan), 101 (magenta), and 280 (orange) days. In all panels, the 10%, 1%, and 0.1% FAP levels are indicated by dotted, dash-dotted, and dashed grey lines, respectively. The vertical black dotted line marks the baseline for each dataset. We highlight the different scale in the y axis in each panel.

In the text
thumbnail Fig. A.2

Posterior distribution of hyperparameters from the uninformative (left) and constrained (right) GPqp models used in Sect. 3.2. Prior distributions are in Table A.2.

In the text
thumbnail Fig. A.3

Reproduction of the RV analyses presented in Fig. 1 from L22. All the models include quadratic and linear terms to account for the long-term trend. (a) GLS periodogram of APF and HIRES datasets. (b) GLS periodogram of the RV residuals after fitting the 101 d signal (vertical magenta line). (c) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d (vertical red line) and 101 d signals. (d) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d, 28.6 d (vertical green line), and 101 d signals. (e) GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d, 28.6 d, 38 d (vertical cyan line), and 101 d 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 in the y axis in each panel.

In the text
thumbnail Fig. A.4

TESS photometry from Sectors 15–19, 21, 22, 24, and 25 along with the best-fit model (see Fig. A.5 for Sectors 41 and 48). Upward-pointing triangles mark the transits for HD 191939 b (red), c (cyan), and d (green). Downward-pointing triangles with error bars mark the expected t0 and ± 1σ uncertainty for the non-transiting planets HD 191939 e (magenta) and g (orange).

In the text
thumbnail Fig. A.5

As in Fig. A.4 but for Sectors 41 and 48.

In the text
thumbnail Fig. A.6

RV results from the joint fit model. Detrended RV time series (top panel) of APF (blue circles), HIRES (green up triangles), CARMENES (red squares), and HARPS-N (orange down triangles) along with the best-fit Keplerian model (black line) and the 3σ confidence interval (shaded grey area). The error bars include the instrumental jitter term added in quadrature.

In the text
thumbnail Fig. A.7

As in Fig. A.6 but considering a Keplerian signal for planet f (see Sect. 3.4).

In the text
thumbnail Fig. A.8

RVs phase-folded to the P and t0 (shown above the panel, P units are days and t0 units are BJD-2 457 000) for planet f along with the best-fit model (black line) and the 1σ, 2σ, and 3σ confidence intervals (shaded grey areas).

In the text
thumbnail Fig. A.9

As in Fig. 6 but including planet f. Planet f parameters are from Table 3.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.