© V. Bourrier et al. 2022

1 Introduction

Precise knowledge of a planet mass and radius is essential to infer its internal structure and the presence of an atmosphere. This is especially relevant for small exoplanets (below ~3 R_⊕), which could encompass a wide range of compositions from mini-Neptunes with volatile H/He envelopes, to ocean planets with water mantle and steam atmospheres, to ultrahot rocky planets with molten lava-rich surfaces and heavyweight envelopes (e.g., Winn et al. 2018; Otegi et al. 2020).

In that respect the HD 3167 system is of particular interest, as it hosts three known planets (Vanderburg et al. 2016; Christiansen et al. 2017; Gandolfi et al. 2017): HD 3167b (P = 0.96 d, R_⊕, M_p = 5.02 ± 0.38 M_⊕), HD 3167d (P = 8.51 d, M_p sin i = 6.90 ± 0.71M_⊕), and HD 3167c (P = 29.84 days, R_⊕, M_⊕). Planets b and c are transiting their nearby (47 pc) and bright (V = 9) K0V star, allowing for detailed measurements of their radius and atmospheric properties. The intermediate planet d is not transiting, but is nonetheless expected to have a low mutual inclination with planet c based on dynamical calculations (Dalal et al. 2019). The orbital architecture of the HD 3167 system is particularly intriguing, because the orbital plane of its innermost planet b is close to the stellar equatorial plane and perpendicular to the orbital planes of the outer planets d and c, on polar orbits around the star (Dalal et al. 2019; Bourrier et al. 2021).

HD 3167 c is a mini-Neptune, that is an exoplanet smaller than Neptune still harboring a substantial volatile envelope of hydrogen and helium, or possibly a large fraction of water (e.g., Mousis et al. 2020). Transit observations in the near-IR with the Wide Field Camera 3 (WFC3) onboard the Hubble Space Telescope (HST), combined with broadband transit measurements with Kepler/K2 and Spitzer/IRAC, revealed molecular absorption in the atmosphere of HD 3167 c (Guilluy et al. 2021; Mikal-Evans et al. 2021).

HD 3167 b belongs to the population of ultra-short period planets (USPs, R < 2R_⊕; P < 1 d). These exotic worlds receive so much energy from their parent star, via irradiation, tides, or even electromagnetic induction (Kislyakova et al. 2017, 2018) that they cannot retain a volatile atmosphere (e.g., Lopez 2017; Winn et al. 2018). The planetary lithosphere is further expected to weaken and melt, leading to the formation of magma oceans and potential volcanic activity, particularly on the irradiated day side (e.g. Schaefer & Fegley 2009; Barnes et al. 2010; Gelman et al. 2011; Léger et al. 2011; Elkins-Tanton 2012). When planetary equilibrium temperatures exceed ~1000 K, outgassing can release massive amounts of dust and metals effluents and sustain a secondary atmosphere around the planet (Rappaport et al. 2012, 2014; Ito et al. 2015). Guenther & Kislyakova (2020) used high-resolution spectroscopy with UVES to search for the absorption lines of a metal-rich envelope during the transit of HD 3167 b. While they were only able to set upper limits on tracers of this envelope (such as sodium and oxygen, Miguel et al. 2011), this does not preclude the possibility that their signatures vary over time (like 55 Cnc e, Ridden-Harper et al. 2016) or that the planet is surrounded by an envelope whose refractory content has a broadband optical signature. Measuring with high precision the radius of HD 3167b, and of USPs in general, in combination with their mass is thus essential to investigate the presence of such exotic atmospheres and gain further insights into their mysterious nature.

The CHEOPS satellite (Benz et al. 2021) was used to characterize, with a high precision, the transits of the ultra-short period planet HD 3167b. We combined these observations with published transit photometry of HD 3167b and c, and with both published and new radial velocity (RV) data of the system (Sect. 2), to perform a global analysis of all data available on the system (Sect. 3), and obtain a complete and refined view of its planets bulk properties (Sect. 4). These properties are used to constrain the internal structures of the two transiting planets in the system in Sect. 5, and to simulate their past dynamical (Sect. 6) and atmospheric (Sect. 7) evolution.

2 Observations and data reduction

2.1 CHEOPS photometry

Transits of HD 3167b were observed with CHEOPS within the frame of the Guaranteed Time Observation (GTO), as part of a sub-program dedicated to measuring with high precision the radius of USPs and better constrain their internal structure. Twelve visits were obtained between 2 August and 14 November 2020. We scheduled three CHEOPS orbits per visit, so that one orbit would cover the transit and two orbits would cover the pre-transit and post-transit phases, allowing us to measure the baseline stellar flux and detrend the light curve. We set an exposure time of 36.5 s with which we expected to reach a precision of about 7% on the transit depth in one visit. We eventually obtained an average precision of 15% per visit due to a higher noise level than anticipated, which we partly attributed to the presence of new hot pixels in the photometric aperture. We note that no transit of planet c was caught by the CHEOPS observations.

Observations were processed with the CHEOPS DRP (Data Reduction Pipeline, version 13.1.0; Hoyer et al. 2020), which performs aperture photometry and provides four light curves extracted with four different aperture radii. The so-called default aperture, with a radius of 25 pixels, had a consistent lowest rms noise level throughout all visits and we selected this data set for our analysis.

2.2 K2 photometry

During the 80-day long campaign 8 of its K2 mission (3 January to 23 March 2016), the Kepler space telescope acquired 30-min cadence photometry of the star HD 3167, from which Vanderburg et al. (2016) discovered the transiting planets b and c.

In addition to the simple aperture photometry (SAP) and the Pre-search Data Conditioning SAP (PDCSAP), the MAST archive¹ provides several high-level science data products (HLSP) based on different photometric extraction techniques. We compared them all and identified the ones leading to light curves with the lowest noise levels: K2SC (Aigrain et al. 2016) and K2SFF (Vanderburg & Johnson 2014). We selected the light curve provided by K2SFF as it preserves the photometric variability in the data. We fitted the data with a joint model of this variability and of the transit light curves, to properly propagate uncertainties throughout all parameters. K2SFF proposes a best solution between photometry for different aperture shapes and sizes. We compared them one by one, confirmed that the proposed best aperture has the highest signal-to-noise ratio, and subsequently used it for our analysis. K2SFF does not provide uncertainties on the photometric points. We thus investigated the error values generated by the other extraction methods and selected the most conservative (largest) ones for our data set, which were computed by K2SC_SAP.

We used the Gaia DR2 catalogue (Gaia Collaboration 2016, 2018) to check for contamination by nearby stars in the aperture. We only considered objects with G-band magnitude differences smaller than nine with HD 3167, and found none within the aperture.

2.3 Spitzer/IRAC photometry

We used Spitzer data from the General Observing (GO) program 13052 (PI: M. Werner) that contain four observations of the HD 3167 system at 4.5 µm (channel 2). These observations covered three transits of planet b (AORs 61072896, 61072640, 68163072) on 22 and 25 October 2016 and 16 October 2018, and one transit of planet c (AOR 61070592) between 31 October and 1 November 2016. The transit of planet c was analyzed in Mikal-Evans et al. (2021).

The data were downloaded from the Spitzer Heritage Archive². We extracted and pre-processed the photometry following a method described in Demory et al. (2016), which relies on the modeling of intra-pixel sensitivity of the IRAC instrument (Ingalls et al. 2016) using the bilinearly-interpolated subpixel sensitivity (BLISS) mapping technique (Stevenson et al. 2012). Our method also includes a correction as a linear function of the full-width at half-maximum (FWHM) of the pixel response function (PRF). An additional log² ramp as a function of time was added for the pre-processing of the transit of planet c. The uncertainties associated with these corrections were propagated to the error bars of the resulting data points. The four resulting de-trended time series were sampled at a cadence of 27 s and we measured a negligible red noise.

2.4 HST photometry

We used five transit observations of HD 3167c collected with HST/WFC3 using the G141 grism configuration (wavelength range 1.1−1.7 µm). These observations are part of the GO program 15333 (PI: I. Crossfield) and they were acquired on 22 May 2018, 20 July 2018, 14 June 2019, 12 August 2019, and 5 July 2020. Each of the five visits covers seven HST orbits.

We used the broadband photometric light curves presented in Mikal-Evans et al. (2021) that we obatined as a three-column file with the time of observation in JD_UTC, the normalized flux, and the normalized uncertainty. The light curves were extracted from the sum of all spectra across the full wavelength range. The resulting data set is made of 69.6-s exposures at a cadence of 111 s. We also downloaded the raw data from the MAST archive³ to have access to the housekeeping parameters.

We converted the time from JD_UTC to BJD_TDB (barycentric Julian date in barycentric dynamical time) using the Python package astropy (Astropy Collaboration 2013, 2018) and assuming that the HST spacecraft is located at the center of the Earth. This approximation leads to a timing error of ±23.08 ms, which we considered negligible. The correction from JD_UTC to BJD_TDB is significant as it leads to correction of up to nearly 10 min.

The light curves feature strong systematics that are typical of HST/WFC3 observations, with one ramp repeatable as a function of HST orbital phase, and one global ramp as a function of time. The flux level also jumps every two points due to the switching from forward to backward scanning of the HST detector between two consecutive exposures. This results in having two mean flux levels that have to be fitted independently with two offsets.

