© S. Dalal et al. 2019

1 Introduction

Obliquity is defined as the angle between the normal angle of a planetary orbit and the rotation axis of the planet host star. It is an important probe for understanding the dynamical history of exoplanetary systems. Solar system planets are nearly aligned and have obliquities lower than 7°, which might be a consequence of their formation from the protoplanetary disk. However, this is not the case for all exoplanetary systems. Various misaligned systems, that is, λ ≿ 30°, including some retrograde (λ ~ 180°, e.g., Hébrard et al. 2008) or nearly polar (λ ~ 90°, e.g., Triaud et al. 2010) orbits have been discovered. These misaligned orbits may result from Kozai migration and/or tidal friction (Nagasawa et al. 2008; Fabrycky & Tremaine 2007; Guillochon et al. 2011; Correia et al. 2011), where the close-in planets migrate as a result of scattering or of early-on interaction between the magnetic star and its disk (Lai et al. 2011), or the migration might be caused later by elliptical tidal instability (Cébron et al. 2011). Another possibility is that the star has been misaligned since the days when the protoplanetary disk was present as a result of inhomogeneous accretion (Bate et al. 2010) or a stellar flyby (Batygin 2012).

Most of the obliquity measurements are available for single hot Jupiters. Some of the smallest planets detected with a Rossiter measurement are GJ 436 b (4.2 ± 0.2 R_⊕) and HAT-P-11 b (4.4 ± 0.1 R_⊕), which are nearly polar (Bourrier et al. 2018; Winn et al. 2010), and 55 Cnc e (1.94 ± 0.04 R_⊕), which is also misaligned (Bourrier & Hébrard 2014), although the latest result has been questioned (López-Morales et al. 2014). Kepler 408 b is the smallest planet with a misaligned orbit among all planets that are known to have an obliquity measurement (Kamiaka et al. 2019). A few obliquity detections have been reported for multi-planet systems such as KOI-94 and Kepler 30 (Hirano et al. 2012; Albrecht et al. 2013; Sanchis-Ojeda et al. 2012), whose planets have coplanar orbits that are aligned withthe stellar rotation.

We study the multi-planet system hosted by HD 3167. This system includes two transiting planets and one non-transiting planet. Vanderburg et al. (2016) first reported the presence of two small short-period transiting planets from photometry. The third planet HD 3167 d was later discovered in the radial velocity (RV) analysis by Christiansen et al. (2017). Gandolfi et al. (2017) found evidence of two additional signals in the RV measurements of HD 3167 with periods of 6.0 and 10.7 days. However, they were unable to confirm the nature of these two signals. Furthermore, Christiansen et al. (2017) did not find any signal at 6 or 10.7 days. The masses of the transiting planets were found to be 5.02 ± 0.38 M_⊕ for HD 3167 b, a hot super-Earth, and 9.80 M_⊕ for HD 3167c, a warm sub-Neptune. The non-transiting planet HD 3167 d with a mass of at least 6.90 ± 0.71 M_⊕ orbits the star in 8.51 days. The two transiting planets have orbital periods of 0.96 days and 29.84 days and radii of 1.70 R_⊕ and 3.01 R_⊕, respectively.We measure the sky-projected obliquity for HD 3167 c, whose larger radius makes it the most favorable planet for the obliquity measurements. Because the period of planet c is longer than that of planet b, the data sampling during a given transit is three times better.

It is difficult to measure the true 3D obliquity, and most methods only access the projection of the obliquity. The sky-projected obliquity for a transiting exoplanet can be measured by monitoring the stellar spectrum during planetary transits. During a transit, the partial occultation of the rotating stellar disk causes asymmetric line profiles that can be detected using different methods such as the Rossiter-McLaughlin (RM) anomaly, Doppler tomography, and the reloaded RM method. These methods use different approaches to retrieve the path of the planet across the stellar disk. This allows us to quantify the systematic errors related to the data analysis method. The RM anomaly takes into account that asymmetry in line profiles induces an anomaly in the RV of the star (Queloz et al. 2000; Hébrard et al. 2008). However, changes in the cross-correlation function (CCF) morphology are not analyzed. Doppler tomography uses the spectral information present in the CCF of the star rather than just their RV centroids. This method entails tracking the full time-series of spectral CCF by modeling the additional absorption line profiles that are superimposed on the stellar spectrum during the planet transit (e.g., Collier Cameron et al. 2010; Bourrier et al. 2015; Crouzet et al. 2017). This model is then subtracted from the CCFs, and the spectral signature of the light blocked by the planet remains. Finally, the reloaded RM technique directly analyzes the local CCF that is occulted by the planet to measure the sky-projected obliquity (e.g., Cegla et al. 2016a; Bourrier et al. 2017). It isolates the CCFs outside and during the transit with no assumptions about the shape of the stellar line profiles.

The amplitude of the RM anomaly is expected to be below 2 m s⁻¹ for HD 3167 c. Detecting such a low-amplitude effect is challenging, therefore we decided to determine the robustness and significance of our results using the three different methods described above. The different methods have their respective advantages and limitations. A combined analysis involving the three complementary approaches therefore provides an obliquity measurement that is more robust against systematic effects that are due to the analysis method.

We measure the sky-projected obliquity of HD 3167 c using the three methods and finally discuss the dynamics of the system. This paper is structured as follows. We describe the spectroscopic observations during the transit in Sect. 2. The detection of spectroscopic transit followed by the data analysis using the RM anomaly, Doppler tomography, and the reloaded RM is presented in Sect. 3. We discuss the obliquity of planets b and d from geometry in Sect. 4. We study the dynamics of the system in Sect. 5 and explore the possible outer companion in Sect. 6. Finally, we conclude in Sect. 7.

2 Observations

We obtained the spectra of HD 3167 during the two transits of planet c on 2016 October 1 and 2017 November 23 with the spectrograph HARPS-N with a total of 35 observations and 24 observations, respectively. HARPS-N, which is located at the 3.58 m Telescopio Nazionale Galileo (TNG, La Palma, Spain), is an echelle spectrograph that allows high-precision RV measurements. Observations were taken with resolving power R = 115 000 with 15 min of exposure time. We used the spectrograph with one fiber on the star and the second fiber on a thorium-argon lamp so that the observation had high RV precision. The signal-to-noise ratio (S/N) per pixel at 527 nm for the spectra taken during the 2016 transit was 56–117 with an average S/N = 87. The 2017 transit was observed in poor weather conditions with S/N values ranging from 34 to 107 with an average S/N = 72. We primarily worked with the 2016 transit data for the reasons explained in Sect. 3.2.3.

