Nearly polar orbit of the sub-Neptune HD 3167 c Constraints on the dynamical history of a multi-planet system

Aims. We present the obliquity measurement, that is, the angle between the normal angle of the orbital plane and the stellar spin axis, of the sub-Neptune planet HD3167 c, which transits a bright nearby K0 star.We study the orbital architecture of this multi-planet system to understand its dynamical history. We also place constraints on the obliquity of planet d based on the geometry of the planetary system and the dynamical study of the system. Methods. New observations obtained with HARPS-N at the Telescopio Nazionale Galileo (TNG) were employed for our analysis. The sky-projected obliquity was measured using three different methods: the Rossiter-McLaughlin anomaly, Doppler tomography, and reloaded Rossiter-McLaughlin techniques. We performed the stability analysis of the system and investigated the dynamical interactions between the planets and the star. Results. HD3167 c is found to be nearly polar with sky-projected obliquity, λ = -97◦± 23◦. This misalignment of the orbit of planet c with the spin axis of the host star is detected with 97% confidence. The analysis of the dynamics of this system yields coplanar orbits of planets c and d. It also shows that it is unlikely that the currently observed system can generate this high obliquity for planets c and d by itself. However, the polar orbits of planets c and d could be explained by the presence of an outer companion in the system. Follow-up observations of the system are required to confirm such a long-period companion.


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 multiplanet 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 with the 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 +1.30 −1.24 M ⊕ for HD 3167 c, 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 skyprojected 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 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 −1 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.

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 signalto-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 −1 with a mean value of 0.9 m s −1 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). 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 χ 2 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.

Detection of a spectroscopic transit
To separate the observation taken during the planet transit, only RVs between the beginning of the ingress (T 1) and end of the egress (T 4) were considered. The photometric values of midtransit (T 0 ), period (P), and transit duration (T 14) of HD 3167 c along with their uncertainties were taken from Christiansen et al. (2017). The total uncertainty of ∼16 min on T 0 , inferred from the respective uncertainties of 15, 6, and 2 min on P, T 1/T 4, and T 0 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 0 , it was not clear whether the remaining 4 RVs were present inside or outside the transit. In the following analysis, T 0 is fixed to the photometric value as the uncertainty on T 0 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 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 χ 2 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 −1 , 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 −1 , 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 −1 , 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 −1 , 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) 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.

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.

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 0 for HD 3167 c 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 −1 ). We performed a grid search for the free parameters and computed χ 2 A28, page 3 of 12 A&A 631, A28 (2019) at each grid point. The contribution from the uncertainties of i p , a/R , and v sin i was also added quadratically to χ 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 χ 2 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 −1 , 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 χ 2 . The 1σ error bars were determined for all five free parameters following the χ 2 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 λ = −92 • +11 −20 , indicating a nearly polar orbit.
The derived v sin i (2.8 +1.9 −1.3 kms −1 ) 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 χ 2 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 χ 2 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 λ = −91 • +7 −16 , which is within the 1 σ uncertainty. The large v sin i obtained here (4.8 ± 2.1 kms −1 ) did not significantly affect the measurement of λ. Because the planetary 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 0 , T 14, 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.

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 0 , 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 −1 , which is small compared to the RM signal from HD 3167 c and the RV measurement accuracy. Lower panel: residuals after the best-fit RM is subtracted. 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 in Sect. 3.2.2. The residual average inside and outside the transit was found to be 0.23 ± 1.29 m s −1 , and 0.39 ± 1.66 m s −1 , 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 we considered the spectroscopic transit to be not significantly detected in the 2017 data and did not considered it for further analysis.

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 A28, page 4 of 12 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 . 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 0 , 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), 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 −1 and the spectra have a resolution of 7.5 km s −1 , 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, where σ Uncorr / √ n bin 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 −1 , which is consistent with the estimate from spectroscopy (v sin i = 1.7 ± 1.1 kms −1 ). 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 0 by varying it within 1 σ error bars. The value of λ remained within the 1 σ uncertainty derived above.

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 −1 and λ = −112.5 • +8.7 −8.5 , which agrees at 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 0 , 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.

Comparison between the three methods
The most commonly used method to estimate sky-projected obliquity using RV measurements is the analysis of the RM A28, page 5 of 12 A&A 631, A28 (2019)   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 the full 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 assumptions of 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.

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 0 the unit vector along the line of sight directed toward Earth and u 1 a unit vector perpendicular to u 0 , 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 i b = 83.4 • +4.6 −7.7 and i c = 89.3 • +0.5 −0.96 (Christiansen et al. 2017). For a planet k (here k stands for either b or d), we define φ k as the angle between u 1 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 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 nontransiting 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 were fulfilled, where R is the stellar radius and a d is the semimajor 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.

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.

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 6 realizations of the HD 3167 system by drawing from the best fit of masses, eccentricities, and semimajor 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 1 . 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 where G k = m k GM S a k (1 − e 2 k ) is the norm of the angular momentum of planet k. Then, we computed the total angular momentum deficit (AMD, Laskar 1997) of the system 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.  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.
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 dynamically excited 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 is about 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.

Interactions of planets and stellar spin
Because the system'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 threevector 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 2 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 2 The characteristic frequencies designate the coupling parameters between the different vectors, as explained in Appendix A. They have the dimension of a frequency, but are not properly speaking the frequencies of the system. Here, we use the terminology introduced in Boué & Fabrycky (2014). 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.

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 b = a √ 1 − e 2 , and its initial inclination I 0 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 0 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 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 0 . 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 0 for I 0 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  (8)). For β KL > 1, the outer companion can destabilize the observed system. measurement of the eccentricities of planets d and c will also help to constrain this scenario better.

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 −1 yr −1 before the Keplerian was fit. However, we detected a linear drift of about 7.6 ± 1.6 m s −1 yr −1 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 0 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.

Conclusion
We used new observations obtained with HARPS-N to measure the obliquity of a sub-Neptune in a multi-planetary system. The three different methods we applied on this challenging dataset A28, page 9 of 12 A&A 631, A28 (2019) 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.