In addition, we carefully checked for the possibility of HD 3167b transiting during these observations and found out that such a double event occurs in visits 3 and 4, and was not reported in previous analysis of these data sets. We therefore included planet b in the analysis of our HST time series.

2.5 Radial velocity data

The RV data analyzed in this work are coming from several instruments. We used data from HARPS programs 097.C-0948 and 098.C-0860, as published in Gandolfi et al. (2017), data from APF/Levy and Keck/HIRES, as published in Christiansen et al. (2017), and data from the HARPS-N GTO programme (Cosentino et al. 2012), as well as from programs A33TAC_15 (PI: D. Gandolfi), CAT16B_61 (PI: H. J. Deeg), A34DDT2 and A36DDT2 (PI: G. Hébrard). The HARPS-N data were partly published in Christiansen et al. (2017), Gandolfi et al. (2017) and Dalal et al. (2019), but we publish here 42 new RV points from the HARPS-N GTO programme.

The HARPS-N RV data were extracted from the instrument raw frames using the latest version of the ESPRESSO data reduction software (DRS, version 2.3.5). Following the work described in Dumusque et al. (2021), the ESPRESSO pipeline has been optimized to work with HARPS-N data. Compared to the HARPS-N DRS version 3.8, the new version of the reduction pipeline, along with the performed optimizations, provide smaller night-to-night variations, estimated to be 0.5 m s⁻¹ rms compared to 0.8 m s⁻¹, and a better stability of the RVs on the long-term, due to a careful selection of Thorium lines used for calibrating the instrument. The new pipeline further extracts RVs from cross-correlation functions computed with improved binary masks (here the G9 mask closer in type to HD 3167), built with weights more representative of the RV information content of each spectral line (Bourrier et al. 2021). We rejected the observation ‘2018-01-08T20-30-55.308’ because it was not possible to correct for the color effect induced by Earth atmospheric diffusion. We removed observation ‘2017-11-14T21-32-32.784’ as well, because the corresponding RV was clearly an outlier of the RV time series, and all stars observed during the same night showed outliers as well. Finally, we removed all observations that were taken at an airmass larger than 1.7, to prevent color dependencies in the RV time series (the ADC corrects for atmospheric extinction up to an airmass of 2.0). This selection gives us a total of 213 HARPN-N RV measurements to analyze.

Merging the data from all four instruments, we obtain a time series of 434 RV data points. To prevent biases induced by the Rossiter-McLaughlin signals from the planets b and c (Dalal et al. 2019; Bourrier et al. 2021), we discarded the 111 data points observed during their transits. The filtered RV data set that we included in our analysis represents a total of 323 data points (39 HARPS points, 102 APF/Levy points, 55 Keck/HIRES points, and 127 HARPS-N points) covering more than 5.3 yr (~1940 days).

In addition to the RV signals, the HARPS and HARPS-N data sets included several stellar activity indicators. The HARPS data included the bisector inverse slope span ( BIS SPAN’) of the cross-correlation function (CCF), the full-width at half maximum (FWHM) of the CCF, and the log . The activity indicators provided with the HARPS-N data were the BIS SPAN’ of the CCF, the FWHM of the CCF, the contrast, the S_MW-index (Lovis et al. 2011), the Hα-index (Gomes da Silva et al. 2011), the Na I lines (Díaz et al. 2007), and the Ca II lines.

3 Global analysis of the system

3.1 Stellar properties

We derived the stellar atmospheric parameters (T_eff, log g, microturbulence, [Fe/H]) using ARES+MOOG, following the same methodology described in Santos et al. (2013); Sousa (2014); Sousa et al. (2021). We used the latest version of ARES⁴ (Sousa et al. 2007, 2015) to measure the equivalent widths (EW) of iron lines on the combined HARPS-N spectrum of HD 3167. We used a minimization process to find ionization and excitation equilibrium and converge to the best set of spectro-scopic parameters. This process makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973). The same method was also applied to a combined spectrum from HARPS observations, providing a completely compatible set of parameters. The stellar abundances [Mg/H] = 0.07 ± 0.03 dex and [Si/H] = 0.00 ± 0.04 dex were derived using the classical curve-of-growth analysis method assuming local thermodynamic equilibrium (e.g., Adibekyan et al. 2012, 2015). The same codes and models were used for the abundance determinations.

We determine the radius of HD 3167 using the infrared flux method (IRFM; Blackwell & Shallis 1977) in a Markov-chain Monte Carlo (MCMC) approach (Schanche et al. 2020). We constructed spectral energy distributions (SEDs) from stellar atmospheric models using the stellar parameters that were derived via the spectral analysis detailed above as priors. These fluxes are compared to observed broadband photometry to derive the apparent bolometric flux, and hence the stellar angular diameter and effective temperature of HD 3167. To achieve this we retrieve data taken from the most recent data releases for the following bandpasses; Gaia G, G_BP, and G_RP, 2MASS J, H, and K, and WISE W1 and W2 (Skrutskie et al. 2006; Wright et al. 2010; Gaia Collaboration 2021) and use stellar atmospheric models from the ATLAS Catalogues (Castelli & Kurucz 2003). The stellar angular diameter is converted to the stellar radius using the offset corrected Gaia EDR3 parallax (Lindegren et al. 2021) from which we obtain R_* = 0.871 ± 0.006 R_⊙.

Stellar mass M_* and age t_* were derived from isochrones starting from T_eff, [Fe/H], and R_*. To make our final estimates more robust we adopted two different stellar evolutionary models, namely PARSEC⁵ v1.2S (Marigo et al. 2017) and CLES (Code Liégeois d’Évolution Stellaire, Scuflaire et al. 2008). In detail, we inferred a first pair of mass and age values by interpolating the input values within pre-computed grids of PARSEC isochrones and tracks through the isochrone placement technique presented in Bonfanti et al. (2015, 2016). To further improve the convergence we also inputted v sin i = 2.41 ± 0.37 km s⁻¹ (Bourrier et al. 2021) into the code to benefit of the synergy between isochrones and gyrochronology as described in Bonfanti et al. (2016). The second pair of mass and age, instead, was inferred by injecting the stellar input values into the CLES code, which retrieves the best-fit output values following the Levenberg-Marquardt minimization scheme (see Salmon et al. 2021, for the details). As thoroughly described in Bonfanti et al. (2021a), we finally merged the two respective pairs of outcomes after carefully checking their mutual consistency through a χ²-based criterion and we obtained M_⊙ and Gyr. Relevant stellar parameters are summarized in Table 1.

3.2 Joint photometry – velocimetry analysis

We performed a joint fit combining all the photometry and velocimetry datasets described in Sect. 2. In the following subsections, we detail how we modeled the planetary signals (transits and RV) consistently in every time series (Sect. 3.2.1), and how we corrected the systematics related to each instrument in the CHEOPS, K2, Spitzer, HST and RV data sets (Sect. 3.2.2 to 3.2.6). Our approach consisted in first performing an analysis of each data set separately to better identify their specificity, and then jointly fitting all the data together (Sect. 3.2.7).

3.2.1 Planetary signals

The transits of planets b and c were modeled using the python package batman (Kreidberg 2015). We selected the quadratic law to describe the effect of the stellar limb darkening, and defined a set of two free coefficients for each of the four instrument passbands. We reduced the number of free parameters by one by using the third Kepler’s law and by fitting for the stellar density ρ_* instead of the normalized semi-major axes a/R_*. With this approach, the planetary properties are fitted while accounting for the stellar unicity. We took advantage of exploiting photometry from four instruments to perform broadband transmission spectroscopy, letting the planet-to-star radii ratios of both transiting planets vary between each passband.

We modeled the RV planetary signals with Keplerian functions, and performed the joint fit of all planets while fitting for the following parameters: the time of inferior conjunction T0, the orbital period P, the combinations of the eccentricity and the argument of periastron e cos ω and e sin ω, and the RV semi-amplitude K. The systemic velocity v_γ was fitted independently for each instrument (see Sect. 3.2.6). Additional free parameters were used for the two transiting planets: the planet-to-star radii ratio k, the orbital inclination i, and a common stellar density ρ_*. This represents a total of five free parameters per planet, with five more for the transit light curves. We made use of a normal prior on the stellar density ρ_* ~ 1.289 ± 0.048 ρ_⊙ that we derived from the stellar properties (see Sect. 3.1).

3.2.2 CHEOPS photometry

We discarded 69 out of 3496 (1.97%) CHEOPS data points flagged by the DRP (DRP ‘EVENT’ flag >0). Among the remaining 3427 data points, we identified 48 (1.40%) outliers by performing a 3σ-clipping visit by visit, and also flagged 121 (3.53%) points with background levels higher than 4.1 × 10⁵ electrons beyond which correlations with the flux start to appear (e.g., Fig. 2 of Deline et al. 2022).