The Data Reduction Software (DRS version 3.7) pipeline was used to extract the HARPS-N spectra and to cross-correlate them with numerical masks following the method described in Baranne et al. (1996) and Pepe et al. (2002). The CCFs obtained were fit by Gaussians to derive the RVs and their uncertainties. We tested different numerical masks such as G2, K0, and K5 and also determined the effect of removing some low S/N spectral orders to obtain the CCFs. These tests were performed to improve the data dispersion after the Keplerian fit. The method of fitting a Keplerian is discussed in detail in Sect. 3.1. Final RVs were obtained from CCFs that we produced using the K5 mask and removing the first 15 blue spectral orders with low S/N.

The resulting RVs with their uncertainties are listed in Table 1 for the 2016 observations and in Table B.1 for the 2017 observations. The typical uncertainties were between 0.6 and 1.5 m s⁻¹ with a mean value of 0.9 m s⁻¹ for the 2016 data. The stellar and planet parameters for HD 3167 that we used were taken from Tables 1 and 5 of Christiansen et al. (2017), except for the value of limb-darkening coefficient (ε), which was taken from Gandolfi et al. (2017).

3 Analysis

3.1 Detection of a spectroscopic transit

Figure 1 displays the RV measurements of HD 3167 during the 2016 transit of planet c. The upper panel shows RVs along with the best-fit RM model found from χ² minimization (discussed in Sect. 3.2), and the lower panel shows residual RVs after the fit. The red dashed line is the Keplerian model for the orbital motion of the three planets. During the transit, the deviation between the Keplerian model and the observed RVs is caused by the RM anomaly.

To separate the observation taken during the planet transit, only RVs between the beginning of the ingress (T1) and end of the egress (T4) were considered. The photometric values of mid-transit (T₀), period (P), and transit duration (T14) of HD 3167 c along with their uncertainties were taken from Christiansen et al. (2017). The total uncertainty of ~16 min on T₀, inferred from the respective uncertainties of 15, 6, and 2 min on P, T1/T4, and T₀ from Christiansen et al. (2017), was taken into account in determining the RVs outside the transit. Thirteen RVs (8 before and 5 after the transit) lay outside the transit, while 18 RVs were present inside the transit. Because of the uncertainty in the observed T₀, it was not clear whether the remaining 4 RVs were present inside or outside the transit. In the following analysis, T₀ is fixed to the photometric value as the uncertainty on T₀ is negligible in our analysis, as shown in Sect. 3.2.2.

The 13 RVs outside the transit were not sufficient for an independent Keplerian model for the three planets. We therefore took the orbital parameters for the three planets to fit the Keplerian from Table 5 of Christiansen et al. (2017) in Eq. (1) as (1)

Here K_i represents the RV semi-amplitude, the true anomaly and eccentricity are denoted by f_i and e_i, respectively, and ω_i is the argument of periastron. Finally, a Keplerian model was fit by minimizing the χ² considering only one free parameter, that is, the systemic velocity γ. The average of the residual RVs that were taken outside the transit was found to be 0.11 ± 0.72 m s⁻¹, in agreement with the expected uncertainties.

After the Keplerian fit, we noted that the average of residual RVs inside the transit was 1.17 ± 0.76 m s⁻¹, showing an indication of an RM anomaly detection. We fit this using the RM model in Sect. 3.2.2. According to Gaudi & Winn (2007), the expected amplitude of the RM anomaly is 1.7 m s⁻¹, which is within the order of magnitude of the deviation from the Keplerian model observed during the transit.

Furthermore, the slope that is visible in RVs within the observation time (8.7 h) was due to the short periodicity of HD 3167 b (P_b = 0.96 day). To compute the mass of HD 3167 b, a Keplerian in the RVs outside the transit was fit in which K_b was kept as a free parameter. K_b was found to be 3.86 ± 0.35 m s⁻¹, corresponding to a planet mass of HD 3167 b of M_b = 5.45 ± 0.50 M_⊕. This is consistent with the measurements of Christiansen et al. (2017) (K_b = 3.58 ± 0.26 m s⁻¹, M_b = 5.02 ± 0.38 M_⊕). K_b was fixed to the more accurate measurement of Christiansen et al. (2017) in the further analysis.

We note that the sky-projected obliquity λ was defined as the angle counted positive from the stellar spin axis toward the orbital plane normal, both projected in the plane of the sky. The sky-projected obliquity was fit using three different methods, as described in the following sections.

Fig. 1 RV measurements of HD 3167 taken on 2016 October 1 as function of time. Upper panel: solid black circles represent the HARPS-N data, the dashed red line indicates the Keplerian fit, and the solid green line depicts the final best fit with the RM effect. Lower panel: black solid circles are the residuals after subtracting the Keplerian, and green solid circles are the residuals after subtracting the best-fit RM model.

3.2 Rossiter-McLaughlin anomaly

The model to fit the RM anomaly is presented in the following section. We applied this model to fit both datasets to measure the sky-projected obliquity.

3.2.1 Model

The method developed by Ohta et al. (2005) was implemented to model the shape of the RM anomaly. These authors derived approximate analytic formulae for the anomaly in RV curves, considering the effect of stellar limb darkening. Following their approach, we adopted a model with five free parameters: γ, λ, the sky-projected stellar rotational velocity v sin i_⋆, the orbital inclination i_p, and the ratio of orbital semi-major axis to stellar radius a∕R_⋆. The values of the radius ratio r_p∕R_⋆, P, T₀ for HD 3167c were fixed to their photometric values (Christiansen et al. 2017), and ε for HD 3167 was fixed to 0.54 (Gandolfi et al. 2017). The parameters i_p and a∕R_⋆ were kept free because their values were poorly constrained from the photometry. Gaussian priors were applied to i_p and a∕R_⋆ as obtained from photometry (Christiansen et al. 2017). We adopted a value of v sin i_⋆ as a Gaussian prior based on the spectroscopy analysis in Christiansen et al. (2017) (v sin i_⋆ = 1.7 ± 1.1 kms⁻¹). We performed a grid search for the free parameters and computed χ² at each grid point. The contribution from the uncertainties of i_p, a∕R_⋆, and v sin i_⋆ was also added quadratically to χ².

3.2.2 2016 dataset

The data taken on 2016 October 1 are the best dataset for the obliquity measurement in terms of data quality and transit sampling. The 2016 data were fit with the Ohta model, and the reduced χ² with 30 degrees of freedom (n) for the best-fit model (RM fit) was found to be 0.95. With the RM fit, the averages of residuals inside and outside the transit were 0.01 ± 0.75 and 0.11 ± 0.72 m s⁻¹, respectively. The uncertainties agree with the expected uncertainties on the RVs (see Col. 3 of Table 1). The best-fit value for each parameter corresponds to a minimum of χ². The 1σ error bars were determined for all five free parameters following the χ² variation as described in Hébrard et al. (2002). The best-fit values together with 1 σ error bars are listed in Table 2. We measured λ = , indicating a nearly polar orbit.

