11institutetext:  Instituto de Astrofísica de Canarias (IAC), 38205 La Laguna,
Tenerife, Spain
11email<EMAIL_ADDRESS>22institutetext: Departamento de Astrofísica, Universidad
de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain. 33institutetext:
Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía
s/n, 18008 Granada, Spain. 44institutetext: Department of Astronomy &
Astrophysics, University of Chicago, Chicago, IL 60637, USA. 55institutetext:
Department of Earth, Atmospheric and Planetary Sciences, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA. 66institutetext: Department
of Physics, and Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA.
77institutetext: Institut de Ciències de l’Espai (CSIC), Campus UAB, c/ de Can
Magrans s/n, E-08193 Bellaterra, Barcelona, Spain. 88institutetext: Institut
d’Estudis Espacials de Catalunya, E-08034 Barcelona, Spain. 99institutetext:
Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4
7AL, United Kingdom. 1010institutetext: Max-Planck-Institut für Astronomie,
Königstuhl 17, 69117 Heidelberg, Germany. 1111institutetext: Institut für
Astrophysik, Georg-August-Universität, Friedrich-Hund-Platz 1, 37077
Göttingen, Germany. 1212institutetext: McDonald Observatory and Center for
Planetary Systems Habitability. 1313institutetext: The University of Texas,
Austin, Texas, USA. 1414institutetext: Dipartimento di Fisica, Universitá di
Torino, via P. Giuria 1, I-10125 Torino, Italy. 1515institutetext: Thüringer
Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany.
1616institutetext: Astronomical Institute of the Czech Academy of Sciences,
Fričova 298, 25165, Ondřejov, Czech Republic. 1717institutetext: Department of
Space, Earth and Environment, Astronomy and Plasma Physics, Chalmers
University of Technology, Chalmersplatsen 4, 412 96 Gothenburg, Sweden.
1818institutetext: Institute of Planetary Research, German Aerospace Center
(DLR), Rutherfordstrasse 2, D-12489 Berlin, Germany. 1919institutetext:
Astrobiology Center, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan.
2020institutetext: National Astronomical Observatory of Japan, 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan. 2121institutetext: Department of Astronomy, The
Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka,
Tokyo, Japan. 2222institutetext: Department of Space, Earth and Environment,
Astronomy and Plasma Physics, Chalmers University of Technology, 412 96
Gothenburg, Sweden. 2323institutetext: Department of Space, Earth and
Environment, Chalmers University of Technology, Onsala Space Observatory,
SE-439 92 Onsala, Sweden. 2424institutetext: Mullard Space Science Laboratory,
University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK.

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

J. Orell-Miquel HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone.HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone. G. Nowak HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone. F. Murgas HD 191939 revisited: New and
refined planet mass determinations, and a new planet in the habitable zone.HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone. E. Palle HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone.HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone. G. Morello HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone. R. Luque HD 191939 revisited: New and
refined planet mass determinations, and a new planet in the habitable zone.HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone. M. Badenas-Agusti HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone.HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone. I. Ribas HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone. M. Lafarga HD 191939 revisited: New and
refined planet mass determinations, and a new planet in the habitable zone.HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone. N. Espinoza HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone.HD 191939
revisited: New and refined planet mass determinations, and a new planet in the
habitable zone. J. C. Morales HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone.HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone. M. Zechmeister HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone. A. Alqasim HD 191939
revisited: New and refined planet mass determinations, and a new planet in the
habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone. W. D. Cochran HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone. D. Gandolfi HD 191939 revisited: New and
refined planet mass determinations, and a new planet in the habitable zone.HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone. E. Goffo HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone.HD 191939
revisited: New and refined planet mass determinations, and a new planet in the
habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone. P. Kabáth HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone. J. Korth HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone. K. W. F. Lam HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone. J. Livingston HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.HD 191939 revisited:
New and refined planet mass determinations, and a new planet in the habitable
zone.HD 191939 revisited: New and refined planet mass determinations, and a
new planet in the habitable zone.HD 191939 revisited: New and refined planet
mass determinations, and a new planet in the habitable zone.HD 191939
revisited: New and refined planet mass determinations, and a new planet in the
habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone. A. Muresan HD 191939
revisited: New and refined planet mass determinations, and a new planet in the
habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone. C. M. Persson HD
191939 revisited: New and refined planet mass determinations, and a new planet
in the habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone. V. Van Eylen HD 191939
revisited: New and refined planet mass determinations, and a new planet in the
habitable zone.HD 191939 revisited: New and refined planet mass
determinations, and a new planet in the habitable zone.

(Received 27 May 2022 / Accepted 16 October 2022)

HD 191939 (TOI-1339) is a nearby (d $=$ 54 pc), bright (V $=$ 9 mag), 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.

###### Key Words.:

stars: individual: HD 191939 – planets and satellites: individual: HD 191939 g
– planets and satellites: individual: HD 191939 d – planetary systems –
techniques: photometric – techniques: radial velocities

##  1 Introduction

Table 1:  Stellar parameters of HD 191939.

Parameter | Value | Reference
---|---|---
Name | HD 191939 | HD
| TIC 269701147 | TESS
| TOI-1339 | TOI
| HIP 99175 | HIP
Coordinates and spectral type
$\alpha$ (J2000) | 20h 08m 05s.755 | Gaia
$\delta$ (J2000) |  $+$66º 51’ 02”.077 | Gaia
Spectral type | G9 V | L22
Parallax and kinematics
$\mu_{\alpha}$ [mas yr-1] | 150.256 $\pm$ 0.044 | Gaia
$\mu_{\delta}$ [mas yr-1] |  $-$63.909 $\pm$ 0.047 | Gaia
Parallax [mas] | 18.706 $\pm$ 0.071 | Gaia
Distance [pc] | 53.48$\pm$ $-$0.20 | BA20
$V_{\mathrm{r}}$ [$\mathrm{km\,s^{-1}}$] |  $-$9.266 $\pm$ 0.002 | Gaia
Magnitudes
$B$ [mag] | 9.720 $\pm$ 0.038 | TYC
$V$ [mag] | 8.97 $\pm$ 0.03 | HIP
$G$ [mag] | 8.7748 $\pm$ 0.0002 | Gaia
$T$ [mag] | 8.292 $\pm$ 0.006 | TESS
$J$ [mag] | 7.597 $\pm$ 0.029 | 2MASS
$H$ [mag] | 7.215 $\pm$ 0.023 | 2MASS
$K$ [mag] | 7.180 $\pm$ 0.021 | 2MASS
Stellar parameters
Radius $R_{\star}$ [$R_{\odot}$] | 0.94 $\pm$ 0.02 | BA20, L22
Mass $M_{\star}$ [$M_{\odot}$] | 0.92 $\pm$ 0.06 | BA20
| 0.81 $\pm$ 0.04 | L22
$\rho_{\star}$ [g cm-3] | 1.56 $\pm$ 0.15 | BA20
| 1.37 $\pm$ 0.11 | L22
$L_{\star}$ [$L_{\odot}$] | 0.69 $\pm$ 0.01 | BA20
| 0.65 $\pm$ 0.02 | L22
$T_{\mathrm{eff}}$ [K] | 5427 $\pm$ 50 | BA20
| 5348 $\pm$ 100 | L22
$\log(g\,[$cm s${}^{-2}])$ | 4.3 $\pm$ 0.1 | L22
Metallicity [Fe/H] |  $-$0.15 $\pm$ 0.06 | L22
Age [Gyr] | 7 $\pm$ 3 | BA20
$\log($R${}^{{}^{\prime}}_{HK})$ |  $-$5.11 $\pm$ 0.05 | L22
$v\,\sin{i}$ [km s-1] | 0.6 $\pm$ 0.5 | BA20
|  $<$2.0 | L22