After validating that planet c never transits during the several CHEOPS observations, we included a transit model for planet b only using the batman python package (Kreidberg 2015). Section 3.2.1 describes in details the joint modeling of planetary signals. We de-trended each visit with a Gaussian process (GP) as a function of the spacecraft roll angle using a Matérn-3/2 kernel from the celerite2 package (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). The hyper-parameter values were the same for all visits, but each visit was fitted independently. We also included a slope in our model for some visits (no. 2, 3, 6, 7, 8, 10, 12) showing a significant linear trend. We quadratically added a jitter term to the error bars of each visit in order to account for underestimation of uncertainties. This leads to a total number of 33 free parameters for the correction of CHEOPS-related systematics, with one flux mean level per visit, one jitter term per visit, a slope for seven visits, and two hyperparameters for the GP model. The de-trended CHEOPS transit of planet b obtained after the joint-fit are shown in Fig. 1 and the individual light curves are in Fig. A.1.

Fig. 1 CHEOPS phase-folded transit of HD 3167 b. Top panel: de-trended transit data points (blue) and the binned data (black) obtained from the joint fit. A sample of 100 transit light curves drawn from the posterior distribution is represented in orange. Lower panel: best-fit residuals.

3.2.3 K2 photometry

We performed a preliminary minimization fit to remove the transits of planets b and c (modeled with the batman package) and the long-term photometric variability (modeled with celerite2 GP and a Matérn-3/2 kernel) in order to identify outlying data points from the residuals. We used a 6.5σ-clipping criterion on the residuals to discard 23 outliers out of the 3448 data points (0.67%). The choice of 6.5σ was motivated as a good trade-off between efficient clipping and avoiding cutting off non-outliers in the noisiest parts of the time series.

From the residuals, we spotted by eye some significant changes in the spread of data points indicating differences in the noise level over the 80-day long observations. We identified three time ranges with the middle one having the lowest apparent noise level (see Fig. 2). We investigated the cause of this phenomenon and found that it correlates well with the frequency at which K2’s thrusters are firing to correct the pointing drift of the spacecraft. We identified precisely the times at which the noise level changes by selecting the times providing the best likelihood among several minimization fits of the light curve. Each fit was performed using the model described before (GP and transit models for planets b and c) and allocating an individual jitter term per time range. We computed the best-fit likelihood for several pairs of times and selected the time pairs with the maximimum best-fit likelihood. The final jump timings were fixed at BJD_TDB times of 2457 406.95 and 2457 436.07, respectively. By considering these three time ranges separately and using individual jitter terms, we aimed at limiting the bias induced by noisy regions due the underestimation of error bars, and at maximising the precision obtained from the low-noise middle region.

For the final fit of K2 time series, we used the same model as for the preliminary fit (batman transit models and Matérn-3/2 GP kernel). The transit models are oversampled with respect to the data cadence by a factor 30 (oversampled cadence of about 1 min), and binned down to the sampling rate of the light curve. This technique allows one to account for distortion effects due to long integration times (e.g., Kipping 2010) with strong effect during ingress and egress especially. The GP model fits for the correlated noise in the data that corresponds to both instrumental systematics and stellar activity. Gandolfi et al. (2017) pointed out the presence of the latter in the K2 data with a significant peak in the periodogram matching the stellar rotation period at ~ 24 days. We placed normal priors on both GP hyper-parameters to help convergence of the fit based on the values obtained when analysing the K2 time series alone.

The prior values are N(370, 70) ppm for the GP amlitude σ_GP, and N(10, 0.5) days for the GP correlation time scale ρ_GP, where N(µ, σ) represents a normal prior of mean µ and variance σ².

Figure 3 shows the best joint-fit de-trended K2 transits of planets b and c. The V-shaped transit of HD 3167 b is due to the long cadence of the observations that averages the sharp ingress and egress with the flat regions outside and inside the transit (e.g., Kipping 2010).

Fig. 2 Noise level time ranges in the K2 time series. Top panel: normalized K2 flux (blue points) with the best-fit model (transits + GP) obtained from minimization. Mid-panel: residuals data after removing best-fit model. Bottom panel: flag indicating when the K2 spacecraft is firing its thrusters to correct pointing drifts. In all three panels, the vertical dotted black lines show the noise jump timings providing the best likelihood.

3.2.4 Spitzer photometry

The Spitzer photometry was pre-processed prior to the joint analysis to correct for instrumental systematics (see Sect. 2.3). During the joint fit, the Spitzer observations were fitted with a transit model of planets b and c, and two additional free parameters per observation to account for the flux offset and the underestimation of uncertainties (white noise jitter term). The resulting transit light curves obtained with Spitzer are shown in Fig. 4.

3.2.5 HST photometry

We started by manually flagging one obvious outlier in the fourth orbit of the fifth transit observation. We also discarded 126 out of 722 (17.45%) points that correspond to the first orbit of every visit, and to the first point of every orbit, following standard approach (e.g., Mikal-Evans et al. 2021).

Given the periodicity of the systematic ramp every HST orbit, we decided to adopt an approach similar to the one typically used with CHEOPS, that is a GP detrending as a function of the spacecraft orbital phase. To properly determine the phase of each data point, we computed a precise orbital period for HST from the housekeeping parameters. We used the spacecraft latitude stored in the jitter files (*jit.fits files) and, for each visit, we fitted the latitude variations with a sine wave. We combined the outcome of the fits and computed a precise orbital period of min.

The GP correction as a function of the orbital phase was performed using a Matérn-3/2 kernel with a single set of hyper-parameters for all visits, but each visit was fitted individually. We let free the period of HST, while using a strong Gaussian prior N(95.23, 0.02) min based on the previously derived value.

The orbital phase of each data point was therefore computed at each iteration as a function of the HST period value and the JD_UTC time. For each visit, we also included two flux mean values (forward and backward scans) and a jitter term to account for underestimation of error bars. We found that the addition of a linear trend as a function of time was necessary for all visits, whereas the addition of a quadratic trend was required for visits 2 and 3 only. This leads to a total of 25 parameters for the correction of HST systematics in the data.

We mentioned in Sect. 2.4 that our careful analysis of the HST data lead to the discovery of serendipitous transits of planet b during visits 3 and 4. Therefore, we added a batman transit model accounting for both planets b and c and the best-fit models derived from the joint analysis are represented in Fig. 5. The data of each individual visit are shown in Fig. B.1.

Fig. 3 K2 phase-folded transits of planets b (top) and c (bottom). Blue points represent the data after detrending for any other signal. The orange curve are transit models with parameter sets randomly drawn from the posterior distribution. Binned data are represented in the top panel (transit of planet b).

Fig. 4 Spitzer phase-folded transits of planets b (top) and c (bottom). Blue and black points are the de-trended and binned data, respectively. The orange shaded area is made of samples drawn from the posterior distribution.

Fig. 5 HST phase-folded transits of planets b (top) and c (bottom). De-trended data corrected for the strong periodic instrumental systematics are shown in blue. Binned data points are represented in black. The orange curves are samples from the posterior distribution of the joint analysis.

3.2.6 Radial velocity

We computed the Generalized Lomb-Scargle (GLS – Ferraz-Mello 1981; Zechmeister & Kürster 2009) periodogram of the selected RV data set after removing an offset (median value of the time series) for each instrument, and we clearly detected the signals from the three planets (see Fig. 6). We also have three other significant peaks. Two of them seem to be induced by the stellar rotation with one peak at the rotation period (~24 days, also present in the K2 data; Gandolfi et al. 2017) and another at half of this value (see Sect. 4.1 for a detailed discussion). The last very significant peak spans a range of possible periods from 70 to 120 days, and it does not match any of the known objects in the system.

We analyzed the GLS periodograms of the different stellar indices available for the HARPS and HARPS-N data. We found significant peaks at ~24 days for the Hα and Na I lines in the HARPS-N data, and nothing around 100 days in either data sets. We also looked at possible correlations between the indices and the RV signals using the Pearson, Kendall and Spearman criteria. We detected significant correlations of the HARPS data with the FWHM (p – values < 6 × 10⁻³) and of the HARPS-N data with the Hα line (p–values < 2 × 10⁻⁶).

We first designed our RV model using three Keplerian functions for planets b, c and d, and a systemic velocity value for each of the four instruments (APF/Levy, Keck/HIRES, HARPS and HARPS-N). We also included a white-noise term (jitter) for each instrument to account for uncertainty underestimation. Based on the correlation analysis, we jointly fitted for linear functions of the FWHM and the Ha line to correct the HARPS and HARPS-N data, respectively. This modeling choice was motivated by the will to minimize the number of free parameters, even though these correlations may not be as strictly linear as we assume (e.g., Boisse et al. 2011; Santos et al. 2014; Collier Cameron et al. 2019). We ran a minimization fit and sampled the parameter space with a MCMC approach. We obtained planetary parameters fully consistent with the values from both Christiansen et al. (2017) and Gandolfi et al. (2017). We computed the periodogram of the residuals and retrieved the very significant peak at ~100 days, with a false alarm probability (FAP) smaller than 10⁻¹⁰.