The derived v sin i_⋆ (2.8 kms⁻¹) from the RM anomaly suggested a 2 σ detection of the spectral transit. In order to properly determine the significance of our RM detection, we performed Fischer’s classical test. The two models considered for the test were a K (only Keplerian) fit and an RM (Keplerian+RM) fit. The χ² for the K and RM fits is 63.55 (n = 34) and 28.76 (n = 30), respectively. A significant improvement was noted for the second model with F = 1.95 (p = 0.03) obtained using an F-test. The improvement to the χ² was attributed to the RM anomaly detection with 97% confidence. We conclude that the spectroscopic transit is significantly detected.

As a test, we applied a similar grid procedure without the spectroscopic constraint on v sin i_⋆ from Christiansen et al. (2017). We obtained λ = , which is within the 1 σ uncertainty. The large v sin i_⋆ obtained here (4.8 ± 2.1 kms⁻¹) did not significantly affect the measurement of λ. Because theplanetary orbit was found to be polar and it transits near the center of the star (b = 0.50 ± 0.32, Christiansen et al. 2017), the corresponding RM anomaly shape did not place a strong constraint on v sin i_⋆. The v sin i_⋆ can be estimated more accurately using the Doppler tomography technique in Sect. 3.3.

Furthermore, the effect of fixed parameters such as r_p∕R_⋆, P, T₀, T14, and K_b on λ was investigated. When these fixed parameters were varied within their 1 σ uncertainty, λ was found to remain within the 1 σ uncertainty derived above.

3.2.3 2017 dataset

Here, we evaluate whether the lower-quality 2017 dataset agrees with the results obtained above using the 2016 dataset. We first determine the observations taken outside the 2017 transit using the same method as explained in Sect. 3.1. After considering uncertainty on T₀, we found that only one RV measurement was taken clearly outside the transit. The scarcity of data and poor data sampling outside the transit and along with the low-quality observations during 2017 transit prevented us from finding a good model for a Keplerian and finally an independent value of λ. Thus the RM model parameters were fixed to the best-fit values from the 2016 transit, and the model derived previously was scaled to the RV level of this epoch. We also realized that during the 2017 transit, HD 3167 b and HD 3167 c transited simultaneously. However, the expected amplitude of the RM anomaly from HD 3167 b is 0.56 m s⁻¹, which is small compared to the RM signal from HD 3167 c and the RV measurement accuracy.

Figure 2 shows the best-fit RM model from Sect. 3.2.2 during the 2017 transit and the residuals after the best-fit RM was subtracted. This fit shows that the 2017 dataset roughly agrees with the results obtained from the RM anomaly fit for the 2016 observations; despite its lower quality, it did not invalidate the results presented inSect. 3.2.2. The residual average inside and outside the transit was found to be 0.23 ± 1.29 m s⁻¹, and 0.39 ± 1.66 m s⁻¹, respectively.The obtained uncertainties were slightly larger than the expected uncertainties on the RVs. The 2017 dataset presented short-term variations in the first half of the transit that could not be due to RM or Keplerian effects. We interpreted them as an artifact due to the bad weather conditions. We achieved no significant improvement from fitting the RM anomaly (F = 0.97, p = 0.44), therefore weconsidered the spectroscopic transit to be not significantly detected in the 2017 data and did not considered it for further analysis.

Fig. 2 RV measurement of HD 3167 taken on 2017 November 23 as a function of time. Upper panel: solid black circles represent the HARPS-N data, the dashed red line indicates the Keplerian fit, and the green line is the over plotted best-fit RM model from the 2016 transit. The blue dotted line marks the transit ingress and egress of planet b. The expected RM amplitude due to the transit of planet b is 0.6 m s−1. Lower panel: residuals after the best-fit RM is subtracted.

3.3 Doppler tomography

Here we present the obliquity measurement we performed on the 2016 dataset using Doppler tomography in order to compare it with the measurement from the RM anomaly technique presented above. When a planet transits its host star, it blocks different regions of the rotating stellar disk, which introduces a Gaussian bump in the spectral lines of the star. This bump can be tracked by inspecting the changes in the CCF, which allows us to measure the obliquity. The stellar rotational speed can also be obtained independent from the spectroscopic estimate by Christiansen et al. (2017). The CCFs obtained from the DRS with the K5 mask were used for this analysis (Sect. 2). Following the approach of Collier Cameron et al. (2010), we considered a model of the stellar CCF, which is the convolution of limb-darkened rotation profile with a Gaussian corresponding to the intrinsic photospheric line profile and instrumental broadening. When the CCFs are fit by the model including the stellar spectrum and the transit signature, some residual fixed patterns appear that are constant throughout the whole night. These patterns, also called “sidelobes” by Collier Cameron et al. (2010), are produced by coincidental random alignments between some stellar lines and the lines in the mask when the mask is shifted to calculate the CCFs. To remove these patterns, we assumed that they do not vary during the night, and we averaged the residuals of the out-of-transit CCFs after subtracting the best fit to the CCFs that was calculated by considering the stellar spectrum alone. We made a tomographic model that depended on the same parameters as the Keplerian plus RM model (Sect. 3.2), and added the local line profile width, s (non rotating local CCF width) expressed in units of the projected stellar rotational velocity (Collier Cameron et al. 2010). The most critical free parameters to fit the Gaussian bump were λ, v sin i_⋆, γ, i_p, a∕R_⋆, and s. Other parameters such as P, r_p ∕R_⋆, T₀, and ε were fixed to the same values as were used for the RM fit.

The following merit function was used to fit the CCFs following Bourrier et al. (2015), (2)

where f_i,j is the flux at the velocity point j in the i_th observed or model CCFs. The error on the CCF estimate was assumed to be constant over the full velocity range for a given CCF. To find the errors σ_i in the CCF profiles, we first used the constant errors, which are the dispersion of the residuals between the CCFs and the best-fit model profiles. As the CCFs were obtained using DRS pipeline with a velocity resolution of 0.25 km s⁻¹ and the spectra have a resolution of 7.5 km s⁻¹, the residuals were found to be strongly correlated. This led to an underestimation of the error bars on the derived parameters. A similar analysis as in Bourrier et al. (2015) was used to retrieve the uncorrelated Gaussian component of the CCFs. The residual variance as a function of data binning size (n_bin) is well represented by a quadratic harmonic combination of a white and red noise component, (3)