HD: Cannon & Pickering (1993); TESS: Stassun et al. (2018); TOI: Guerrero et
al. (2021); HIP: van Leeuwen (2007); Gaia: Gaia Collaboration et al. (2018),
Soubiran et al. (2018); L22: Lubin et al. (2022); BA20: Badenas-Agusti et al.
(2020); TYC: Høg et al. (2000); 2MASS: Cutri et al. (2003). Parallax is
corrected for a systematic offset of $+$0.082$\pm$0.033 mas, as described in
Stassun & Torres (2018).

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al., 2015) is a
NASA-sponsored space telescope whose original mission was a two-year full-sky
survey searching for transiting planets. One of TESS main scientific goals is
to look for small planets ($R_{\mathrm{P}}<4\,R_{\oplus}$) 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 d, respectively. Their radii are very similar, ranging
3.16–3.42 $R_{\oplus}$. 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$\pm$0.9 $M_{\oplus}$ and 7.2$\pm$1.4
$M_{\oplus}$ for planets b and c, respectively. L22 only presented an upper
limit of 5.8 $M_{\oplus}$ ($2\sigma$ 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 d and a minimum mass of 108$\pm$3
$M_{\oplus}$. 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 d. 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 $\sim$300 d.

##  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, 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
(MAST111https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html).

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

Figure 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 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 $\mu$m with a resolving
power of $\mathcal{R}$ = 94 600, and the near-infrared (NIR) channel, which
goes from 0.96 to 1.71 $\mu$m with a resolving power of $\mathcal{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
serval222https://github.com/mzechmeister/serval (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$\alpha$, Na D1 and Na D2, and Ca II IRT
line indices. We also used the RACCOON
code333https://github.com/mlafarga/raccoon (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 Figure 1).

####  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 ($\mathcal{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 seconds, 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
tool444Available at http://ia2-harps.oats.inaf.it:8000., 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$\alpha$ 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 $\mathrm{m\,s^{-1}}$, with a
mean value of 1.07 $\mathrm{m\,s^{-1}}$. In the analysis presented here, we
used the relative RVs measured from HARPS-N spectra with serval (shown in
Figure 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 Figure 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
$P_{\mathrm{rot}}/\sin{i}$ = 79$\pm$66 d, 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)
periodogram555https://github.com/mzechmeister/GLS. 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 Figure 9.
None of the indices present a significant peak at either the periodicities of
the known planets or at the $\sim$300 d 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 d
that is not detected in the HARPS-N indices, and CARMENES Ca II IRT1 and
HARPS-N H$\alpha$ display significant periodicities near 100 d and 200 d,
respectively. However, none of those peaks have a counterpart in the RVs
periodogram analyses.

###  3.2 Radial velocity analysis

Figure 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 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. $(f)$
GLS periodogram of the RV residuals after simultaneously fitting the 8.8 d,
28.6 d, 38 d, 101 d, and 300 d (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.

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 d) and f (¿1700 d).

We analysed the planetary signals in the RVs computing the GLS periodogram and
modelling the detected signals with
juliet666https://juliet.readthedocs.io/en/latest/index.html (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., 2014; 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
$\Delta\ln{Z}=\ln{Z_{\mathrm{M2}}}-\ln{Z_{\mathrm{M1}}}>5$, then the M2 model
is strongly preferred statistically over the M1 model (Trotta, 2008). If
$\Delta\ln{Z}\leq 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
($\Delta\ln{Z}\sim 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 photometric-only 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 ($\sigma$) and a systemic velocity
($\gamma$) term. The procedure described below to fit the periodicities
detected in the periodograms is illustrated in Figure 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 ($\Delta\ln{Z}>5$) over an
only linear or only quadratic term model. The conspicuous signal at 101 d 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 $\sim$300 d (FAP $<$ 1%). As in Fig. 1 from L22, the signature of planet d
at 38.35 d 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 $\sim$300 d which slightly
increased its significance until FAP $\sim$ 0.1%. The $\sim$300 d peak is
already visible and significant after removing planet e.

(f) The periodogram after fitting the $\sim$300 d signal is flat without
significant peaks. We refer to the $\sim$300 d signal as HD 191939 g. We
analyse this signal in detail below.

Table 2:  Comparative between Bayesian log-evidences ($\Delta\ln{Z}$) and
planet semi-amplitudes ($K_{\mathrm{p}}$) for the different explored models.
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_{\mathrm{p}}$ from the joint fit in Sect. 3.3 for illustrative purposes.

Model | $\Delta\ln{Z}$ |  $K_{\mathrm{b}}$ ($P=8.8$ d) |  $K_{\mathrm{c}}$ ($P=28.6$ d) |  $K_{\mathrm{d}}$ ($P=38$ d) |  $K_{\mathrm{e}}$ ($P=101$ d) |  $K_{\mathrm{g}}$ ($P\simeq 300$ d)
---|---|---|---|---|---|---
4 planets | 0.0 | 3.66$\pm$0.25 | 1.75$\pm$0.25 | 0.57$\pm$0.24 | 17.7$\pm$0.3 | –
GPexponential | 3.1 | 3.62$\pm$0.30 | 1.76$\pm$0.30 | 0.58$\pm$0.28 | 17.8$\pm$0.3 | –
GPMatern | 4.2 | 3.58$\pm$0.25 | 1.93$\pm$0.32 | 0.67$\pm$0.33 | 17.8$\pm$0.5 | –
GPqp (a) | 9.2 | 3.55$\pm$0.26 | 1.96$\pm$0.34 | 0.70$\pm$0.35 | 17.8$\pm$0.5 | –
5 planets (a) | 7.7 | 3.56$\pm$0.25 | 1.91$\pm$0.26 | 0.61$\pm$0.25 | 17.6$\pm$0.3 | 1.3$\pm$0.3
GPqp (b) | 10.0 | 3.54$\pm$0.26 | 1.96$\pm$0.33 | 0.72$\pm$0.35 | 17.8$\pm$0.5 | –
5 planets (b) | 9.4 | 3.55$\pm$0.25 | 1.91$\pm$0.26 | 0.61$\pm$0.25 | 17.6$\pm$0.3 | 1.3$\pm$0.3
Joint fit | 3.56$\pm$0.24 | 1.93$\pm$0.24 | 0.61$\pm$0.13 | 17.75$\pm$0.15 | 1.53$\pm$0.23

777 (a) Using an uninformative prior for the periodic parameter. (b) Using a
normal prior for the periodic parameter centred at 300 d.

To crosscheck our results, we also analysed the RV dataset using Exo-
Striker888https://github.com/3fon3fonov/exostriker (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 $\sim$300 days,
which is the signal we named HD 191939 g. Furthermore, the signal at 17.7 d
found by L22 (Sect. 8 therein) is observed in the residuals in Fig. 2 but at a
very low significance level (FAP$\gg$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 $\sim$300 d (see
Fig. 11). 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 $\sim$300d signal. Here, we used 362 RVs with a baseline of 677
days, which cover two complete periods.

To unveil the nature of the $\sim$300 d 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 d.
The priors and posteriors of the hyperparameters used in the GP models are
shown in Table 6. Table 2 presents the Bayesian log-evidence for the explored
models and the measured $K$ for each of the fitted signals. The derived $K$s
for the different models are consistent within errors, ensuring that the
planet mass is not model-dependent or affected by the $\sim$300 d 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
$\Delta\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.
10). 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.

###  3.3 Joint fit

Table 3:  Parameters and $1\sigma$ uncertainties for the juliet joint fit
model for HD 191939 planetary system. Priors and description for each
parameter are presented in Table 7. The adopted stellar properties used to
derive the planetary parameters are the ones from L22 in Table 1. The model to
estimate planet f limits is explained in Sect. 3.4.

Parameter | Planet b | Planet c | Planet d | Planet e | Planet g | Planet f
---|---|---|---|---|---|---
$P$ [d] | 8.8803256 (47) | 28.579743 (45) | 38.353037 (60) | 101.12$\pm$0.13 | 284${}^{+10}_{-8}$ | ¿2200
$t_{0}$ (a) | 2443.54236${}^{+0.00025}_{-0.00023}$ | 2440.5491${}^{+0.0011}_{-0.0007}$ | 2433.906${}^{+0.0005}_{-0.0009}$ | 2348.12$\pm$0.50 | 2385${}^{+13}_{-10}$ | –
$K$ [$\mathrm{m\,s^{-1}}$] | 3.56${}^{+0.21}_{-0.24}$ | 1.93${}^{+0.23}_{-0.24}$ | 0.61$\pm$0.13 | 17.73${}^{+0.14}_{-0.16}$ | 1.53${}^{+0.21}_{-0.23}$ | ¿36
ecc | 0.031${}^{+0.010}_{-0.011}$ | 0.034${}^{+0.034}_{-0.013}$ | 0.031${}^{+0.018}_{-0.012}$ | 0.031${}^{+0.008}_{-0.016}$ | 0.030${}^{+0.025}_{-0.011}$ | –
$\omega$ (deg) | 5${}^{+40}_{-35}$ |  $-$90${}^{+220}_{-40}$ | 15${}^{+90}_{-160}$ |  $-$130${}^{+180}_{-16}$ | 18$\pm$65 | –
$r_{1}$ | 0.739${}^{+0.016}_{-0.020}$ | 0.754${}^{+0.031}_{-0.020}$ | 0.626${}^{+0.035}_{-0.033}$ | – | – | –
$r_{2}$ | 0.03319${}^{+0.00032}_{-0.00017}$ | 0.03118${}^{+0.00027}_{-0.00034}$ | 0.02916${}^{+0.00034}_{-0.00021}$ | – | – | –
Derived parameters
$p={R}_{\mathrm{p}}/{R}_{\star}$ | 0.03319${}^{+0.00032}_{-0.00017}$ | 0.03118${}^{+0.00027}_{-0.00034}$ | 0.02916${}^{+0.00034}_{-0.00021}$ | – | – | –
$b=(a/{R}_{\star})\cos{i_{\mathrm{p}}}$ | 0.610${}^{+0.025}_{-0.030}$ | 0.630${}^{+0.050}_{-0.030}$ | 0.440$\pm$0.050 | – | – | –
$a/{R}_{\star}$ | 18.36${}^{+0.50}_{-0.33}$ | 40.0${}^{+1.1}_{-0.7}$ | 48.7${}^{+1.3}_{-0.9}$ | 92.9${}^{+2.5}_{-1.7}$ | 186.0${}^{+5.0}_{-5.5}$ | ¿730
$i_{\mathrm{p}}$ (deg) | 88.10${}^{+0.14}_{-0.10}$ | 89.10${}^{+0.06}_{-0.08}$ | 89.49${}^{+0.05}_{-0.08}$ | – | – | –
$t_{\mathrm{T}}$ [h] | 2.939${}^{+0.036}_{-0.051}$ | 4.162${}^{+0.24}_{-0.096}$ | 5.36${}^{+0.19}_{-0.11}$ | – | – | –
$R_{\mathrm{p}}$ [${R}_{\oplus}$] | 3.410$\pm$0.075 | 3.195$\pm$0.075 | 2.995$\pm$0.070 | – | – | –
$M_{\mathrm{p}}$ [${M}_{\oplus}$] (b) | 10.00$\pm$0.70 | 8.0$\pm$1.0 | 2.80$\pm$0.60 | ¿112.2$\pm$4.0 | ¿13.5$\pm$2.0 | ¿660
$\rho_{\mathrm{p}}$ [$\mathrm{g\,cm^{-3}}$] | 1.40${}^{+0.15}_{-0.13}$ | 1.35$\pm$0.20 | 0.57$\pm$0.13 | – | – | –
$g_{\mathrm{p}}$ [$\mathrm{m\,s^{-2}}$] | 8.4${}^{+0.8}_{-0.7}$ | 7.7${}^{+1.1}_{-1.0}$ | 3.1$\pm$0.7 | – | – | –
$a_{\mathrm{p}}$ [AU] | 0.0804${}^{+0.0025}_{-0.0023}$ | 0.1752${}^{+0.0055}_{-0.0050}$ | 0.2132${}^{+0.0065}_{-0.0061}$ | 0.407$\pm$0.012 | 0.812$\pm$0.028 | ¿3.2
$T_{\mathrm{eq}}$ [K] (c) | 880$\pm$20 | 600$\pm$13 | 540$\pm$11 | 390$\pm$8 | 278$\pm$6 | ¡125
${S}$ [$\mathrm{S}_{\oplus}$] | 100$\pm$7 | 21.0$\pm$1.4 | 14.3$\pm$1.0 | 3.9$\pm$0.3 | 0.99$\pm$0.08 | –
Model Parameters
| Stellar density |  |
$\rho_{\star}$ [kg m-3] | 1485${}^{+120}_{-80}$ |  |
| Photometry parameters |  |
$\mu_{\textit{TESS}}$ (ppm) | $-$26.2$\pm$1.0 |  |
$\sigma_{\textit{TESS}}$ (ppm) | 152.2${}^{+1.9}_{-1.3}$ |  |
$q_{1,\textit{TESS}}$ | 0.304${}^{+0.051}_{-0.052}$ |  |
$q_{2,\textit{TESS}}$ | 0.4${}^{+0.18}_{-0.12}$ |  |
| RV parameters |  |
Slope $\dot{\gamma}$ [$\mathrm{m\,s^{-1}\,d^{-1}}$] | $0.1902^{+0.0022}_{-0.0020}$ |  |
Curve $\ddot{\gamma}$ [$\mathrm{m\,s^{-1}\,d^{-2}}$] | $-2.27^{+0.07}_{-0.08}\times 10^{-5}$ |  |
$\gamma_{\mathrm{APF}}$ [$\mathrm{m\,s^{-1}}$] | $-6.1^{+1.6}_{-1.4}$ |  |
$\sigma_{\mathrm{APF}}$ [$\mathrm{m\,s^{-1}}$] | $3.48^{+0.35}_{-0.38}$ |  |
$\gamma_{\mathrm{HIRES}}$ [$\mathrm{m\,s^{-1}}$] | $-1.4^{+1.7}_{-1.6}$ |  |
$\sigma_{\mathrm{HIRES}}$ [$\mathrm{m\,s^{-1}}$] | $1.95^{+0.16}_{-0.15}$ |  |
$\gamma_{\mathrm{CARMENES}}$ [$\mathrm{m\,s^{-1}}$] | $4.5\pm$1.4 |  |
$\sigma_{\mathrm{CARMENES}}$ [$\mathrm{m\,s^{-1}}$] | $4.24^{+0.33}_{-0.26}$ |  |
$\gamma_{\mathrm{HARPS-N}}$ [$\mathrm{m\,s^{-1}}$] | $-5.0^{+1.4}_{-1.5}$ |  |
$\sigma_{\mathrm{HARPS-N}}$ [$\mathrm{m\,s^{-1}}$] | $1.76^{+0.15}_{-0.22}$ |  |

999(a) Central time of transit ($t_{0}$) units are BJD $-$ 2 457 000. (b) The
masses for planets e, g, and f are a lower limit
($M_{\mathrm{p}}\sin{i_{\mathrm{p}}}$) because they are only detected in the
RV data. (c) Equilibrium temperatures were calculated assuming zero Bond
albedo.

Figure 3:  Photometry data phase-folded to the period $P$ and central time of
transit $t_{0}$ (shown above each panel, $t_{0}$ 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.

Figure 4:  RV data phase-folded to the period $P$ and central time of transit
$t_{0}$ (shown above each panel, $t_{0}$ 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\sigma$
confidence interval (shaded grey area). The error bars include the
instrumental jitter term added in quadrature.

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}_{\mathrm{p}}/{R}_{\star}$), 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 ($\rho_{\star}$) from
L22 in Table 1 to set a normal prior for $\rho_{\star}$. We fixed the dilution
factor to 1 based on BA20 and we added a relative flux offset ($\mu$) and a
jitter term ($\sigma$) 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 $\alpha$ $=$
$1.52$ and $\beta$ $=$ $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 ($\Delta\ln{Z}$ = $\ln{Z_{\rm
ecc}}-\ln{Z_{\rm no\leavevmode\nobreak\ ecc}}>8$) and the results from both
models are coincident within their uncertainties.

The priors used in the joint fit are listed in Table 7. 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\sigma$. 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 Figures 3 and 4, respectively. Photometric
and RV time series along with the best-fit models are shown in Figures 12 and
13, and Figure 14, 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 ($\Delta
F\simeq({R}_{\mathrm{p}}/{R}_{\star})^{2}$) and transit duration
${t}_{\mathrm{T}}$ of planets e and g with their predicted radii (see Sect.
4.1). HD 191939 e would produce a flux decrease of $\sim$1.6% over $\sim$8 h.
According to the derived ephemeris, a transit of planet e was expected during
Sector 16 (see Fig. 12) 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 ($\sim$1300 ppm). However, the transit
would span $\sim$12 h due to its large period. Planet g was expected to
transit at some point in Sectors 18 and 19 (see Fig. 12). Although there is no
clear evidence of a complete transit, ingress, or egress during these sectors,
the $1\sigma$ uncertainty for the transit midpoint ($\sim$16 d) comprises
observing gaps where the transit could have happened. These observing gaps
represent 18% of the $\pm\,1\sigma$ 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.

###  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 ${M}_{\mathrm{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 d, 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. 15) and the phase-
folded RVs for planet-candidate f(Fig. 16) also confirm that HD 191939 f
properties are not constrained. The period of $\gtrsim$2200 d is in agreement
with the 1700–7200 d period range, and the semi-amplitude of $\gtrsim$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.

##  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 284${}^{+10}_{-8}$
d. We derived a minimum mass of 13.5$\pm$2.0 ${M}_{\oplus}$ 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\sigma$ 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}_{\oplus}$ up to 1.5 ${M}_{\mathrm{J}}$, and 1
${R}_{\oplus}$ up to 1.5 ${R}_{\mathrm{J}}$, along with theoretical
composition models of Zeng et al.
(2019)101010https://lweb.cfa.harvard.edu/~lzeng/planetmodels.html.

Figure 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}_{\oplus}$ up to 1.5 ${M}_{\mathrm{J}}$ and 1
${R}_{\oplus}$ up to 1.5 ${R}_{\mathrm{J}}$, 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 $\pm
1\sigma$ mass regions of HD 191939 g and e. Temperate planets with
$T_{\mathrm{eq}}$ = 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.1 HD 191939 g: a new planet

Figure 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 $T_{\mathrm{eq}}$ = 273–373 K. HD 191939
is marked at the centre with a black star. A diagram including planet f is
shown in Fig. 17.

With an orbital period of $\sim$280 d and a minimum mass of 13.5
${M}_{\oplus}$, 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 ($\sim$5–20 ${M}_{\oplus}$) and well-
determined masses, namely: HD 31527 d (P$\simeq$274 d, 16$\pm$3
${M}_{\oplus}$; Mayor et al., 2011), HD 10180 g (P$\simeq$602 d, 21$\pm$3
${M}_{\oplus}$; Lovis et al., 2011), GJ 3138 d (P$\simeq$257 d, 10$\pm$2
${M}_{\oplus}$; Astudillo-Defru et al., 2017), Barnard b (P$\simeq$232 d,
3.2$\pm$0.4 ${M}_{\oplus}$; Ribas et al., 2018), GJ 273 d (P$\sim$414 d,
11$\pm$4 ${M}_{\oplus}$) and e (P$\sim$542 d, 9$\pm$4 ${M}_{\oplus}$; Tuomi et
al., 2019), and GJ 687 c (P$\simeq$727 d, 16$\pm$4 ${M}_{\oplus}$; 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 d, $<$4.6 ${M}_{\oplus}$), e
(P$\simeq$369 d, 12$\pm$5 ${M}_{\oplus}$), and f ($P$=542 d, 12$\pm$3
${M}_{\oplus}$; 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 d, 6.4$\pm$0.8 ${M}_{\oplus}$; Ofir et al.,
2014), KOI-1783 c ($P$ = 284 d, 15.0${}^{+4.3}_{-3.6}$ ${M}_{\oplus}$;
Vissapragada et al., 2020), and Kepler-90 g ($P$ = 210 d, 15$\pm$1
${M}_{\oplus}$; 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 implementation111111https://github.com/chenjj2/forecaster. The code
predicted planetary radii for planets g and e of 3.7${}^{+1.5}_{-1.1}$
${R}_{\oplus}$ and 12.9${}^{+4.7}_{-3.9}$ ${R}_{\oplus}$, 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}$ $\mathrm{g\,cm^{-3}}$ and
0.30${}^{+0.70}_{-0.16}$ $\mathrm{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 $\sim$3.4 ${R}_{\oplus}$ 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 $\sim$17
${R}_{\oplus}$, 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_{\mathrm{eq}}$ is
compatible with the presence of liquid water ($T$ $\in$ $[273,373]$ K). With a
semi-major axis of $\sim$0.82 AU, planet g is in the outer edge of the HZ,
which we set at $\sim$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 $\sim$13.5
${M}_{\oplus}$ are above the Earth-like composition line, supporting the idea
that HD 191939 g is likely a gaseous planet. Moreover, its $1\sigma$ 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 $\rho_{\mathrm{p}}$ = 0.15$\pm$0.05 $\mathrm{g\,cm^{-3}}$ (Liang et al.,
2021) and $\rho_{\mathrm{p}}$ = 0.09$\pm$0.02 $\mathrm{g\,cm^{-3}}$ (Santerne
et al., 2019), respectively.

###  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 $\rho_{\mathrm{d}}$ = 0.57$\pm$0.13
$\mathrm{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_{\mathrm{eq}}$
$\simeq$ 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$\pm$0.14 ${R}_{\oplus}$, $\rho_{\mathrm{p}}$ = 0.53$\pm$0.15
$\mathrm{g\,cm^{-3}}$, $T_{\mathrm{eq}}$ $\simeq$ 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 $\rho_{\mathrm{p}}$ = 0.09 and 1.43 $\mathrm{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-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 $\sim$3 ${R}_{\oplus}$ and/or a
$T_{\mathrm{eq}}$ $\sim$ 500–600 K in the database of exoplanet atmospheric
observations of
ExoAtmospheres121212http://research.iac.es/proyecto/exoatmospheres/index.php.
Only the warm sub-Neptune GJ 1214 b (2.74$\pm$0.05 ${R}_{\oplus}$,
$T_{\mathrm{eq}}$ $\simeq$ 596 K; Cloutier et al., 2021) fitted our
conditions, although it is a denser planet ($\rho_{\mathrm{p}}$ =
2.20$\pm$0.17 $\mathrm{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).

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

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 $\mu$m), NIRSpec-G395M
(2.88–5.20 $\mu$m), and MIRI-LRS (5–12 $\mu$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
$\mathcal{R}$ $\sim$ 50, and of 27 ppm for the MIRI-LRS bins with sizes of
0.1–0.2 $\mu$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 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
($\mathcal{R}$ $\sim$ 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 1 000$\times$ solar metallicity or a hazy one
with 100$\times$ 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.

Table 4: Transmission spectroscopy amplitudes at low resolution ($\mathcal{R}$
$\sim$ 170) for various models of the HD 191939 d atmosphere and nominal
wavelength ranges of HST and JWST instrumental modes. These include clear or
hazy atmospheres with scaled solar metallicities. We adopted the formalisms of
Lee et al. (2013) for the hazy models, where $\alpha$ denotes the particle
size in $\mu$m, the mixing ratio is $\chi_{c}$ = 10-12, and the extinction
coefficient is $Q_{0}$ = 40. Lines in bold correspond to the synthetic spectra
shown in Figure 7.

Haze | Met. | Spectroscopic amplitudes (ppm)
---|---|---
| ($[M/H]_{\oplus}$) | 1.1–1.7 (a) | 0.6–2.8 (b) | 2.88–5.20 (c) | 5–12 (d) $\mu$m
Clear | 1 | 460 | 664 | 860 | 774
Clear | 10 | 490 | 665 | 670 | 643
Clear | 100 | 195 | 348 | 294 | 219
Clear | 1000 | 38 | 77 | 59 | 37
$\alpha=$ 0.05 | 1 | 201 | 400 | 658 | 597
$\alpha=$ 0.05 | 100 | 144 | 213 | 276 | 206
$\alpha=$ 0.10 | 1 | 174 | 693 | 380 | 400
$\alpha=$ 0.10 | 100 | 57 | 181 | 250 | 190
131313Nominal 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.

###  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 d, Pc=28.58 d, Pd=38.35 d). 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 ($P_{\mathrm{e}}$ $=$ 101
d, $P_{\mathrm{g}}$ $\simeq$ 280 d). 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_{\mathrm{J}}$ in a very
long-period orbit ($\gg$3000 d). 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}_{\oplus}$ and a
period of 259 d. 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
Archive141414https://exoplanetarchive.ipac.caltech.edu/ database looking for
systems with 2–4 intermediate planets (2–25 ${M}_{\oplus}$ or 2–8
${R}_{\oplus}$) interior to a gas-giant planet (¿50 ${M}_{\oplus}$ or ¿8
${R}_{\oplus}$) 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).

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

KOI-94 has two inner planets with masses of 10.5$\pm$4.6 ${M}_{\oplus}$ and
¡21.3 ${M}_{\oplus}$, and then a warm Saturn-like planet (106$\pm$11
${M}_{\oplus}$, 11.2$\pm$1.1 ${R}_{\oplus}$) and an outer planet slightly more
massive than the inner ones (35${}^{+18}_{-28}$ ${M}_{\oplus}$, 6.6$\pm$0.6
${R}_{\oplus}$), but without evidence for a long-period companion like HD
191939 f. KOI-94 planets are gaseous with low densities ranging between
$\sim$0.35 and 1 $\mathrm{g\,cm^{-3}}$ (except for KOI-94 b; 10.1$\pm$5.5
$\mathrm{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 (Pb $=$ 3.7 d, Pc $=$ 10.4 d, Pd $=$ 22.3 d), with
the giant planet also close to 5:2 with the outermost planet KOI-94 e (Pe $=$
54 d). 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 Figure 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 $\sim$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 $\sim$0.82 au (HD 191939 g semi-major
axis), which is comparable to the Venus orbital distance of 0.72 au.

##  5 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 $\sim$2 years, allowed the detection of a new non-
transiting planetary signal (HD 191939 g) with a period of $\sim$280 d and a
minimum mass of $M_{\mathrm{p}}\sin i$ $=$ 13.5$\pm$2.0 ${M}_{\oplus}$. 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\sigma$ level of significance, for which only an upper limit was known.

We determine a mass for HD 191939 d of 2.80$\pm$0.60 ${M}_{\oplus}$, leading
to a mean bulk density of $\rho_{\mathrm{d}}$ $=$ 0.57 $\mathrm{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 d, 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 research has made use of data obtained from or
tools provided by the portal exoplanet.eu of The Extrasolar Planets
Encyclopaedia. JK gratefully acknowledge the support of the Swedish National
Space Agency (SNSA; DNR 2020-00104). PK is 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 acknowledge 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
acknowledge the inspiring discussions with Yess, Alejandro, his colleagues,
and friends.

## References

  * Abel et al. (2011) Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2011, Journal of Physical Chemistry A, 115, 6805
  * Abel et al. (2012) Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2012, J. Chem. Phys., 136, 044319
  * Agúndez et al. (2012) Agúndez, M., Venot, O., Iro, N., et al. 2012, A&A, 548, A73
  * Al-Refaie et al. (2021) Al-Refaie, A. F., Changeat, Q., Waldmann, I. P., & Tinetti, G. 2021, ApJ, 917, 37
  * Alam et al. (2022) Alam, M. K., Kirk, J., Dressing, C. D., et al. 2022, ApJ, 927, L5
  * Ambikasaran et al. (2014) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2014
  * Astudillo-Defru et al. (2017) Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017, A&A, 602, A88
  * Badenas-Agusti et al. (2020) Badenas-Agusti, M., Günther, M. N., Daylan, T., et al. 2020, AJ, 160, 113
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373
  * Bean et al. (2010) Bean, J. L., Miller-Ricci Kempton, E., & Homeier, D. 2010, Nature, 468, 669
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Caballero et al. (2016) Caballero, J. A., Guàrdia, J., López del Fresno, M., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9910, Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, 99100E
  * Cabrera et al. (2014) Cabrera, J., Csizmadia, S., Lehmann, H., et al. 2014, ApJ, 781, 18
  * Cannon & Pickering (1993) Cannon, A. J. & Pickering, E. C. 1993, VizieR Online Data Catalog, III/135A
  * Chachan et al. (2020) Chachan, Y., Jontof-Hutter, D., Knutson, H. A., et al. 2020, AJ, 160, 201
  * Chen & Kipping (2017) Chen, J. & Kipping, D. 2017, ApJ, 834, 17
  * Cloutier et al. (2021) Cloutier, R., Charbonneau, D., Deming, D., Bonfils, X., & Astudillo-Defru, N. 2021, AJ, 162, 174
  * Cosentino et al. (2012) Cosentino, R., Lovis, C., Pepe, F., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84461V
  * Cosentino et al. (2014) Cosentino, R., Lovis, C., Pepe, F., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, 91478C
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Dragomir et al. (2019) Dragomir, D., Teske, J., Günther, M. N., et al. 2019, ApJ, 875, L7
  * Espinoza (2018) Espinoza, N. 2018, Research Notes of the American Astronomical Society, 2, 209
  * Espinoza et al. (2019) Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262
  * Espinoza et al. (2022) Espinoza, N., Pallé, E., Kemmer, J., et al. 2022, AJ, 163, 133
  * Feng et al. (2020) Feng, F., Shectman, S. A., Clement, M. S., et al. 2020, ApJS, 250, 29
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Fletcher et al. (2018) Fletcher, L. N., Gustafsson, M., & Orton, G. S. 2018, ApJS, 235, 24
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Angus, R., & Ambikasaran, S. 2017, AJ, 154, 220
  * Fortney et al. (2013) Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gao & Zhang (2020) Gao, P. & Zhang, X. 2020, ApJ, 890, 93
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39
  * Günther et al. (2019) Günther, M. N., Pozuelos, F. J., Dittmann, J. A., et al. 2019, Nature Astronomy, 3, 1099
  * Hadden & Lithwick (2017) Hadden, S. & Lithwick, Y. 2017, AJ, 154, 5
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, ed. G. Chiozzi & J. C. Guzman, 99133E
  * Jontof-Hutter et al. (2014) Jontof-Hutter, D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, ApJ, 785, 15
  * Kasper et al. (2020) Kasper, D., Bean, J. L., Oklopčić, A., et al. 2020, AJ, 160, 258
  * Kempton et al. (2018) Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152
  * Kreidberg (2015) Kreidberg, L. 2015, Publications of the Astronomical Society of the Pacific, 127, 1161
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69
  * Lafarga et al. (2020) Lafarga, M., Ribas, I., Lovis, C., et al. 2020, A&A, 636, A36
  * Lee et al. (2013) Lee, J.-M., Heng, K., & Irwin, P. G. J. 2013, ApJ, 778, 97
  * Liang et al. (2021) Liang, Y., Robnik, J., & Seljak, U. 2021, AJ, 161, 202
  * Libby-Roberts et al. (2020) Libby-Roberts, J. E., Berta-Thompson, Z. K., Désert, J.-M., et al. 2020, AJ, 159, 57
  * Lovis et al. (2011) Lovis, C., Ségransan, D., Mayor, M., et al. 2011, A&A, 528, A112
  * Lubin et al. (2022) Lubin, J., Van Zandt, J., Holcomb, R., et al. 2022, AJ, 163, 101
  * Luque et al. (2022) Luque, R., Fulton, B. J., Kunimoto, M., et al. 2022, A&A, 664, A199
  * Luque et al. (2021) Luque, R., Serrano, L. M., Molaverdikhani, K., et al. 2021, A&A, 645, A41
  * Masuda (2014) Masuda, K. 2014, ApJ, 783, 53
  * Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv e-prints, arXiv:1109.2497
  * Mills et al. (2019) Mills, S. M., Howard, A. W., Weiss, L. M., et al. 2019, AJ, 157, 145
  * Morello et al. (2021) Morello, G., Zingales, T., Martin-Lagarde, M., Gastaud, R., & Lagage, P.-O. 2021, AJ, 161, 174
  * Morris et al. (2017) Morris, R. L., Twicken, J. D., Smith, J. C., et al. 2017, Kepler Data Processing Handbook: Photometric Analysis, Kepler Science Document KSCI-19081-002
  * Ofir et al. (2014) Ofir, A., Dreizler, S., Zechmeister, M., & Husser, T.-O. 2014, A&A, 561, A103
  * Ohno & Tanaka (2021) Ohno, K. & Tanaka, Y. A. 2021, ApJ, 920, 124
  * Orell-Miquel et al. (2022) Orell-Miquel, J., Murgas, F., Pallé, E., et al. 2022, A&A, 659, A55
  * Petit dit de la Roche et al. (2020) Petit dit de la Roche, D. J. M., van den Ancker, M. E., & Miles-Paez, P. A. 2020, Research Notes of the American Astronomical Society, 4, 231
  * Quinn et al. (2019) Quinn, S. N., Becker, J. C., Rodriguez, J. E., et al. 2019, AJ, 158, 177
  * Quirrenbach et al. (2014) Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91471F
  * Quirrenbach et al. (2020) Quirrenbach, A., CARMENES Consortium, Amado, P. J., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114473C
  * Ribas et al. (2018) Ribas, I., Tuomi, M., Reiners, A., et al. 2018, Nature, 563, 365
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Santerne et al. (2019) Santerne, A., Malavolta, L., Kosiarek, M. R., et al. 2019, arXiv e-prints, arXiv:1911.07355
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000
  * Soubiran et al. (2018) Soubiran, C., Jasniewicz, G., Chemin, L., et al. 2018, A&A, 616, A7
  * Southworth (2011) Southworth, J. 2011, MNRAS, 417, 2166
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Stassun et al. (2018) Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102
  * Stassun & Torres (2018) Stassun, K. G. & Torres, G. 2018, ApJ, 862, 61
  * Stumpe et al. (2014) Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985
  * Thorngren et al. (2016) Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64
  * Trifonov (2019) Trifonov, T. 2019, The Exo-Striker: Transit and radial velocity interactive fitting tool for orbital analysis and N-body simulations
  * Trifonov et al. (2020) Trifonov, T., Tal-Or, L., Zechmeister, M., et al. 2020, A&A, 636, A74
  * Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71
  * Tsiaras et al. (2018) Tsiaras, A., Waldmann, I. P., Zingales, T., et al. 2018, AJ, 155, 156
  * Tuomi et al. (2019) Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, arXiv e-prints, arXiv:1906.04644
  * Van Eylen & Albrecht (2015) Van Eylen, V. & Albrecht, S. 2015, ApJ, 808, 126
  * Van Eylen et al. (2019) Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, AJ, 157, 61
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
  * Vanderburg et al. (2016) Vanderburg, A., Becker, J. C., Kristiansen, M. H., et al. 2016, ApJ, 827, L10
  * Vissapragada et al. (2020) Vissapragada, S., Jontof-Hutter, D., Shporer, A., et al. 2020, AJ, 159, 108
  * Vogt et al. (2014) Vogt, S. S., Radovan, M., Kibrick, R., et al. 2014, PASP, 126, 359
  * Weiss et al. (2013) Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14
  * Wolfgang et al. (2016) Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19
  * Xie et al. (2016) Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Science, 113, 11431
  * Yu et al. (2021) Yu, X., He, C., Zhang, X., et al. 2021, Nature Astronomy, 5, 822
  * Zechmeister & Kürster (2009) Zechmeister, M. & Kürster, M. 2009, A&A, 496, 577
  * Zechmeister et al. (2018) Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Science, 116, 9723

##  Appendix A Additional figures and tables

Table 5:  Prior and posterior distributions for PDC-SAP detrending fit. Prior labels $\mathcal{F}$, $\mathcal{N}$, and $\mathcal{J}$ represent fixed, normal, and Jeffrey’s distributions, respectively.  Parameter | Prior | Posterior
---|---|---
$D_{\textit{TESS}}$ | $\mathcal{F}(1)$ | –
$\mu_{\textit{TESS}}$ (ppm) | $\mathcal{N}(0,0.1)$ | 20${}^{+110}_{-100}$
$\sigma_{\textit{TESS}}$ (ppm) | $\mathcal{J}(10^{-6},10^{6})$ | 135${}^{+1}_{-2}$
$\sigma_{\textit{GP}}$ (ppm) | $\mathcal{J}(10^{-6},10^{6})$ | 970${}^{+55}_{-35}$
$\rho_{\textit{GP}}$ [d] | $\mathcal{J}(10^{-3},10^{3})$ | 10.2$\pm$0.5

Figure 9:  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.

Figure 10:  Posterior distribution of hyperparameters from the uninformative (left) and constrained (right) GPqp models used in Sect. 3.2. Prior distributions are in Table 6.  Figure 11:  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.  Table 6:  Prior and posterior distributions for the RV GP models explored in Sect. 3.2. Prior labels $\mathcal{U}$, $\mathcal{N}$, and $\mathcal{J}$ represent uniform, normal, and Jeffrey’s distributions, respectively. Corner plots for the quasi-periodic model hyperparameters are shown in Fig. 10.  Hyperparameter | Prior | Posterior
---|---|---
Exponential model
$\sigma$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(0,100)$ | 3.3${}^{+0.8}_{-0.7}$
$\tau$ [d] | $\mathcal{J}(0,10^{3})$ | 1.3${}^{+0.8}_{-0.3}$
Matern model
$\sigma$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(0,100)$ | 1.7${}^{+0.5}_{-0.4}$
$\rho$ [d] | $\mathcal{J}(0,10^{3})$ | 15${}^{+27}_{-11}$
Quasi-periodic models
B [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(0,100)$ | 4.02${}^{+2.7}_{-1.4}$
C [$\mathrm{m\,s^{-1}}$] | $\mathcal{J}(10^{-3},10)$ | 0.11 ${}^{+1.85}_{-0.11}$
L [d] | $\mathcal{J}(0.1,10^{3})$ | 13 ${}^{+131}_{-10}$
Prot [d] | $\mathcal{N}(300,50)$ | 297$\pm$44
B [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(0,100)$ | 4.2${}^{+2.8}_{-1.5}$
C [$\mathrm{m\,s^{-1}}$] | $\mathcal{J}(10^{-3},10)$ | 0.15 ${}^{+2.17}_{-0.14}$
L [d] | $\mathcal{J}(0.1,10^{3})$ | 11 ${}^{+70}_{-7}$
Prot [d] | $\mathcal{U}(150,600)$ | 390${}^{+130}_{-145}$
Table 7:  Priors for each parameter used in the juliet joint fit model for HD
191939 planetary system. Prior labels $\mathcal{F}$, $\mathcal{U}$,
$\mathcal{N}$, $\mathcal{B}$, and $\mathcal{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) are both explained in Section 3.3.

Parameter | Planet b | Planet c | Planet d | Planet e | Planet g
---|---|---|---|---|---
$P$ [d] | $\mathcal{N}(8.8803,0.0005)$ | $\mathcal{N}(28.5795,0.0005)$ | $\mathcal{N}(38.3531,0.0005)$ | $\mathcal{N}(101.0,2.0)$ | $\mathcal{N}(300.0,20.0)$
$t_{0}$ (a) | $\mathcal{N}(2443.5414,0.005)$ | $\mathcal{N}(2440.5455,0.005)$ | $\mathcal{N}(2433.9098,0.005)$ | $\mathcal{N}(2347.5,2.0)$ | $\mathcal{N}(2400.0,50.0)$
$K$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(0,25)$ | $\mathcal{U}(0,25)$ | $\mathcal{U}(0,25)$ | $\mathcal{U}(0,25)$ | $\mathcal{U}(0,25)$
ecc | $\mathcal{B}(1.52,29)$ | $\mathcal{B}(1.52,29)$ | $\mathcal{B}(1.52,29)$ | $\mathcal{B}(1.52,29)$ | $\mathcal{B}(1.52,29)$
$\omega$ (deg) | $\mathcal{U}(-180,180)$ | $\mathcal{U}(-180,180)$ | $\mathcal{U}(-180,180)$ | $\mathcal{U}(-180,180)$ | $\mathcal{U}(-180,180)$
$r_{1}$ | $\mathcal{U}(0,1)$ | $\mathcal{U}(0,1)$ | $\mathcal{U}(0,1)$ | – | –
$r_{2}$ | $\mathcal{U}(0,1)$ | $\mathcal{U}(0,1)$ | $\mathcal{U}(0,1)$ | – | –
Model Parameters
Stellar density
$\rho_{\star}$ [kg m-3] | $\mathcal{N}(1370,150)$ |
Photometry parameters
$D_{\textit{TESS}}$ | $\mathcal{F}(1)$ |
$\mu_{\textit{TESS}}$ (ppm) | $\mathcal{N}(0,0.1)$ |
$\sigma_{\textit{TESS}}$ (ppm) | $\mathcal{J}(0.1,1000)$ |
$q_{1,\textit{TESS}}$ | $\mathcal{U}(0,1)$ |
$q_{2,\textit{TESS}}$ | $\mathcal{U}(0,1)$ |
RV parameters
Intercept $\gamma$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{F}(-300)$ |
Slope $\dot{\gamma}$ [$\mathrm{m\,s^{-1}\,d^{-1}}$] | $\mathcal{U}(-1,1)$ |
Curve $\ddot{\gamma}$ [$\mathrm{m\,s^{-1}\,d^{-2}}$] | $\mathcal{U}(-0.1,0.1)$ |
$\gamma_{\mathrm{APF}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(-10,10)$ |
$\sigma_{\mathrm{APF}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{J}(0.1,10)$ |
$\gamma_{\mathrm{HIRES}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(-10,10)$ |
$\sigma_{\mathrm{HIRES}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{J}(0.1,10)$ |
$\gamma_{\mathrm{CARMENES}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(-10,10)$ |
$\sigma_{\mathrm{CARMENES}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{J}(0.1,10)$ |
$\gamma_{\mathrm{HARPS-N}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{U}(-10,10)$ |
$\sigma_{\mathrm{HARPS-N}}$ [$\mathrm{m\,s^{-1}}$] | $\mathcal{J}(0.1,10)$ |

151515 (a) Central time of transit ($t_{0}$) units are BJD $-$ 2 457 000.

Figure 12:  TESS photometry from Sectors 15–19, 21, 22, 24, and 25 along with
the best-fit model (see Fig. 13 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 $t_{0}$ and
$\pm\,1\sigma$ uncertainty for the non-transiting planets HD 191939 e
(magenta) and g (orange). Figure 13:  As in Fig. 12 but for Sectors 41 and 48.
Figure 14:  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\sigma$ confidence interval (shaded
grey area). The error bars include the instrumental jitter term added in
quadrature.  Figure 15:  As in Fig. 14 but considering a Keplerian signal for
planet f (see Sect. 3.4).  Figure 16:  RVs phase-folded to the $P$ and $t_{0}$
(shown above the panel, $P$ units are days and $t_{0}$ units are
BJD$-2\,457\,000$) for planet f along with the best-fit model (black line) and
the $1\sigma$, $2\sigma$, and $3\sigma$ confidence intervals (shaded grey
areas).  Figure 17:  As in Fig. 6 but including planet f. Planet f parameters
are from Table 3.