We investigated the possible source of this signal by first looking at indicators of stellar activity at those periods but found no significant signatures. We also compared the periodograms of each spectrograph to search for potential discrepancy that could mean an instrumental origin for the long-period signal (see Fig. 7). We found that both HARPS-N and Keck/HIRES data have power in this regime (FAP < 10⁻⁴). We note that there is also a hint of signal in the APF/Levy time series even though it is less significant. The HARPS data covers only 128 days in total with a poor sampling over this baseline and thus the individual periodogram does not feature any peak around 100 days. The presence of power in several data sets strongly suggested that the signal was not induced by instrumental systematics and might actually have a planetary origin. To further consolidate this hypothesis, we compared the phases of the long-period signal measured by each spectrograph by looking at the RV time series phase-folded on the detected period. Figure 8 shows the best result of the joint fit performed in this work where one can visually validate the consistency of the signal phase through all instruments. In the view of these different outcomes, we rejected the stellar or instrumental origins to explain the residual signal and we interpreted it as the RV signature of a fourth planet, whose presence was suggested by Dalal et al. (2019).

We ran another fit to the data including an additional Keplerian model with a uniform prior on the period spanning a large range from 60 to 200 days. The resulting semi-amplitude of the Doppler signature produced by the new planet was significant by more than 9σ and the residual periodogram was not featuring significant peaks anymore. The comparison of the 3-planet and 4-planet models using the Bayesian and Akaike Information Criteria (BIC and AIC) clearly favored the inclusion of the new planet (ΔBIC < −57 and ΔAIC < −68). Therefore, the final RV model used in the joint fit with the photometry included a Keplerian model for each of the three known planets b, c and d, and another for the new long-period planet e. We fitted the systemic velocity and a jitter term for each instrument, and added linear correction of the HARPS and the HARPS-N data as functions of the FWHM and the Hα line, respectively. We used ten parameters in total to fit the instrument-related effects (offset, noise and decorrelation).

The extracted planetary signals fitted from the joint analysis are represented in Fig. 9. The gap visible in the RV data phase-folded on the orbital period of HD 3167 b is due to the removal of in-transit points to avoid being affected by the Rossiter-McLaughlin effect.

Fig. 6 Generalized Lomb-Scargle periodogram of the RV data. The colored triangles at the top represent several periods of interest: orbital periods of the planets b, c and d, expected rotation period of the star (~24 days), and one year. The full triangles show the main periods while the empty ones are the first three harmonics of each period. The horizontal grey lines highlight the False Alarm Probabilities (FAP) of 10−3, 10−4 and 10−5.

Fig. 7 Generalized Lomb-Scargle periodogram of the RV data for each instrument. For comparison purposes, each periodogram is normalized by the power corresponding to a False Alarm Probability (FAP) of 10−4, which is highlighted by the horizontal dotted line.

Fig. 8 Radial-velocity signals of the planet e measured by each instrument and phase-folded on the best-fit orbital period (P = 102.09 days). Binned data are shown in black.

3.2.7 Joint MCMC fit

We used the MCMC algorithm emcee (Foreman-Mackey et al. 2013) to explore the parameter space simultaneously for all the parameters defining the systematic correction of each instrument and the planetary properties from the transit light curves and the radial-velocity signals. We had a total of 120 parameters and, at each MCMC iteration, we computed the global log-probability by summing the log-probability obtained with each data set (CHEOPS, K2, Spitzer, HST, RV).

The MCMC run was initiated by estimating the best-fit parameters and their uncertainties on each data set independently. We used the posterior distribution to generate the first-guess parameter sets. The 1280 chains of our MCMC joint fit started with a burn-in phase of 145 000 iterations, and then sampled the parameter space with 260 000 steps. We kept one iteration every 2000 to reduce the effect of the chain autocorrelation. We checked the convergence of the chains by visually inspecting the trace plots and validated it based on the Gelman-Rubin criterion (Gelman & Rubin 1992).

Fig. 9 Phase-folded radial-velocity data of the four planets of the HD 3167 system. The color code of the data points highlights the instrument with which they were observed (APF, HARPS, HARPS-N, HIRES). The black points show the binned data. The multiple black curves are samples randomly drawn from the posterior distribution.

Fig. 10 Probability map of the spot latitude attributed to the K2 photometry and RV modulations as a function of the relative differential rotation rate. The colorscale represents the χ2 probability and the dotted lines highlight the 1, 2 and 3σ confidence intervals. The vertical dashed lines shows the 1σ confidence interval of the expected α value derived from Eq. (2) of Balona & Abedigamba (2016).

4 Revision of the system properties

4.1 Star

The Rossiter-McLaughlin analysis performed by Bourrier et al. (2021) revealed variations in the contrast of the stellar lines occulted by HD 3167b and c along their respective transit chords. The authors could explain these variations through a latitudinal dependence of the stellar line contrast, which allowed them to constrain the inclination of the star with respect to the line of sight ( or , the two configurations being degenerate) and thus the true equatorial velocity ( km s⁻¹). Assuming that this scenario is correct, we combined our stellar radius (Sect. 3.1) with the stellar inclination and projected rotational velocity from Bourrier et al. (2021) to update the true equatorial period ( d).

In addition to the planetary signals, we observe consistent peaks in the periodograms of the K2 photometry and RV data at about ~24 days, which likely trace spots on the rotating stellar surface. The difference between this period, and that derived independently from Rossiter-McLaughlin analysis for P_eq, could indicate that the star is rotating differentially. During their lifetime, spots would spend on average more time in a region located at higher latitudes, rotating more slowly than the equator. Under this hypothesis we can estimate the spot location by assuming a solar-like law for the stellar differential rotation: (1)

with αr the relative differential rotation rate between equator and pole, and θ the stellar latitude. We computed the value of P(θ) from a Gaussian fit to the periodogram peak in the K2 and RV data sets, yielding P_K2 = 23.4 ± 2.2 days and P_RV = 24.1 ± 1.2 days. The close agreement between these periodic signals from different datasets give us confidence that they arise from stellar modulation, rather than instrumental variability. We then built a χ² map of the stellar latitude as a function of α, comparing the theoretical P_eq/P(θ) ratio with the measured values (Fig. 10). The relative differential rotation rate for HD 3167 can be estimated independently from Eq. (2) of Balona & Abedigamba (2016), which is based on photometric modulations in a wide sample of Kepler stars. Using T_eff = 5300 ± 73 K and Ω_eq = 2π/P_eq, we derive α = 0.179 ± 0.028. Within 2σ, this value is consistent with the measurements of P_eq and P(θ) for stellar latitudes ≳50deg. Spots on HD 3167 would thus be located closer to its poles than on the Sun, where spots appear at latitudes of ~35° at the beginning of a new cycle and converge toward the equator over 11 yr.

4.2 Planets

We derived the parameter values and uncertainties from the posterior distribution of the MCMC run. We computed the median and the 68.3% confidence interval for each of the fitted parameters (see Table 2). We combined the MCMC chains to calculate the values and uncertainties on a series of useful parameters that were not fitted directly (Table 3). Some of these derived parameters were obtained using the stellar mass or radius, and we propagated the uncertainties on these parameters by drawing random values from normal distributions based on their estimated values: M_* = 0.852 ± 0.026 M_⊙ and R_* = 0.871 ± 0.006 R_* (Sect. 3.1).

All the fitted and derived parameter values are consistent with the ones reported by Christiansen et al. (2017) and Gandolfi et al. (2017). The inclusion of the CHEOPS, HST and Spitzer data sets allows us to improve significantly the precision on the orbital periods of planets b, c, and d by factors of ~40, > 50 and ~ 17, respectively. We obtain a better precision on HD 3167 b and c planet-to-star radii ratios in the K2 passband analyzed in Christiansen et al. (2017) and Gandolfi et al. (2017), and we improve by more than a factor two the absolute planetary size thanks to the smaller uncertainty on the stellar radius. We also reduce the errors on the absolute and minimum masses of b, c and d thanks to the improved RVs reduction and additional datapoints. These improvements on the planets mass and radius lead to an overall reduction of the uncertainty on the bulk densities of HD 3167 b and c by more than a factor three (see Fig. 11).

We detect the fourth planet HD 3167 e with a semi-amplitude significance > 8σ and a minimum mass of M_⊕. This planet was hinted in the RV data previously available to Dalal et al. (2019), as a 0.03 M_J outer companion with an orbital period of 78 days that could explain the peculiar orbital architecture of the system. The order of magnitude of both the mass and the period we derive matches well their original estimates. We note that the orbital period of this new planet has a very peculiar marginalized posterior distribution. Indeed, the uncertainty of about 0.5 days listed in Table 2 is dominated by the main mode of the distribution. However, the MCMC solution does also explore other possible orbital periods that are less likely but nevertheless span a large range from 79 to 125 days (see Fig. 12). We explain this distribution and its invariance with respect to the time of inferior conjunction T_0,e by the fact that the data are strongly unevenly sampled. Most of the data points (94%) were taken during the first year of observation (over about 200 days) and the remaining 6% are HARPS-N data spread over four years. Therefore, T_0,e is strongly constrained by the bulk distribution of the first year that lasts less than two periods of planet e. The multi-modal distribution of P_e reflects the uneven sampling of the RV time series by highlighting the periods that best match the scattered data. The other orbital parameters of HD 3167 e are well defined and not correlated with P_e.

The orbits of the four planets are fully consistent with circular configurations with upper limits (99.73% confidence) on their eccentricities of e_b < 0.11, e_c < 0.15, e_d < 0.45 and e_e < 0.61 for planets b, c, d and e, respectively.