where is the intrinsic uncorrelated noise and σ_Corr is the constant term characterizing the correlation between the binned pixels. We adopted Gaussian priors for i_p and a∕R_⋆ from photometry (Christiansen et al. 2017).

The planet transit was clearly detected in the CCF profiles, as shown in Fig. 3. The v sin i_⋆ was found to be 2.1 ± 0.4 m s⁻¹, which is consistent with the estimate from spectroscopy (v sin i_⋆ = 1.7 ± 1.1 kms⁻¹). The sky-projected obliquity was measured to be λ = −88° ± 15°, which is in accordance with the result from the RM analysis (see Sect. 3.2.2). Table 2 lists the best-fit values together with 1 σ error bars.

We also performed a test to check the effect of the fixed parameter T₀ by varying it within 1 σ error bars. The value of λ remained within the 1 σ uncertainty derived above.

3.4 Reloaded Rossiter-McLaughlin technique

We applied the reloaded RM technique (Cegla et al. 2016a; Bourrier et al. 2018) to the HARPS-N observations of HD 3167 c. CCFs computed with the K5 mask (Sect. 2) were first corrected for the Keplerian motion of the star induced by the three planets in the system (calculated with the properties from Christiansen et al. 2017). The CCFs outside of the transit were co-added to build a master-out CCF, whose continuum was normalized to unity. The centroid of the master-out CCF, derived with a Gaussian fit, was used to align the CCFs in the stellar rest frame. The continuum of all CCFs was then scaled to reflect the planetary disk absorption by HD 3167 c, using a light curve computed with the batman package (Kreidberg 2015) and the properties from Christiansen et al. (2017). Residual CCFs were obtained by subtracting the scaled CCF from the master-out (Fig. 4).

No spurious features are observed in the residual CCFs out of the transit. Within the transit, the residual RM spectrally and spatially resolve the photosphere of the star along the transit chord. The average stellar lines from the planet-occulted regions are clearly detected and were fit with independent Gaussian profiles to derive the local RVs of the stellar surface. We used a Levenberg-Marquardt least-squares minimization, setting flux errors on the residual CCFs to the standard deviation in their continuum flux. Because the CCFs are oversampled in RV, we kept one in four points to perform the fit. All average local stellar lines were well fit with Gaussian profiles, and their contrast was detected at more than 3 σ (using the criterion defined by Allart et al. 2017). The local RV series was fit with the model described in Cegla et al. (2016a) and Bourrier et al. (2017), assuming solid-body rotation for the star. We sampled the posterior distributions of v sin i_⋆ and λ using the Markov chain Monte Carlo (MCMC) software emcee (Foreman-Mackey et al. 2013), assuming uniform priors. Best-fit values were set to the medians of the distributions, with 1 σ uncertainties derived by taking limits at 34.15% on either side of the median. The best-fit model shown in Fig. 4 corresponds to v sin i_⋆ = 1.9 ± 0.3 km s⁻¹ and λ = −112.5, which agreesat better than 1.4 σ with the results obtained from the RM and Doppler tomography (Sects. 3.2 and 3.3). The error bars on λ are small because i_p and a∕R_⋆ were fixed in this particular analysis. However, when i_p, T₀, and a∕R_⋆ were varied within their 1 σ uncertainty, λ did not vary significantly and remained within 1 σ uncertainty. The best-fit values with their 1 σ uncertainties are listed in Table 2.

Fig. 3 Maps of the time-series CCFs as a function of RV relative to the star (in abscissa) and orbital phase (in ordinate). The dashed vertical white lines are marked at ± v sin i⋆, and first and fourth contact of transit is indicated by white diamonds. Upper panel: map of the transit residuals after the model stellar profile was subtracted. The signature of HD 3167 c is the moderately bright feature that is visible from ingress to egress. Middle panel: transiting signature of HD 3167 c using the best-fit model, obtained with λ = − 88°. Lower panel: overall residual map after the model planet signature was subtracted.

3.5 Comparison between the three methods

The most commonly used method to estimate sky-projected obliquity using RV measurements is the analysis of the RM anomaly. However, the RM method does not exploit the full spectral CCF. In some extreme cases, the classical RM method can introduce large biases in the sky-projected obliquity because of asymmetries in the local stellar line profile or variations in its shape across the transit chord (Cegla et al. 2016b). The Doppler tomography method is less affected than the RM anomaly method because it explores thefull information in the CCF. However, a bias in the obliquity measurements can also be introduced by assuming a constant, symmetric line profile and ignoring the effects of the differential rotation. Results from the reloaded RM technique suggest that the bias is not significant here. The reloaded RM technique does not make prior assumptionsof the local stellar line profiles and allows a clean and direct extraction of the stellar surface RVs along the transit chord. This results in an improved precision on the obliquity, albeit under the assumption that the transit light-curve parameters (in particular the impact parameter and the ratio of the planet-to-star radius) are known to a good enough precision to be fixed. In the present case, we might thus be underestimating the uncertainties on λ with this method.

The sky-projected obliquities measured by all three methods agree to better than 1.4 σ. This confirms that the spectroscopic transit in the 2016 data is significantly detected and suggests that the corresponding obliquity measurement is not reached by strong systematics that would be due to the method. Combining the λ values from all three methods, we estimated the sky-projected obliquity for HD 3167 c to be λ = −97° ± 23°, after taking into account both the systematic and statistical errors. We adopted this conservative value in our final obliquity measurement.

As discussed in Sect. 3.2.2, the stellar rotation speed was poorly constrained by the RM method. However, the v sin i_⋆ more accurately measured from Doppler tomography and the reloaded RM technique was consistent with the measurements of Christiansen et al. (2017). The v sin i_⋆ from three methods was also found to be within 1 σ. Furthermore, the two photometric parameters a∕R_⋆ and i_p also agreed within their uncertainties for the RM and Doppler tomography methods. The systemic velocity γ is slightly different in each case because a different definition was employed in each method.

Fig. 4 Upper panel: map of the residual CCF series as a function of orbital phase (in abscissa) and RV in the stellar rest frame (in ordinate). Colors indicate flux values. The four vertical dashed black lines show the times of transit contacts. The in-transit residual CCFs correspond to the average stellar line profiles from the regions that are occulted by HD 3167 c across the stellar disk. The solid black line is the best-fit model to the local RVs of the planet-occulted regions (λ = −112.5°), assuming solid-body rotation for the star (v sin i⋆ = 1.89 km s−1). Lower panel: RVs of the stellar surface regions occulted by the planet (blue points), best fit with thesolid black line (same as in the upper panel). The gray area corresponds to the 1 σ envelope of the best fit, derived from the MCMC posterior distributions. The dashed red line shows a model obtained with the same stellar rotational velocity, but an aligned orbit (λ = 0°). This highlights the large orbital misalignment of HD 3167 c.