We provide average transit depths and planetary radii for planets b and c, obtained from the merged distributions over the four available instrumental passbands (CHEOPS, K2, HST/WFC3/G141, Spitzer/IRAC/Ch2). To further characterize the system, we allowed the radii of the two transiting planets to vary independently in those passbands. We measure consistent radii for the inner planet b, which is expected from an USP planet unable to retain any volatile atmosphere. However, we note a significant difference (>3.5σ) between the radius obtained for planet c in the HST/WFC3/G141 (λ ~ 1.4 µm) and the Spitzer/IRAC/Ch2 (λ ~ 4.5 µm) passbands (see Fig. 13). This difference could arise from broadband variations in the optical depth of the planet atmosphere linked to its chemical composition and physical structure (Guilluy et al. 2021; Mikal-Evans et al. 2021).

4.3 Stellar companions

In order to check for stellar companions lying within the environment of HD 3167, high angular resolution optical speckle interferometric imaging was performed. HD 3167 was observed on 2021 June 28 UT using the ‘Alopeke speckle instrument on Gemini North (Scott et al. 2021). ‘Alopeke provides simultaneous speckle imaging in two narrow bands (562 and 832 nm) with output data products including a reconstructed image and robust contrast limits on companion detections (Howell et al. 2011, 2016). The night had clear skies and good seeing (<1.0 arcsec) during the observations. As shown in Fig. 14, we detect no stellar companions which are brighter than two delta-magnitudes within 0.1″ and no companions brighter than five to 8.5 magnitudes within the angular separation limits of 0.1″ – 1.2″. Using a distance of d = 47 pc for HD 3167, these angular and luminosity limits on stellar companions correspond to main sequence stellar types of K6V (at 0.94 au) and M2.5V to M4.5V between 4.7 and 56.4 au.

5 Planetary internal structures

Using the derived stellar and planetary properties (in particular the mass and average transit depths reported in Table 3), we computed the internal structure of both transiting planets using a Bayesian analysis, and following the method described in Leleu et al. (2021). We recall here the two main elements in this method: the assumed priors and the forward model.

Our forward model computes the radius of planets as a function of some hidden parameters: mass of the solid Fe/S core, fraction of Fe in the core, mass of the silicate mantle and its composition (Si, Mg and Fe molar ratios), mass of the water layer, mass of the gas envelope (composed in this model of pure H/He), equilibrium temperature of the planet, and age (which is supposed to be the same as the age of the star). We assume in our model that the Si/Mg/Fe ratio in the bulk planet is the same as in the star (Dorn et al. 2015, Thiabaud et al. 2015). Note that recently, Adibekyan et al. (2021) have shown that these ratio are indeed correlated but may not follow a 1-to-1 correlation. Including this in the model is the subject of future work. Regarding the priors, the core, mantle and water mass fraction (relative to the non-gas part) follow a uniform prior (subject to the constraint that they add up to one), whereas the mass fraction of the H/He layer follows a prior which is uniform in log. We finally note that the gaseous (H/He) part of the planet does not influence, in our model, the ‘non-gas’ part of the planet (core, mantle and water layer). This means that the innermost layers of the planet are not modified by the potential compression and thermal isolation effect from the gas envelope.

Figure 15 shows the resulting internal structure of both planets presented as corner plots and summarized in Table 4. Planet c hosts a substantial gaseous envelope, weighing a little less than 0.2 M_⊕, whereas its fraction of water is unconstrained. We emphasize that this result depends on the assumed priors. In particular, the resulting planetary model would be more gas rich and less water-rich if the H/He layer followed a uniform prior. One of our main findings is that planet b mass and radius seem to be inconsistent with a pure iron-core and silicate-mantle structure whose composition would reflect the Fe/Si/Mg ratio in the star. Indeed, the density of the planet is smaller than what would be expected for such a model. Since our model assumes the inner layers of the planet to be unaffected by the influence of the gas envelope, for planet b we underestimate the temperature of the ‘non-gas’ part of the planet. If we increase the temperature of the core and mantle layers in our model, we do observe an increase in radius of up to 2% for pure iron-core and silicate-mantle structures matching the Fe/Si/Mg ratio of the star. However, this effect alone is not enough to explain the observed radius of the planet. We hence conclude that a light element must be present in the planet.

Our fit converges toward a negligible mass fraction of gas, which is expected considering that the intense irradiation of this USP would lead a H/He atmosphere to be lost extremely fast. However, the mass fraction of water for our model of HD 3167 b is quite well constrained and non zero. With an equilibrium temperature in excess of 1600 K (Table 3) any water layer would be made of steam, which has been shown to be much more resilient to atmospheric loss (e.g., Lopez et al. 2012). It should be kept in mind that our model assumes a fully differentiated planet. It is possible that water is mixed with a magma ocean covering HD 3167b, in which case its actual mass fraction of water would be reduced compared to the one we derive (see Dorn & Lichtenberg 2021). More detailed internal structure model accounting for this mixing, and for the existence of a dust- and metal-rich envelope, are required to better constrain the true nature of this planet.

Fig. 11 Mass-radius diagram. The grey point represent all the exo-planets from the Extrasolar Planets Encyclopaedia (http://www.exoplanet.eu/) as of 1 March 2022, with masses and radii known with a precision better than 30%. HD 3167 b and c are represented in the different instrument passbands analyzed in this work (CHEOPS. K2, HST/WFC3/G141 and Spizer/IRAC/Ch2). The black dash-dotted lines indicate two iso-density profiles matching the planets b and c. The colored solid and dashed lines are indicative chemical compositions as computed by Zeng et al. (2019).

Fig. 12 Corner plot of the posterior distribution of the planetary parameters of HD 3167 e. The orbital period shows a very peculiar distribution that is not significantly correlated to any other parameter.

Fig. 13 Posterior distributions of the planet-to-star radii ratios of planets b and c measured in the four instrument passbands.

Fig. 14 5σ contrast curves of HD 3167 from Gemini North/’Alopeke at 562nm (blue) and 832 nm (red). The insert shows our reconstructed 832 nm speckle image.

6 Dynamical evolution

With planet b aligned with its host star (Bourrier et al. 2021) and the distant planet c on a polar transit (Dalal et al. 2019), the HD 3167 system is particularly rich and interesting for dynamical studies. Its planets have wide orbital separations and are far from mean motion resonances, so that their dynamics is fully secular. Moreover the age of the system suggests that its orbital configuration is dynamically stable.

The dynamical analysis by Dalal et al. (2019) suggested that planet b could have stayed aligned with the equatorial plane of the star, which has since been confirmed by Bourrier et al. (2021). They further showed that planets c and d have a low mutual inclination, and that the three planets known at that time could not, by themselves, have caused the polar orbit of planet c. Dalal et al. (2019) thus proposed that the orbits of planets c and d could have been tilted with respect to the star due to a massive outer companion, whose existence we have confirmed in the present study (Sect. 4.2). It is thus natural to investigate whether this planet e could indeed explain the polar orbit of planet c.

To gain some insights onto the dynamics of the system, we consider an analytical framework describing the precession of the orbits (Boué & Laskar 2006). Following Boué & Fabrycky (2014), we compute the characteristic frequencies v^k/j that represents the relative influence of the body k over the direction of the angular momentum of body j. We refer to Dalal et al. (2019, Sect. 5.2.) for a more precise description of the model used here. To summarize, if v^k/j ≪ v^j/k, the angular momentum direction of j is almost constant while the angular momentum of k precesses around.

We study the dynamics of the system at two different epochs. First, in its current configuration, with a stellar rotation period of about 17 d. Second, right after the system’s formation, when the star was rotating fast and its spin had a stronger influence onto planet b. The different characteristic frequencies⁶ in these two configurations are summarized in Table 5. While the precession frequencies ruling the planet interactions are the same in both settings, there is a significant change for the interaction between the star’s spin and the planets as we can see on the frequencies v^b/S and v^S/b. The change is not only due to the shorter rotation period of the star early on, but also because the second fluid Love number k₂ can be significantly larger for a fast rotating star (Becker et al. 2020). We adopt a value of k₂ = 0.18 for the fast rotating star, an order of magnitude larger than the expected value for HD 3167 today.

Fig. 15 Corner plot showing the results on the interior composition models of HD 3167b (left) and HD 3167c, (right). The vertical dashed lines and the ‘error bars’ given at the top of each column represent the 5 and 95% percentiles. The internal structure parameters are the mass fraction of the core and the water layer with respect to the solid planet, the molar fraction of Si and Mg in the mantle, the molar fraction of Fe in the inner core and the logarithm of the gas mass in Earth masses.

Fig. 16 Inclination evolution in the current system configuration, assuming the planets c, d and e are close to coplanar.

Fig. 17 Evolution of the inclinations and eccentricities assuming the current star properties and a planet e initially misaligned with the system.

6.1 System stability

In the present system configuration we have v^b/d, v^d/b ≫ v^S/b ≫ v^b/S, which indicates that planet b’s orbit precesses around the angular momentum of the outer planets and that the star plays a negligible role. Moreover, the stellar spin is dynamically unaffected by the planets. We confirm this hypothesis by running an N-body simulation using the integrator WHFast (Rein & Tamayo 2015) from the library Rebound (Rein & Liu 2012). We include relativistic corrections as well as the influence of the stellar J₂ using the library Reboundx (Tamayo et al. 2019). The stellar spin is fixed and along the z-axis. For this particular simulation, the initial orbits are assumed circular, the planets c, d and e are assumed coplanar. We use an initial condition compatible with the 3D configuration determined by Bourrier et al. (2021), i_b = 30°, Ω_b = 100° and i_d,c,e = 100°, Ω_d,c,e = 0°. As a result, the mutual inclination between planet b and the rest of the system is i_bc = 103°. This approximate initialization is sufficient because we only want to illustrate the typical dynamics at play. In reality, there is likely a non zero inclination between planet d, c, and e since planets d and e are not transiting. We integrate the system over 100 kyr and plot on Fig. 16 the planetary inclinations with respect to the star, as well as the mutual inclination between planet b and c. During this short integration, we observe no evolution of the eccentricities. We see that the mutual inclination between b and c, as well as the inclination of the outer planets are almost constant while planet b orbital plane precesses around the orbital plane of the outer planets. We conclude that, while planet b is not strongly coupled with the star today, a primordial misalignment of the outer planets with respect to the star and planet b leads to a stable configuration, compatible with the observations. However, early in the history of the system, we have v^S/b ≫ v^b/d, v^d/b, v^b/S which means that planet b could have been coupled with the star instead of the outer planets. If planets c and d gained their large obliquities early on, planet b would have stayed close to the stellar equator.

6.2 Spin-orbit misalignment by planet e

We now investigate whether planet e could be the cause of the polar orbits of planets c and d. The characteristic frequency v^e/c is not negligible, which suggests that planet e could tilt planets c and d if it was originally inclined with respect to the star. Boué & Fabrycky (2014) have determined that an external companion can tilt a planetary system as a whole if the coefficient (2)

where . However, planet e is too close to c to tilt the system without triggering Kozai-Lidov oscillations for planet d and c (Kozai 1962; Lidov 1962). Indeed, we have β_KL = 0.88 ± 0.13. As a result, a large mutual inclination of planet e with respect to c excites the orbital eccentricities, eventually leading to the system destruction. We run a numerical simulation starting with the inner planet b, d and c on coplanar, circular orbits within the star equatorial plane and a planet e on an orbit tilted by 80°. The simulation lasts 200 kyr and we plot on Fig. 17 the eccentricities and inclinations as a function of time. As expected planets c, d and e enter Kozai-Lidov oscillations and the eccentricities grow to values close to 0.5, which is excluded by the observations (at 2σ, e < 0.3 for planets d and e and e < 0.08 for c) and is enough to trigger the dynamical instability of the system. Moreover, in that scenario planet b would remain in the plane of planets d and c, which confirms that the outer system had to get misaligned early-on when the coupling between planet b and the star was stronger. Planet e thus cannot explain the polar orbit of planet c.

6.3 Spin-orbit misalignment by an outer companion

While a misaligned planet e leads to the instability of the system, we explore whether a more distant companion could explain the present configuration. Such an hypothetical companion should lie in a range of masses and semi-major axes that verify the Kozai-Lidov condition (3)

while being able to tilt the system as a whole. We plot on Fig. 18 the ratio of the precession frequency v^co./pl.∙/v^S/pl. that represents the relative influence of a companion onto the outer planets with respect to the influence of the star onto the system as a function of the semi-minor axis and the mass of a potential companion. We also plot the levels of β_KL,ce = 1 and β_KL,ce = 0.1. A companion can tilt the outer planets if condition (Eq. (3)) is met and precession frequencies verify v^co./pl.≫ v^S/pl.. Additionally, the companion has a significant influence on the inclination of the system if its angular momentum is larger than the system’s angular momentum. Otherwise, the system orbital plane remains unchanged and the companion orbital plane precesses around. We plot this theoretical constraint in Fig. 18.

We further included in Fig. 18 the constraints derived from our direct imaging measurements (Sect. 4.3). We converted luminosity differences with the K0V type star HD3167 into spectral types for various separations, and then used Table 5 from Pecaut & Mamajek (2013) to assign mean masses to these spectral types. We also assumed circular orbits to estimate the semi-minor axes (which is the most conservative assumption). While the masses adopted for a given dwarf subtype are tentative, it provides an approximate upper limit on the companion mass. The direct imaging constraints impose that it orbits within ~30au from the star and cannot be more massive than about 0.1 solar mass. Finally we included in Fig. 18 the stringent constraints from our RV dataset. The constraints change little with the unknown orbital inclination of the companion unless it is seen nearly pole-on. The RV constraints rule out most, if not all, the configurations where an outer companion could lead to a significant misalignment of the system. We conclude that the polar orbits of the outer planets are most likely due to an early misalignment during the system formation that did not rely on a companion still present in the system. Mechanisms that do not rely on binary companion or secular interactions have been proposed such as magnetic coupling between the young star and the disk (Lai et al. 2011; Foucart & Lai 2011; Romanova et al. 2021).

Fig. 18 Range of masses and semi-minor axis that allows a companion to tilt the outer planets of the system while preserving its mutual inclinations. The system is not destroyed for companions below the line βKL,ce = 1 and tilting is possible for a frequency ratio larger than 1. The companion has a significant influence on the system if its angular momentum is larger than the system’s angular momentum. Otherwise, the companion’s orbit precesses around the system angular momentum and does not induce large obliquity for the system. The green region is the 5σ limit set by direct imaging constraints. The purple and blue curve are the 5σ limit set by RVs, assuming a circular orbit and inclinations of 90 and 35°.

7 Atmospheric evolution

To constrain the atmospheric evolution of the planets in the HD 3167 system and the stellar rotation history of the host star, we employed the tool Planetary Atmospheres and Stellar RoTation RAtes (PASTA; Bonfanti et al. 2021b). PASTA uses the measured system parameters and the present-day atmospheric mass fractions determined by the internal structure modeling (see Sect. 5) to return posterior probability distributions for the initial atmospheric mass-fraction of each planet, further constraining the history of the stellar rotation rate. Because of the need of an estimate of the present-day atmospheric mass fraction, or at least of a radius measurement, to constrain the evolution of the planetary atmosphere, this tool in its full capability can only be employed on the two transiting planets HD 3167 b and c. Since the present-day stellar rotation rate is not well defined, we employed a uniform prior ranging between 15 and 20 days.

7.1 Transiting planets HD 3167 b and c

Planet HD 3167 b orbits very close to its host star, resulting in it having being subject to large amounts of X-ray and extreme ultraviolet (XUV) irradiation, particularly in the early phases of its evolution. PASTA predicts that this planet has lost its primary H/He-dominated atmosphere at some point in the past, and thus the code is unable to constrain the initial atmospheric mass-fraction, resulting in a uniform posterior distribution (left panel of Fig. 19).

For planet HD 3167 c, PASTA prefers evolutionary tracks for which atmospheric mass loss did not play a significant role. This is represented by the posterior distribution of the initial atmospheric mass fraction that peaks around the present-day value (right panel of Fig. 19). However, the results also indicate that evolutionary tracks characterized by significant mass loss, though less likely, are not completely excluded. In case the host star was a particularly fast rotator when it was young, it would indeed have emitted a significant amount of XUV radiation (Sanz-Forcada et al. 2011). Therefore, we explored this possibility more thoroughly. Unfortunately, the characteristics of HD 3167 c, which is the only transiting planet in the system still holding its primordial H/He-dominated atmosphere, do not enable PASTA to constrain the rotation history of the star. This is represented by a rather flat posterior distribution of the stellar rotation rate after 150 Myr that we use as proxy to illustrate the evolution of the stellar rotation rate (Fig. 20). The rotation rate distribution of stars with a mass comparable to that of HD 3167 and member of open clusters with ages of about 150 Myr is bimodal (e.g., Johnstone et al. 2015): one peak represents fast rotators at a rotation rate close to one day, while the other peak represents moderate rotators at roughly six days. The observed distribution is shown in Fig. 20. We have performed additional runs with PASTA, imposing priors on the stellar rotation rate at 150 Myr corresponding to Gaussian fits to each of the two peaks of the distribution. Assuming the host star was a fast rotator when it was young, the relative occurrence of evolutionary tracks presenting some significant atmospheric loss increases compared to when not constraining the stellar rotation type at all. However, we still obtain a posterior distribution of the initial atmospheric mass fraction which peaks close to the present-day atmospheric mass fraction, indicating that most likely atmospheric loss has not played a major role in the evolution of this planet independently of the evolutionary history of the stellar rotation rate.

7.2 Non-transiting planets HD3167 d and e

Since HD 3167 d and e are not transiting, it is not possible to estimate their current atmospheric mass fraction, and thus it is not possible to use PASTA to infer the initial atmospheric mass fraction. To estimate the current atmospheric content of both planets, we start from assuming that the planets accreted a primordial H/He-dominated atmosphere of (Mordasini 2020) (4)

where M_env,0 is the envelope-mass, M_c the core-mass, and a the planetary orbital separation. We used PASTA to compute the atmospheric evolution of HD 3167 d and e starting it with the atmospheric mass fraction given by Eq. (4). The simulations also require an estimate of the evolution of the rotation rate of the host star. Since PASTA has been unable to constrain the stellar rotation history using planets b and c, we further assumed a value of the rotation rate of the host star at 150 Myr of 5.44 days. This value corresponds to the mean of the distribution of stellar rotation rates of young open cluster stars with mass comparable to HD 3167 (Fig. 20).

As HD 3167 d orbits quite close to its host star and has a low mutual inclination with the transiting HD 3167c (Sect. 6), we assumed the measured lower mass limit M_d sin(i) to be a good approximation of the core mass. Through Eq. (4), we then estimate an initial atmospheric mass fraction of 0.029. PASTA’s evolution simulation predicts for this planet to have lost all of its primordial H/He-dominated envelope via photo-evaporation. The atmospheres of planets b and d therefore seem to have had very similar evolutionary paths.

HD 3167 e on the contrary orbits significantly further away than HD 3167 c, for which we already determined that hydrodynamic mass loss was only important if the star would have been a very fast rotator. Therefore, we assumed for the measured lower mass limit M_e sin(i) to be a good approximation of the initial total mass of the planet after the dispersion of the protoplanetary nebula. With these assumptions, Eq. (4) leads to an initial atmospheric mass fraction of 0.147. Using the estimates of the initial total mass and atmospheric mass fraction, estimates on the core-mass and envelope-mass after the dispersion of the nebula can be provided. Assuming an earth-like density for the core, this results in an estimate of the average density of the planet at the beginning of its atmospheric evolution. This average density then changes due to atmospheric loss as the evolution progresses. As expected, PASTA’s evolution simulation predicts no significant mass loss for this planet, resulting in a Saturn-like density of about 0.67 g cm⁻³. This value is however heavily dependent on Eq. (4), as it estimates the initial atmospheric mass fraction and therefore the initial core- and envelope-mass.

Fig. 19 Posterior distributions of the initial atmospheric mass fractions for planets HD 3167 b and c derived by PASTA. The light-blue line represents the distribution of the estimated present-day atmospheric mass-fraction. The orange horizontal lines indicate the uninformative prior distributions.

Fig. 20 Posterior distribution (blue line) of the stellar rotation rate of HD 3167 after 150 Myrs derived by PASTA. The purple area represents the highest posterior density (HPD) interval of the distribution. The black line represents the distribution of the stellar rotation rate of young open cluster stars with mass comparable to that of HD 3167 based on the collection of data provided by Johnstone et al. (2015).

7.3 Comparison with previous studies

Kubyshkina et al. (2019) applied an earlier version of PASTA and Bonfanti et al. (2021b) the same version of PASTA used here on the HD 3167 system considering the system parameters available at the time. Bonfanti et al. (2021b) focused on the two transiting planets b and c, while Kubyshkina et al. (2019) did also investigate the atmospheric evolution of the non-transiting planet d. For the two close-in planets b and d the two previous studies agree with our conclusion that the planets should have lost all of their primary H/He-dominated envelopes.

Both Kubyshkina et al. (2019) and Bonfanti et al. (2021b) ran their fits considering the planetary radii, instead of the current atmospheric mass fractions, which led them to assume slightly larger atmospheric mass fractions compared to what was used here. Furthermore, both studies used significantly longer present-day stellar rotation rates, with Kubyshkina et al. (2019) using a prior peaking at roughly 25 days, while Bonfanti et al. (2021b) used an even larger peak value of about 50 days, based on gyrochronological considerations. Bonfanti et al. (2021b) obtained that the star was likely to be a slow rotator, which might be related to the large value of the current stellar rotation period they considered. Along the lines of our results, they concluded that planet c has most likely retained most of its primary H/He-dominated atmosphere and, therefore, has not undergone significant mass loss. In contrast, Kubyshkina et al. (2019) concluded that the young star was a moderate-to-fast rotator, in agreement with our hypothesis.

8 Conclusions

We performed a joint analysis of transit photometry (CHEOPS, K2, HST/WFC3, Spitzer/IRAC) and radial velocimetry (HARPS-N, HARPS, APF/Levy, Keck/HIRES) to refine the bulk and orbital properties of the planets orbiting the bright and nearby star HD 3167. New CHEOPS photometry and RV measurements were added to the published datasets, which were re-analyzed using improved techniques.

We first revised the stellar age, mass and radius. The discrepancy found between the activity signal measured in the K2/RV data and the equatorial period derived from RM analysis of HD 3167b and c can be tentatively attributed to stellar differential rotation (at a rate of ~18%) and spots located at higher latitudes than 50°.

We confirmed the RV drift measured by Dalal et al. (2019) as HD 3167e, a fourth non-transiting planet with minimum mass of ~10 M_⊕ and a 102 day period. This discovery sheds a new light on the peculiar orbital architecture of the HD 3167 system. The present-day system is dynamically stable, with the orbits of its four planets consistent with circular configurations, but HD 3167b orbits close to the stellar equatorial plane while HD 3167d and c are on nearly coplanar, polar orbits. Using an analytical approach to investigate the secular dynamical history of the system, we showed that the tilting of the outer planets must have happened early in the system history or planet b, initially coupled with the star, would have become coupled with planets d-c and followed their misalignment. Yet planet e cannot explain the polar orbits of planets d-c without destroying the system through Kozai-Lidov oscillations, which further implies that these three planets have low mutual inclinations. We finally explored the possibility that an unknown, more distant companion tilted the outer planetary system. Our analytical estimates, combined with constraints from velocimetry and direct imaging data (Gemini North/’Alopeke), however rule out this possibility unless the companion is a massive sub-solar body with a highly inclined orbit, or a star in the birth cluster that later unbound from the system.

Our revision of the system properties improves the precision on the transiting planets radii by a factor two. When combined with the refined planet masses, this allows us to reduce the uncertainty on their density by a factor three, which is of particular interest to our understanding of their nature. Our internal structure retrievals show that the low density of HD 3167b cannot be explained by a pure iron-core and silicate-mantle, suggesting that the planet contains a substantial fraction of lighter elements. It could be water, mixed with a magma ocean and/or in a steam atmosphere resilient to evaporation, or it could be a more exotic envelope made of dust and metals. In contrast we find that HD 3167c is a mini-Neptune hosting a substantial volatile envelope, with a gaseous mass fraction of ~0.2 M_⊕ or larger if the planet is water-poor. The different passbands used for transit observations allowed us to search for broadband spectral variations in these two planet radii. We measure consistent sizes for HD 3167b, as expected from the absence of a volatile envelope. In contrast we measure significant spectral variations in the size of HD 3167c, with a smaller radius in the infrared. These results strengthen the interest and amenability of HD 3167c for follow-up observations at all wavelengths, to better determine its atmospheric structure and catalog its chemical content. Transit follow-up at high spectral resolution in the ultraviolet-optical domains, or phase curve measurements in the infrared, would also be of interest to disentangle the possible scenarios for the interior of HD 3167b. We emphasize that our data analysis improves the precision on the planetary orbital periods by more than one order of magnitude, which will greatly help future transit follow-up.

Finally, we use atmospheric simulations to bring additional insight into the history and nature of the HD 3167 planets. Due to its strong irradiation, HD 3167b lost any primordial volatile envelope shortly after its formation. In contrast, we find that atmospheric loss did not play a significant role in the evolution of HD 3167c, regardless of the stellar evolutionary history, so that its present-day atmosphere may still trace its primordial composition. With reasonable assumptions on the current atmospheric mass fraction of HD 3167d-c, we further find that planet d likely lost all of its atmosphere through photo-evaporation while planet e was unaffected and retains a substantial gaseous envelope.

To summarize, our revised picture of the HD 3167 system (Fig. 21) consists in:

• HD 3167: an old K-type star, initially a moderate-to-fast rotator

• HD 3167b: a transiting ultra-short period planet with a heavyweight envelope, initially coupled with the star and thus still orbiting near the stellar equatorial plane

• HD 3167d: a non-transiting super-Earth with no gaseous envelope, which followed a similar atmospheric evolution as planet b but a similar dynamical evolution as planet c,e

• HD 3167c: a massive transiting mini-Neptune, which likely kept its primordial envelope and was tilted its present polar orbit early in the system history

• HD 3167e: a non-transiting planet that likely followed the same atmospheric and dynamical evolution as planet c, which implies that it orbits in a nearby plane and is similar in mass.

In-depth characterization of the HD 3167 system and comparison with other multi-planet systems will shed more light on its origins and evolution.

Fig. 21 Mass-period (top panel) and mass-radius (bottom) diagrams of the exoplanet population in the Earth-Neptune range. HD 3167 planets are shown as green disks (HD 3167d and e are positioned at their measured minimum mass).

Acknowledgements

We thank the referee for their appreciative and constructive review. We warmly thank Thibault Kuntzer for his early analysis of the HD3167 photometry, and Tom Mikal-Evans for sharing the HST/WFC3 broadband photometry reduced in Mikal-Evans et al. (2021). This work has been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). D. E. acknowledges financial support from the Swiss National Science Foundation for project 200021_200726. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project SPICE DUNE, grant agreement No 947634; project Four Aces, grant agreement 724427; project SCORE, grant agreement No 851555). CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. A. C. C. acknowledges support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. S. G. S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). A. C. C. and T. G. W. acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. Y. A., M. J. H. and J. A. E. acknowledge the support of the Swiss National Fund under grant 200020_172746. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia Maria de Maeztu-Centro de Astrobiología (INTA-CSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract No. 4000124370. S. C. C. B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004. X. B., S. C., D. G., M. F. and J. L. acknowledge their role as ESA-appointed CHEOPS science team members. A. Br. was supported by the SNSA. This project was supported by the CNES. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. L. D. is an F. R. S.-FNRS Postdoctoral Researcher. This work was supported by FCT - Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 - Programa Operacional Competitividade e Internacionalizacão by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, PTDC/FIS-AST/28987/2017 & POCI-01-0145-FEDER-028987, O. D. S. D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. B. -O. D. acknowledges support from the Swiss National Science Foundation (PP00P2-190080). R. D. H. is funded by the UK Science and Technology Facilities Council (STFC)’s Ernest Rutherford Fellowship (grant number ST/V004735/1). M. F. and C. M. P. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18). D. G. gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 “Gaseous or rocky? Unveiling the nature of small worlds”. M. G. is an F.R.S.-FNRS Senior Research Associate. S. H. gratefully acknowledges CNES funding through the grant 837319. K. G. I. is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. This work was granted access to the HPC resources of MesoPSL financed by the Region Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervized by the Agence Nationale pour la Recherche. M. L. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. PM acknowledges support from STFC research grant number ST/M001040/1. GSc, GPi, IPa, LBo, VNa and RRa acknowledge the funding support from Italian Space Agency (ASI) regulated by “Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C”. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). I. R. I. acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B- C33, as well as the support of the Generalitat de Catalunya/CERCA programme. Gy. M. Sz. acknowledges the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a PRODEX Institute Agreement between the ELTE Eötvös Loránd University and the European Space Agency (ESA-D/SCI-LE-2021-0025), the Lendület LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. V.V.G. is an F.R.S-FNRS Research Associate. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Some of the observations in the paper made use of the High-Resolution Imaging instrument ‘Alopeke obtained under Gemini LLP Proposal Number: GN/S-2021A-LP-105. ‘Alopeke was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. Alopeke was mounted on the Gemini North (and/or South) telescope of the international Gemini Observatory, a program of NSF’s OIR Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).