4 Obliquity of planets b and d from geometry

The spectroscopic transit observations gave constraints only on the obliquity of planet c. Although planet b is also transiting, the low amplitude for the RM signal during the transit precludes measuring its obliquity with the present data. However, because both planets b and c are transiting planets, the mutual inclination can be constrained.

We denote by u₀ the unit vector along the line of sight directed toward Earth and u₁ a unit vector perpendicular to u₀, that is, in the plane of the sky (see Fig. 5). The orbital planes of planets b and c are characterized by the perpendicular unit vectors u_b and u_c. The inclination of their orbits, i_b and i_c, is constrained to be and (Christiansen et al. 2017). For a planet k (here k stands for either b or d), we define ϕ_k as the angle between u₁ and the projection of u_k on the plane of the sky (this is equivalent to the longitude of the ascending node in the plane of the sky). With these definitions, the mutual inclination between the planets b and c, i_bc, is given by (4)

With cos i_b and cos i_c uniformly distributed within their 1 σ error bars and assuming that ϕ_b and ϕ_c are uniformly distributed between 0 and 2π, we calculated the probability distribution of i_bc (Fig. 6). The probability distribution was found to be close to a uniform distribution, except that it is low for i_bc below 10° and above 170°. Based on geometry, no information on the obliquity of planet b can therefore be derived from our measurement of the obliquity of planet c. We note that in the case of two non-transiting planets, the probability distribution of i_bc would peak around 90°, as shown by the dotted line in Fig. 6.

Planet d would transit if the condition (5)

were fulfilled, where R_⋆ is the stellar radius and a_d is the semi-major axis of planet d. Because planet d does not transit, the mutual inclination between planets c and d must be at least 2.3°.

As a result, the obliquity of planets b and d cannot be constrained well from the geometry of the planetary system alone. We place constraints on the obliquity of planet d from the dynamics of the planetary system in Sect. 5.

Fig. 5 Pictorial representation of the reference angles and the unit vectors. uS corresponds to the direction of the stellar spin.

Fig. 6 Probability distribution of the mutual inclination between the planets b and c (solid line). For comparison, the dotted line shows the probability distribution when neither planet transits.

5 Dynamics

We study the dynamics of the system to investigate the interactions between planets and stellar spin which could explain the polar orbit of planet c. We also perform the Hill stability analysis to set bounds on the obliquity of planet d in the following section.

5.1 Planet mutual inclinations

While the available observations were unable to geometrically constrain the mutual inclination of the planets, a bound is given by the stability analysis of the system. Short-period planets with an aligned orbit such as KELT-24 b and WASP-152 b (Rodriguez et al. 2019; Santerne et al. 2016), or with an misaligned orbit such as Kepler-408 b and GJ436 b (Kamiaka et al. 2019; Bourrier et al. 2018) have been detected. The obliquity distribution of short-period planets is not clear. However, because planet b is close to the star, its orbit is most likely circular and its inclination is governed by the interaction with the star, as shown in Appendix A.2. The exact inclination of planet b is not important from a dynamical point of view, and it is safe to neglect the influence of planet b when the stability of the system is investigated. We focus here on the outer pair of planets to constrain the system and study the simplified system that is only composed of the star and the two outer planets. Our goal is to determine the maximum mutual inclination between planets d and c such that the outer pair remains Hill-stable (Petit et al. 2018; Marchal & Bozis 1982). We first created 10⁶ realizations of the HD 3167 system by drawing from the best fit of masses, eccentricities, and semi-major axis distributions given by Christiansen et al. (2017). To each of these copies of the system, we set the mutual inclination between planets c and d with uniformly spaced values of between 0° and 90°.

We assumed that the orbit of planet b is in the invariant plane, that is, the plane perpendicular to the angular momentum vector of the whole system¹. As a result, we computed the inclinations i_c and i_d with respect to the invariant plane because the projection of the angular momentum onto the invariant plane gives (6)

where is the norm of the angular momentum of planet k. Then, we computed the total angular momentum deficit (AMD, Laskar 1997) of the system (7)

and we determined whether the pair d–c is Hill-stable. To do so, we compared the AMD to the Hill-critical AMD of the pair (Eq. (30), Petit et al. 2018). We plot in Fig. 7 the proportion of the Hill-stable system binned by mutual inclination i_dc. We also plot the proportion of the Hill-stable pair for a system with circular orbit.

We observe that for an inclination i_dc below 21°, the system is almost certainly Hill-stable. This means that for any orbital configuration and masses that are compatible with the observational constraints, the system will be long-lived with this low mutual inclination. We emphasize that long-lived configurations with higher mutual inclination than 21° exist. Christiansen et al. (2017) gave the example of Kozai-Lidov oscillations with initial mutual inclinations of up to 65°. However, the choice of initial conditions is fine-tuned because of the circular orbits (a configuration that is rather unlikely for such dynamicallyexcited systems).

When we assume that the stellar spin is aligned with the total angular momentum of the planets, the planet obliquity corresponds to the planet inclination with respect to the invariant plane. When we assume i_dc < 21°, the maximum obliquity of planet c isabout 9°. Even if the mutual inclination i_dc = 65°, the obliquity only reaches 32°. Thus, the observed polar orbit shows that the stellar spin cannot be aligned with the angular momentum of the planet.

From Sect. 4 and the previous paragraphs, we deduce that the most likely value for i_dc is between 2.3° and 21°. Because the mutual inclination of planets c and d is low, we can conclude that planet d is also nearly polar.

Fig. 7 Probability of the pair d–c to be Hill-stable as a function of the mutual inclination of d and c, assuming planet b is within the invariant plane. The masses, semi-major axes, and eccentricity are drawn from the best-fit distribution (Christiansen et al. 2017). The dashed curve corresponds to a system where every planet is on a circular orbit.

5.2 Interactions of planets and stellar spin

Because thesystem’s eccentricities and mutual inclinations are most likely low to moderate, we considered the interaction between the stellar spin and the planetary system. In particular, we investigated whether the motion of the planets can effectively tilt the star up to an inclination that could explain the polar orbit of planets c and d. The currently known estimate of Christiansen et al. (2017) of the stellar rotation period is 27.2 ± 7 days, but the period may have slowed down by a factor 10 (Bouvier 2013). In order to investigate the evolution of the obliquity that could have occurred in earlier stages in the life of the system, we studied the planet-star interaction as a function of the stellar rotation period.