Appendix A CHEOPS observations

Fig. A.1 CHEOPS observations of HD 3167 b. The raw data extracted by the DRP is shown in blue and the best-fit model is represented by the orange line. The transit model is shown in green with an offset for better visibility.

Appendix B HST observations

Fig. B.1 HST observations of HD 3167 c. The left column shows the raw data in blue with the best model in orange for the 5 visits. The right column represents the detrended data with the transit model in orange. Visit 3 and 4 are affected by a serendipitous transit of HD 3167 b.

References

All Tables

All Figures

	Fig. 1 CHEOPS phase-folded transit of HD 3167 b. Top panel: de-trended transit data points (blue) and the binned data (black) obtained from the joint fit. A sample of 100 transit light curves drawn from the posterior distribution is represented in orange. Lower panel: best-fit residuals.
In the text

Fig. 2 Noise level time ranges in the K2 time series. Top panel: normalized K2 flux (blue points) with the best-fit model (transits + GP) obtained from minimization. Mid-panel: residuals data after removing best-fit model. Bottom panel: flag indicating when the K2 spacecraft is firing its thrusters to correct pointing drifts. In all three panels, the vertical dotted black lines show the noise jump timings providing the best likelihood.

In the text

	Fig. 3 K2 phase-folded transits of planets b (top) and c (bottom). Blue points represent the data after detrending for any other signal. The orange curve are transit models with parameter sets randomly drawn from the posterior distribution. Binned data are represented in the top panel (transit of planet b).