To do so, we applied the framework of the integrable three-vector problem to the star and the angular momenta of planets d and c (Boué & Laskar 2006; Boué & Fabrycky 2014; Correia 2015). This model gives both the qualitative and quantitative behavior of the evolution of three vectors that represent different angular momenta directions u_S, u_d, and u_c under their mutual interactions. We describe the model in Appendix A.

As shown in Boué & Fabrycky (2014), the mutual interactions of the three vectors can be described by comparing the different characteristic frequencies² of the system ν^d∕S, ν^S∕d, ν^d∕c, and ν^c∕d with the expressions given in Eq. (A.5). The frequency ν^j∕k represents the relative influence of the body j over the motion of u_k. In other words, if ν^k∕j ≪ ν^j∕k, u_j is almost constant while u_k precesses around. We here neglect the interactions between the star and planet c versus the interaction between the star and planet d because they are smaller by two orders of magnitude.

Because it is coupled with the star, planet b acts as a bulge on the star that enhances the coupling between the orbits of the outer planets and the star (see Appendix A.2). We limit our study to the configuration where the strongest coupling occurs, that is, when the orbit of planet b lies within the stellar equatorial plane. The influence of planet b modifies the characteristic frequencies ν^d∕S and ν^S∕d, as we show in Appendix A.2. The model is valid in the secular approximation if the eccentricities of planets d and c remain low such that G_d and G_c are constant. Boué & Laskar (2006) showed that the motion is quasi-periodic. It is possible to give the maximum spin-orbit angle of planet c as a function of the initial inclination of planet d.

Using the classification of Boué & Fabrycky (2014), we can determine the maximum misalignment between u_S and u_c as a function of the initial inclination between u_S and u_d. We plot the frequencies (cf. Eq. (A.5)) as a function of the stellar rotation period in Fig. 8. We merged the curves that represent ν^d∕c and ν^c∕d into ν^dc because the two terms are almost equal.

We are in a regime where (ν^d∕c ~ ν^c∕d) ≫ (ν^d∕S, ν^S∕d) and the orbital frequencies dominate the interactions with the star. For the shorter periods we have (ν^d∕c ~ ν^c∕d ~ ν^S∕d) ≫ ν^d∕S. Nonetheless, in both cases the dynamics are purely orbital, however, meaning that the star acts as a point mass and is never coupled with the orbits of the outer planets. It is not possible for planets c and d to reach a high mutual inclination with the stellar spin axis starting from almost coplanar orbits or an even moderate inclination. When planet b is misaligned with the star, it is even harder for the planets to tilt the star.

We conclude that even if the star has had a shorter period in the past, it is unlikely that the currently observed system can by itself generate such a high obliquity for planets c and d. However, high initial obliquities are almost conserved, which means that the observed polar orbits are possible under the assumptions made, even though they are not explained by this scenario.

Fig. 8 Characteristic frequencies defined in Eq. (A.5) as a function of stellar period. The current estimated stellar rotation period is marked with a vertical dashed line. The two terms νd∕c and νc∕d are merged into a single curve νdc because they are almost equal.

5.3 System tilt due to an unseen companion

We now assume that while the system only presents moderate inclinations, a distant companion on an inclined orbit exists. We consider the configurations that can cause the system to be tilted with respect to the star. Once again, we used the framework of Boué & Fabrycky (2014). We considered the vectors u_S, u, and u′ that give the direction of the stellar spin axis S, the total angular momentum of the planetary system G, and the angular momentum of the companion G′, respectively. The outer companion is described by its mass m′, its semi-major axis a′, its semi-minor axis , and its initial inclination I₀ with respect to the rest of the system, which is assumed to be nearly coplanar or to have moderate inclinations. Moreover, we assumed that G is initially aligned with S, while the companion is highly inclined with respect to the planetary system, that is to say, I₀ is larger than 45° up to 90°. According to Boué & Fabrycky (2014), all interactions between planets cancel out because we only consider the dynamics of their total angular momentum G.

As in the previous part, we can compare the different characteristics frequencies of the system ν^pla∕S, ν^S∕pla, ν^comp∕pla, and ν^pla∕comp of expression given in Eqs. (A.6) and (A.7). The companion effectively tilts the planetary system as a single body if its influence on planet c is weaker than the interaction between planets d and c. In the other case, planet c will enter Lidov-Kozai oscillations, which can lead to the destabilization of the system through the interactions with planet d. Boué & Fabrycky (2014) reported the limit at which the outer companion starts to perturb the planetary system and excites the outer planet through Kozai-Lidov cycles. They explained that if the coefficient β_KL is defined as (8)

and verifies β_KL ≪ 1, the companion’s influence does not perturb the system and tilts it as a whole.

We plot in Fig. 9 the frequencies ν^pla∕S,  ν^S∕pla,  ν^comp∕pla, and ν^pla∕comp as a function of β_KL and observe different regimes. In the first regime, we have β_KL < 0.1 and ν^pla∕comp ≪ ν^comp∕pla ≪ ν^S∕pla. The influence of the companion is too weak to change the obliquity of the planetary system. For β_KL > 1, we have (ν^pla∕comp ~ ν^comp∕pla) ≫ ν^S∕pla, in which regime the system obliquity can reach I₀. However, the companion destabilizes the orbit of planet c, which can lead to an increase in eccentricity and mutual inclination between the planets. For 0.1 ≲ β_KL ≲ 1, we remark that ν^pla∕S ≪ (ν^pla∕comp ~ ν^comp∕pla ~ ν^S∕pla). According to Boué & Fabrycky’s classification, the maximum possible inclination between the star and the planet, that is, their obliquity, is almost twice I₀ for I₀ ≲ 80°. In this regime, an unseen companion can explain the observed polar orbits of HD 3167 c and d.

We conclude that some stable configurations with an additional outer companion may explain the high obliquity of planets c and d. We further discuss the possible presence of outer companion signals in the existing RV data in Sect. 6. Accurate measurement of the eccentricities of planets d and c will also help to constrain this scenario better.

Fig. 9 Characteristic frequencies defined in Eq. (A.6) as a function of βKL (see Eq. (8)). For βKL > 1, the outer companion can destabilize the observed system.

6 Outer companion

To find the possible signatures of an outer companion, we performed two different tests on the RV data from Christiansen et al. (2017), which cover a span of five months. First we obtained the residual RV after we removed the Keplerian signal caused by the three planets. In the analysis performed by Christiansen et al. (2017), the linear drift was fixed to 0 m s⁻¹yr⁻¹ before the Keplerian was fit. However, we detected a linear drift of about 7.6 ± 1.6 m s⁻¹yr⁻¹ in the residual velocities. When we assume a circular orbit for the outer companion, this linear drift corresponds to a period of at least 350 days and a mass of at least 0.1 M_Jup. A body like this has a β_KL ≃ 0.08, which makes it unlikely that it is able to incline planets c and d with respect to the star.

Second, we generated the periodogram of the RV before and after we removed the known periodic signals of the three planets using the Lomb-Scargle method, as shown in Fig. 10. In addition to the detected planets, two other peaks at 11 days and 78 days were found at a false-alarm probability (FAP) higher than 0.1% in the Fourier power. The peak at 11 days was an alias caused by the concentration of the sampling around lunar cycles, as explained in Christiansen et al. (2017). In the lower panel of Fig. 10, no peak at 11 days was detected, but the peak at 78 days was persistent in both periodograms. The peak around 20 days in the lower panel may be caused by stellar rotation, and the peak around one day was an alias due to data sampling. When we assume a circular orbit with a period of 78 days, this corresponds to a mass of at least 0.03 M_Jup for the outer companion, which gives β_KL ≃ 0.5. This potential outer companion might explain the high obliquity of HD 3167 c if its initial inclination I₀ was high enough.

We found possible indications of an additional outer companion in the system. Additional RV observations of HD 3167 on a long time span are necessary to conclusively establish its presence and determine its orbital characteristics, and thus confirm (or refute) our hypothesis.

7 Conclusion

We used new observations obtained with HARPS-N to measure the obliquity of a sub-Neptune in a multi-planetary system. Thethree different methods we applied on this challenging dataset agree, which means that the sky-projected obliquity we measured is reliable. We report a nearly polar orbit for the HD 3167 c with λ ~ − 97° ± 23°. The measurements of λ from RM anomaly, Doppler tomography, and reloaded RM technique agree at better than 1.4 σ standard deviation with this value. The v sin i_⋆ from the three methods also agree within their uncertainties. To our knowledge, we are the first to apply these three methods and compared them to the spectroscopic observation of a planetary transit.

These observations are a valuable addition to the known planetary obliquity sample, extending it further beyond hot Jupiters. Several small-radius multi-planet systems with aligned spin-orbits such as Kepler 30 (Sanchis-Ojeda et al. 2012)and with a misaligned spin-orbit such as Kepler 56 (Huber et al. 2013) have been reported. Additionally, single small exoplanets with high-obliquity measurement such as Kepler 408 (Kamiaka et al. 2019) and GJ436 (Bourrier et al. 2018) have also been reported. Some of the misalignments might be explained by the presence of an outer companion in the system. One particularly interesting planetary system is Kepler 56, in which two of its transiting planets are misaligned with respect to the rotation axis of their host star. This misalignment was explained by the presence of a massive non-transiting companion in the system (Huber et al. 2013). A third planet in the Kepler 56 system was later discovered by Otor et al. (2016). This supported the finding of Huber et al. (2013). Similarly, the misalignment in HAT-P-11 b may be explained by the presence of HAT-P-11 c (Yee et al. 2018).

Our dynamical analysis of the system HD 3167 places constraints on the obliquity of planet d. We cannot determine the obliquity of planet b with the current data and information about the system. The Hill-stability criterion shows that the orbits of planets c and d are nearly coplanar, so that both planets are in nearly polar orbits. The interactions of the planets with the stellar spin cannot satisfactorily explain the polar orbits of planets c and d. We postulate that an additional unseen companion exists in the system. This might explain the polar orbits of planets c and d. Indications for additional outer companions are present in the available RV dataset. Continued RV measurements of HD 3167 on a longer time span might reveal the outer companion and confirm our speculation.

Fig. 10 Upper panel: Lomb-Scragle periodogram of the RV of HD 3167 taken from Christiansen et al. (2017). The black dashed line represents the false-alarm probability at 0.1%, and the three vertical red dashed lines correspond to the periods of the three planets that are currently detected around HD 3167. Lower panel: Lomb-Scragle periodogram of the RV data after the three known periodic signals are removed.

Acknowledgements

A.P. thanks G. Boué for the helpful discussions during the study. A.C. acknowledges support from the CFisUC strategic project (UID/FIS/04564/2019), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal. V.B. acknowledges support by the Swiss National Science Foundation (SNSF) in the frame of the National Centre for Competence in Research PlanetS, and has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Four Aces; grant agreement No 724 427). S.D. acknowledges Vaibhav Pant for insightful discussions.

Appendix A Details on the three-vector model

A.1 Generic three-vector problem

The three-vector problem (Boué & Laskar 2006; Boué & Fabrycky 2014) studies the evolution of the direction of three angular momenta that are represented by the unit vectors u_k for k = 1, 2, 3 in equations (A.1)

The constants ν^k∕j are called the characteristic frequencies and represent the relative influence of the body k over the evolution of j. Their expression depends on the considered problem. The three-vector problem is integrable (Boué & Laskar 2006), and the solution is quasi-periodic with two different frequencies. Given an initial state where two vectors are aligned and a third is misaligned, it is possible to compute the maximum inclination between the two initially aligned vectors as a function of the initial inclination with the third (Boué & Fabrycky 2014). The maximum inclination depends on the characteristic frequencies, and the different cases have been classified in Sect. 5.3 of Boué & Fabrycky (2014).

A.2 Influenceof planet b

In Sects. 5.1 and 5.2, we claimed that the inclination dynamics of planet b are most likely governed by the star and only influence planets d and c through a modification of the planet–star coupling. We present here the justification for this assumption as well as details on the expressions of the coupling constants.

We first focus on the three-vector problem (u_S, u_b, and u_d). For now, we neglect the effect of planet c because we focus on the dynamics of planet b. Following Boué & Laskar (2006) and Boué & Fabrycky (2014), the characteristic frequencies that appear in Eq. (A.1) are expressed as (A.2)

where G_k is the angular momentum of planet k, S = C_Sω_S is the angular momentum of the stellar rotation, with C_S the stellar moment of inertia and (A.3)

Here α_Sj represents the coupling between the star and planet j, and β_jk is the Laplace-Lagrange coupling between planets j and k (we assume a_j < a_k). We also define (A.4)

the gravitational quadrupole coefficient (Lambeck 1988), where k₂ is the second fluid Love number of the star and ω_S is the stellar rotation speed. For the numerical values of k₂ and C_S, we use Landin et al. (2009). For a star of mass 0.85 M_⊙, we have k₂ = 0.018 and .

Fig. A.1 Characteristic frequencies defined in Eq. (A.2) as a function of the stellar period. The current stellar rotation period is marked with a vertical dashed line. νb∕S dominates for most of the considered frequencies.