In the text

	Fig. 4 Spitzer phase-folded transits of planets b (top) and c (bottom). Blue and black points are the de-trended and binned data, respectively. The orange shaded area is made of samples drawn from the posterior distribution.
In the text

	Fig. 5 HST phase-folded transits of planets b (top) and c (bottom). De-trended data corrected for the strong periodic instrumental systematics are shown in blue. Binned data points are represented in black. The orange curves are samples from the posterior distribution of the joint analysis.
In the text

Fig. 6 Generalized Lomb-Scargle periodogram of the RV data. The colored triangles at the top represent several periods of interest: orbital periods of the planets b, c and d, expected rotation period of the star (~24 days), and one year. The full triangles show the main periods while the empty ones are the first three harmonics of each period. The horizontal grey lines highlight the False Alarm Probabilities (FAP) of 10−3, 10−4 and 10−5.

In the text

	Fig. 7 Generalized Lomb-Scargle periodogram of the RV data for each instrument. For comparison purposes, each periodogram is normalized by the power corresponding to a False Alarm Probability (FAP) of 10−4, which is highlighted by the horizontal dotted line.
In the text

	Fig. 8 Radial-velocity signals of the planet e measured by each instrument and phase-folded on the best-fit orbital period (P = 102.09 days). Binned data are shown in black.
In the text

	Fig. 9 Phase-folded radial-velocity data of the four planets of the HD 3167 system. The color code of the data points highlights the instrument with which they were observed (APF, HARPS, HARPS-N, HIRES). The black points show the binned data. The multiple black curves are samples randomly drawn from the posterior distribution.
In the text

Fig. 10 Probability map of the spot latitude attributed to the K2 photometry and RV modulations as a function of the relative differential rotation rate. The colorscale represents the χ2 probability and the dotted lines highlight the 1, 2 and 3σ confidence intervals. The vertical dashed lines shows the 1σ confidence interval of the expected α value derived from Eq. (2) of Balona & Abedigamba (2016).

In the text

Fig. 11 Mass-radius diagram. The grey point represent all the exo-planets from the Extrasolar Planets Encyclopaedia (http://www.exoplanet.eu/) as of 1 March 2022, with masses and radii known with a precision better than 30%. HD 3167 b and c are represented in the different instrument passbands analyzed in this work (CHEOPS. K2, HST/WFC3/G141 and Spizer/IRAC/Ch2). The black dash-dotted lines indicate two iso-density profiles matching the planets b and c. The colored solid and dashed lines are indicative chemical compositions as computed by Zeng et al. (2019).

In the text

	Fig. 12 Corner plot of the posterior distribution of the planetary parameters of HD 3167 e. The orbital period shows a very peculiar distribution that is not significantly correlated to any other parameter.
In the text

	Fig. 13 Posterior distributions of the planet-to-star radii ratios of planets b and c measured in the four instrument passbands.
In the text

	Fig. 14 5σ contrast curves of HD 3167 from Gemini North/’Alopeke at 562nm (blue) and 832 nm (red). The insert shows our reconstructed 832 nm speckle image.
In the text

Fig. 15 Corner plot showing the results on the interior composition models of HD 3167b (left) and HD 3167c, (right). The vertical dashed lines and the ‘error bars’ given at the top of each column represent the 5 and 95% percentiles. The internal structure parameters are the mass fraction of the core and the water layer with respect to the solid planet, the molar fraction of Si and Mg in the mantle, the molar fraction of Fe in the inner core and the logarithm of the gas mass in Earth masses.

In the text

	Fig. 16 Inclination evolution in the current system configuration, assuming the planets c, d and e are close to coplanar.
In the text

	Fig. 17 Evolution of the inclinations and eccentricities assuming the current star properties and a planet e initially misaligned with the system.
In the text

Fig. 18 Range of masses and semi-minor axis that allows a companion to tilt the outer planets of the system while preserving its mutual inclinations. The system is not destroyed for companions below the line βKL,ce = 1 and tilting is possible for a frequency ratio larger than 1. The companion has a significant influence on the system if its angular momentum is larger than the system’s angular momentum. Otherwise, the companion’s orbit precesses around the system angular momentum and does not induce large obliquity for the system. The green region is the 5σ limit set by direct imaging constraints. The purple and blue curve are the 5σ limit set by RVs, assuming a circular orbit and inclinations of 90 and 35°.

In the text

	Fig. 19 Posterior distributions of the initial atmospheric mass fractions for planets HD 3167 b and c derived by PASTA. The light-blue line represents the distribution of the estimated present-day atmospheric mass-fraction. The orange horizontal lines indicate the uninformative prior distributions.
In the text

Fig. 20 Posterior distribution (blue line) of the stellar rotation rate of HD 3167 after 150 Myrs derived by PASTA. The purple area represents the highest posterior density (HPD) interval of the distribution. The black line represents the distribution of the stellar rotation rate of young open cluster stars with mass comparable to that of HD 3167 based on the collection of data provided by Johnstone et al. (2015).

In the text

	Fig. 21 Mass-period (top panel) and mass-radius (bottom) diagrams of the exoplanet population in the Earth-Neptune range. HD 3167 planets are shown as green disks (HD 3167d and e are positioned at their measured minimum mass).
In the text

	Fig. A.1 CHEOPS observations of HD 3167 b. The raw data extracted by the DRP is shown in blue and the best-fit model is represented by the orange line. The transit model is shown in green with an offset for better visibility.
In the text

	Fig. B.1 HST observations of HD 3167 c. The left column shows the raw data in blue with the best model in orange for the 5 visits. The right column represents the detrended data with the transit model in orange. Visit 3 and 4 are affected by a serendipitous transit of HD 3167 b.
In the text