Independently of the stellar rotation speed, we have α_Sd ∕α_Sb ≤ 0.04. We therefore neglect the terms depending on α_Sd in this analysis. As a result, we can directly apply the results of the analysis reported by Boué & Fabrycky (2014), with the four characteristic frequencies ν^b∕S,  ν^S∕b,  ν^b∕d, and ν^d∕b.

We plot the frequencies ν^k∕j as a function of the stellar period in Fig. A.1. We average the frequencies in each point by randomly drawing the orbital elements from the best fit.

For the considered range of the stellar revolution period, ν^S∕b dominates all other frequencies, and it becomes comparable to ν^d∕b for the current rotation rate. Using the regime classification of Boué & Fabrycky (2014), we can determine the maximum misalignment between u_S and u_b as a function of the initial inclination between u_S and u_d.

For a faster-rotating star (i.e., a younger star), we have ν^S∕b ≫ (ν^b∕S, ν^d∕b, ν^b∕d), in which regime no significant misalignment of planet b can be achieved. As a result, planet b is completely coupled with the star and remains within its equator even if the other planets are mutually inclined. The current rotation rate leads to the so-called Laplace regime where (ν^S∕b ~ ν^d∕b) ≫ (ν^b∕S, ν^b∕d), in which the plane of planet b oscillates between the stellar equatorial plane and the plane of planet d. However, if the initial mutual inclination between planet b and the stellar equator is low, planet b remains close to the stellar equator.

We simplify the problem by considering that planet b is coupled to the star and modifies the stellar precession coupling constant α_{S k} for planets d and c. The modification of the coupling constant can be found in Boué & Laskar (2006, Eq. (129))³. While the expression was derived for a planet in the presence of a satellite, it remains valid in our case. The expression is a generalization of the approximations for close satellites (Tremaine 1991) and far satellites (d’Alembert 1749). We denote with the modified coupling constant to include the effect of planet b when it is considered as a bulge on the star.

A.3 Coupling constants for the interactions of the planet with the star

In Sect. 5.2 we considered the three vectors u_S, u_d, and u_c. The coupling frequencies that appear in Eq. (A.1) for this particular problem are given by (A.5)

where β_bd is defined in Eq. (A.3) and is the coupling between the star and planet k, modified to take the influence of planet b into account, as explained in Appendix A.2. We also simplified Eqs. (A.1) by neglecting over because independently of the stellar rotation period.

A.4 Coupling constants for the problem with a companion

The characteristic frequencies that govern the evolution of the inclination of the planets under the influence of an outer companion as explained in Sect. 5.3 are given by (A.6)

where (A.7)

We neglect the interaction between the star and the companion and as a result disregard the corresponding characteristic frequencies.

Appendix B Radial velocity data

References

All Tables

All Figures

Fig. 1 RV measurements of HD 3167 taken on 2016 October 1 as function of time. Upper panel: solid black circles represent the HARPS-N data, the dashed red line indicates the Keplerian fit, and the solid green line depicts the final best fit with the RM effect. Lower panel: black solid circles are the residuals after subtracting the Keplerian, and green solid circles are the residuals after subtracting the best-fit RM model.

In the text

Fig. 2 RV measurement of HD 3167 taken on 2017 November 23 as a function of time. Upper panel: solid black circles represent the HARPS-N data, the dashed red line indicates the Keplerian fit, and the green line is the over plotted best-fit RM model from the 2016 transit. The blue dotted line marks the transit ingress and egress of planet b. The expected RM amplitude due to the transit of planet b is 0.6 m s−1. Lower panel: residuals after the best-fit RM is subtracted.

In the text

Fig. 3 Maps of the time-series CCFs as a function of RV relative to the star (in abscissa) and orbital phase (in ordinate). The dashed vertical white lines are marked at ± v sin i⋆, and first and fourth contact of transit is indicated by white diamonds. Upper panel: map of the transit residuals after the model stellar profile was subtracted. The signature of HD 3167 c is the moderately bright feature that is visible from ingress to egress. Middle panel: transiting signature of HD 3167 c using the best-fit model, obtained with λ = − 88°. Lower panel: overall residual map after the model planet signature was subtracted.

In the text

Fig. 4 Upper panel: map of the residual CCF series as a function of orbital phase (in abscissa) and RV in the stellar rest frame (in ordinate). Colors indicate flux values. The four vertical dashed black lines show the times of transit contacts. The in-transit residual CCFs correspond to the average stellar line profiles from the regions that are occulted by HD 3167 c across the stellar disk. The solid black line is the best-fit model to the local RVs of the planet-occulted regions (λ = −112.5°), assuming solid-body rotation for the star (v sin i⋆ = 1.89 km s−1). Lower panel: RVs of the stellar surface regions occulted by the planet (blue points), best fit with thesolid black line (same as in the upper panel). The gray area corresponds to the 1 σ envelope of the best fit, derived from the MCMC posterior distributions. The dashed red line shows a model obtained with the same stellar rotational velocity, but an aligned orbit (λ = 0°). This highlights the large orbital misalignment of HD 3167 c.

In the text

	Fig. 5 Pictorial representation of the reference angles and the unit vectors. uS corresponds to the direction of the stellar spin.
In the text

	Fig. 6 Probability distribution of the mutual inclination between the planets b and c (solid line). For comparison, the dotted line shows the probability distribution when neither planet transits.
In the text

	Fig. 7 Probability of the pair d–c to be Hill-stable as a function of the mutual inclination of d and c, assuming planet b is within the invariant plane. The masses, semi-major axes, and eccentricity are drawn from the best-fit distribution (Christiansen et al. 2017). The dashed curve corresponds to a system where every planet is on a circular orbit.
In the text

	Fig. 8 Characteristic frequencies defined in Eq. (A.5) as a function of stellar period. The current estimated stellar rotation period is marked with a vertical dashed line. The two terms νd∕c and νc∕d are merged into a single curve νdc because they are almost equal.
In the text

	Fig. 9 Characteristic frequencies defined in Eq. (A.6) as a function of βKL (see Eq. (8)). For βKL > 1, the outer companion can destabilize the observed system.
In the text

Fig. 10 Upper panel: Lomb-Scragle periodogram of the RV of HD 3167 taken from Christiansen et al. (2017). The black dashed line represents the false-alarm probability at 0.1%, and the three vertical red dashed lines correspond to the periods of the three planets that are currently detected around HD 3167. Lower panel: Lomb-Scragle periodogram of the RV data after the three known periodic signals are removed.

In the text

	Fig. A.1 Characteristic frequencies defined in Eq. (A.2) as a function of the stellar period. The current stellar rotation period is marked with a vertical dashed line. νb∕S dominates for most of the considered frequencies.
In the text