Obliquity Constraints on the Planetary-mass Companion HD 106906 b

Observations of HD 106906 b with the Immersion Grating Infrared Spectrometer (IGRINS; Yuk et al. 2010; Park et al. 2014) on the Gemini South telescope (Mace et al. 2018) were completed 2020 February 04, 07, 08, 09 UT as part of program GS-2020A-Q-135 (PI: Bryan). These observations simultaneously covered H and K bands from ∼1.45–2.52 μm. We observed the K = 15.5 mag companion with individual image exposure times of 1528 s. All observations were taken with a slit orientation perpendicular to the 307.3° position angle (PA) between the host star and the companion (angular separation 71) in order to prevent a flux gradient across the slit. On 2020 February 04 UT we acquired three pairs of AB nodded exposures, amounting to three epochs of observation in ∼2.5 hr of on-source integration time. On the nights of 2020 February 07, 08, and 09 UT ABBA-nodded exposures were combined, producing three additional epochs. In total, there were six epochs of observation from four nights in early 2020 February.

We reduce all data with the IGRINS Pipeline Package (PLP; Lee & Gullikson 2016). The package uses AB pairs of slit-nodded spectra to sky subtract and then optimally extract the target flux. Wavelength calibration is carried out using both OH sky emission and telluric absorption from the A0V star observed immediately before or after the target. The A0V star also serves as a telluric standard, and when the target spectrum is divided by the A0V spectrum the target is corrected for both telluric absorption and the instrument profile. The final product of the PLP is a wavelength-calibrated spectrum of the target star, with flux in counts, and the corresponding signal-to-noise spectrum. While reduced spectra were produced across the wavelength range ∼1.45–2.52 μm, we only consider the K-band spectra in subsequent analyses given the low signal to noise of the H-band spectrum. This companion is brightest in K band (1.85–2.52 μm), although we find that some K-band orders are also unusable due to a low signal-to-noise ratio (S/N).

Instrumental resolution and rotation both produce line broadening, and these two sources of broadening are degenerate. It is thus important to accurately measure the resolution in order to accurately measure $v\sin i$. Using observed standard star spectra, we first select four IGRINS orders spanning wavelengths 2.293–2.325, 2.236–2.267, 2.182–2.212, and 2.105–2.135 μm. For each of the six epochs of data, we use the molecfit routine, which simultaneously fits a telluric model and an instrumental profile defined by a single Gaussian kernel, to the spectrum (Kausch et al. 2015; Smette et al. 2015). We take the median of these 24 instrumental resolution measurements as the instrumental resolution to use in our analysis, and the standard deviation to be the uncertainty. From this analysis we find R = 48,599 ± 1514. We also check whether the resolution changes significantly within an order. We select one of the epochs taken on UT 2020 February 04, and divide each of the four orders into five parts (each ∼400 pixels across). We find that the resulting instrumental resolution values within each order are consistent with each other and with the global resolution measurement.

In the wavelength-calibrated and telluric-corrected output spectra from the IGRINS reduction routine, we remove artifacts from strong sky lines that manifested as spikes in the data. In addition, we find that the short wavelength end of each spectrum contains less flux as the instrument blaze falls off. We cut the leftmost 20–70 pixels off of each reduced spectral order, where the cutoff value grew with increasing order number (decreasing wavelength).

With these spectra of HD 106906 b, we want to measure the amount of line broadening due to the rotational velocity. To do so, we calculate the "data" cross-correlation function (CCF) between the observed spectrum and a model atmosphere, where the model has been broadened to the instrumental resolution. We use an atmospheric model from the Sonora model grid (Marley et al. 2018, 2021; C. Morley et al. 2021, in preparation). These models are calculated assuming that the atmosphere is in radiative–convective and chemical equilibrium, following the approach of Marley et al. (1999), Saumon & Marley (2008), and Morley et al. (2012), with updated chemistry and opacities as described in (Marley et al. 2018; Marley et al. 2021, in preparation.). We assume T_eff = 1820 K and $\mathrm{log}(g)$ = 4.0 for HD 106906 b, following measurements of T_eff = 1820 ± 240 K and $\mathrm{log}({L}_{\mathrm{bol}}/{L}_{\odot })$ = −3.65 ± 0.08 using medium resolution spectra from VLT/SINFONI (Daemgen et al. 2017), and converting $\mathrm{log}({L}_{\mathrm{bol}}/{L}_{\odot })$ and system age to $\mathrm{log}(g)$ using hot-start evolutionary models (Burrows et al. 1997). The Sonora models generated for HD 106906 b have solar metallicity and solar carbon-to-oxygen ratio (C/O), and include silicate, iron, and corundum clouds with a sedimentation efficiency f_sed = 2 as described in Ackerman & Marley (2001).

We compare this "data" CCF to "model" CCFs, where each model is calculated by cross correlating a model atmosphere broadened by the instrumental line profile with that same model additionally broadened by some rotation rate and offset by a radial velocity (RV). We perform this comparison in a Bayesian framework using MCMC to fit for three free parameters: $v\sin i$, RV, and instrumental resolution. We use uniform priors on $v\sin i$ and RV. For the instrumental resolution, we use a Gaussian prior with a mean of 48,599 and standard deviation of 1514 to account for uncertainties on the measured resolution.

The log-likelihood function used in our MCMC framework is given by

$$\begin{eqnarray}&&\mathrm{log}L=\sum _{i=1}^{n}-0.5{\left(\displaystyle \frac{{m}_{i}-{d}_{i}}{{\sigma }_{i}}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where d is the "data" CCF, m is the "model" CCF, and σ_i is the CCF error at position i. We calculate uncertainties on the "data" CCF using the jackknife resampling technique. In this case, uncertainties are given by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{jackknife}}^{2}=\frac{(n-1)}{n}\displaystyle \sum _{j=1}^{n}{({x}_{j}-x)}^{2},\end{eqnarray} \tag{ 2 }$$

where n is the total number of samples. We define a sample as one epoch of data–there are six for HD 106906 b. x_j is the "data" CCF calculated using all epochs of data except the jth epoch, and x is the "data" CCF calculated using all epochs of data.

Before undertaking a joint fit of multiple orders to determine $v\sin i$, we first consider each order in K band individually. We compute "data" CCFs for 22 orders spanning wavelengths 1.85–2.42 μm (the first three spectral orders extended to longer wavelengths than those covered by our models), and determine the significance of the peak (if present) by calculating the ratio of the peak height to the standard deviation of the CCF outside the central peak. We exclude orders with peaks with less than 5σ significance. This cut leaves us with 15 orders running from 1.99–2.42 μm (excluded orders have significant telluric features and lower S/N spectra). We fit each order individually to get independent estimates for $v\sin i$, and find that they are consistent within their uncertainties. We then perform a joint fit, calculating a "data" CCF using all 15 orders and fitting models across that entire 1.99–2.42 μm wavelength range. The measured projected rotation rate for HD 106906 b is $v\sin i$ = 9.5 ± 0.2 km s⁻¹ (see Figures 1–3 for reference).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now consider how our modeling assumptions could impact the measured $v\sin i$. First we test our choice of T_eff and $\mathrm{log}(g)$. From Daemgen et al. (2017) we have T_eff = 1820 ± 240 K, and we converted $\mathrm{log}({L}_{\mathrm{bol}}/{L}_{\odot })$ to $\mathrm{log}(g)$ = 4.0 ± 0.5 using hot-start evolutionary models (Burrows et al. 1997). We take the 1σ errors on these values and generate four models: (1580 K, 3.5 dex), (1580 K, 4.5 dex), (2060 K, 3.5 dex), and (2060 K, 4.5 dex). We calculate new $v\sin i$ values with each of these models to test the possible impact of the measured uncertainties on our adopted T_eff and $\mathrm{log}(g)$, and found that all $v\sin i$ values were consistent with our original measurement at the ≤2σ level (see Table 1).

Table 1. Model Tests and Resulting $v\sin i$'s

Model	$v\sin i$
Original	9.53 ± 0.24 km s−1
1580 K, 3.5 dex	9.88 (+0.25 −0.24)
1580 K, 4.5 dex	8.85 ± 0.22
2060 K, 3.5 dex	9.89 (+0.23 −0.26)
2060 K, 4.5 dex	9.20 (+0.22 −0.21)
0.25× solar C/O	9.63 (+0.24 −0.23)
0.5× solar C/O	9.69 (+0.28 −0.24)
1.5× solar C/O	8.73 (+0.25 −0.20)
0.1×P	9.87 (+0.20 −0.22)
10×P	8.84 (+0.33 −0.28)

Download table as:  ASCIITypeset image

Another assumption we make when generating atmospheric models is a C/O. In our original model we assume a solar (0.54) C/O value, and here we test three additional ones: 0.25 × solar, 0.5 × solar, and 1.5 × solar. When we implement these models in our MCMC framework, we find that resulting $v\sin i$ values are consistent with the original value to ≤0.5σ for the subsolar C/Os, and differ by 2.3σ for the 1.5 × solar model (see Table 1). While not significant, this tentative offset in $v\sin i$ due to higher C/O suggests that abundance assumptions can become important for $v\sin i$'s given our measurement precision of 0.2 km s⁻¹.

Finally, we explore the impact that uncertainties in pressure broadening can have on the measured rotational velocity. Pressure broadening is a degenerate effect along with instrumental broadening and rotational line broadening–higher-pressure broadening with the same instrumental resolution leads to less rotational line broadening and a correspondingly smaller $v\sin i$. To test our pressure broadening assumptions, we run two models with modified molecular opacities, where molecular cross sections were 10× and 0.1× the actual pressure for the whole profile. This simulates a scenario where pressure broadening parameters that are used to create the molecular cross sections are off by an order of magnitude. Collision-induced opacity of hydrogen and helium is treated separately for all models, using the standard pressure for each layer. When we recalculate $v\sin i$ values using these new models, we find that for the 0.1× model (which underpredicts the amount of pressure broadening) the resulting rotation rate is only 1.0σ away from the original value, and the 10× model (which overpredicts pressure broadening) produces a $v\sin i$ that is 1.7σ lower (see Table 1). Both values are consistent with the original rotation rate measurement.

Periodic features in substellar light curves can be produced by cloud patchiness or planetary-scale waves, which manifest as longitudinal bands produced by zonal circulation (Apai et al. 2017, 2021). Lower-gravity objects typically have higher variability amplitudes and higher intrinsic-variability rates (Metchev et al. 2015; Vos et al. 2019). This observed variability appears to be impacted by viewing geometry–more highly inclined objects (i.e., closer to pole on) have more attenuated brightness changes (Vos et al. 2017).

The photometric rotation period for HD 106906 b was published by Zhou et al. (2020). The authors used the Hubble Space Telescope Wide Field Camera 3 (WFC3) near-IR channel in time-resolved direct imaging mode to observe HD 106906 b in three bands: F127M, F139M, and F153M. Using techniques such as two-roll differential imaging and hybrid point-spread function modeling yielded ∼1% precision in the light curves across all three bands. Fitting the light curve with a sinusoid results in a period of 4.1 ± 0.3 hr and an amplitude of 0.49% ± 0.12%.

Zhou et al. (2020) present several caveats to the rotation period measurement. First, they find only marginal evidence of variability, with a significance of 2.7σ. Because of this tentative detection, when fitting the light curves the authors applied a strict sinusoidal shape to the photometric modulations, and could not investigate whether the light curve could have multiple peaks. Previously, Apai et al. (2017) found that for 3 L/T transition brown dwarfs with high signal-to-noise data and extremely long baselines (>1 yr), the power spectra of their light curves produced peaks at both the full rotation period of the object as well as half the rotation period. A more recent analysis of long-baseline high-S/N photometry of Luhman 16 A and B found a similar result, with peaks at both the full and half rotation period (Apai et al. 2021). While this raises the question of whether the detected 4 hr rotation period of HD 106906 b is the full or half rotation period, both higher quality and theoretical light curves of brown dwarfs show that the full rotation period dominates the signal in the power spectrum for a given light curve (Zhang & Showman 2014; Apai et al. 2017, 2021). In addition, periodicity in the light curves of Jupiter and Neptune correspond to their full rotation periods (Karalidi et al. 2015; Simon et al. 2016; Ge et al. 2019). We thus assume that the full rotation period of HD 106906 b is 4.1 ± 0.3 hr.

Another caveat to this rotation period measurement is that photometric modulations for HD 106906 b are only detected in the bluest band (F127M), and not in the other two bands. However, for the majority of substellar objects, rotational modulations are wavelength dependent and have higher amplitudes at shorter (bluer) wavelengths (e.g., Apai et al. 2013; Yang et al. 2015; Zhou et al. 2016, 2019). Assuming a similar wavelength dependence for the photometric modulations in HD 106906 b as measured for 2M1207 b (Zhou et al. 2016), the closest spectral type young companion analog to HD 106906 b with detected modulation, the authors predict that the modulation amplitude for the redder bands would have been too small for them to observe. The detection of modulation in only the bluest band is therefore consistent with low overall amplitude variability and wavelength dependent modulations.

Given the wide projected separation of HD 106906 b (737 au), detecting orbital motion and placing constraints on orbital parameters requires long-baseline precision astrometry (e.g., Bowler et al. 2020). Recently, Nguyen et al. (2021) detected orbital motion of this companion using 14 yr of astrometry measurements between 2004 and 2017. All data were taken with HST using a combination of the Advanced Camera for Surveys, the Space Telescope Imaging Spectrograph, and the WFC3. By cross registering background star locations in the HST images with the Gaia astrometric catalog, Nguyen et al. (2021) calculated astrometry to high precision (at the subpixel level). To obtain constraints on the orbital parameters of HD 106906 b, the authors used the open-source Python package orbitize! ³ (Blunt et al. 2020) to perform an orbit fit to the assembled astrometric measurements. The line-of-sight orbital inclination i_o from these fits is ${56}_{-21}^{+12}$°. In this paper we use the posterior distribution shown in Figure 10 of Nguyen et al. (2021) when incorporating i_o​ in our analyses.

Scattered light imaging of the debris disk in the HD 106906 system with both GPI and SPHERE constrained the morphology of the disk (Kalas et al. 2015; Lagrange et al. 2016). In Lagrange et al. (2016), the disk detection with SPHERE was modeled using the GRATER code (Augereau et al. 1999) as an optically thin, inclined ring centered on the host binary. The authors assumed a dust density distribution that peaks at radius r₀ and has a power-law slope of α_in inside of r₀ and α_out outside of r₀. In addition to r₀ and α_out, fitted model parameters include the inclination of the disk i_d , the position angle (PA), a scaling factor to match the total flux of the disk, and the Henyey-Greenstein coefficient g, which quantifies how anisotropic the scattering is. With this modeling, the authors find a disk inclination i_d = 853 ± 01. In Kalas et al. (2015), the authors estimate a disk inclination of i_d ∼ 85° by assuming the disk is circular and translating the disk aspect ratio from their fitted semimajor and minor axes to line-of-sight inclination. Because it is unclear whether the disk is rotating in a prograde or retrograde fashion, there is a degeneracy in i_d and Ω_d pairs, where Ω_d is the PA of the ascending node. Thus angles (i_d = 85°, Ω_d = 104°) and (i_d = 95°, Ω_d = 284°) are equally likely. In this paper, we assume a bimodal distribution for i_d , with i_d = 853 ± 01 and i_d = 947 ± 01 defining the two Gaussian distributions.

We combine P_rot and our measurement of $v\sin i$ to determine the line-of-sight spin axis inclination of the companion i_p . However, simply computing the inclination as:

$$\begin{eqnarray}&&\sin i=\left(\displaystyle \frac{v\sin i}{2\pi R/{P}_{\mathrm{rot}}}\right)\end{eqnarray} \tag{ 3 }$$

does not account for correlations between relevant parameters (Masuda & Winn 2020). For example, v and $v\sin i$ are not statistically independent given that $v\sin i$ is always less than v. We therefore follow the method described in detail in Masuda & Winn (2020) and summarized for the application to HD 106906 b below.

Given v as the equatorial rotational velocity and u = $v\sin i$ as the projected rotation rate, we have the following two likelihood functions:

$$\begin{eqnarray}&&{L}_{v}(v)=p({d}_{v}| v)\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{L}_{u}(u)=p({d}_{u}| u)\end{eqnarray} \tag{ 5 }$$

where d_v and d_u are the data sets from which these likelihood functions are calculated. In our case, L_u (u) is the probability distribution for $v\sin i$ that we determined from our high-resolution spectra, a Gaussian with peak location 9.5 km s⁻¹ and standard deviation 0.2 km s⁻¹. L_v (v) is the probability distribution for v, which we calculate using v = 2π R/P_rot. For the radius R we calculate the effective blackbody radius:

$$\begin{eqnarray}&&R=\sqrt{\displaystyle \frac{L}{4\pi {\sigma }_{b}{T}_{\mathrm{eff}}^{4}}},\end{eqnarray} \tag{ 6 }$$

where L is the bolometric luminosity log (L_bol/L_⊙) = −3.65 ± 0.08, σ_b is the Stefan–Boltzmann constant, and T_eff is the effective temperature T_eff = 1820 ± 240 K (Daemgen et al. 2017). This yields a radius of ${1.49}_{-0.45}^{+0.37}$ R_Jup. We produce a probability distribution for v by incorporating uncertainties on P_rot, log (L_bol/L_⊙), and T_eff in a Monte Carlo fashion.

With L_u (u) and L_v (v) in hand, Masuda & Winn (2020) specify two key assumptions:

1. d_v and d_u are independent, so the likelihood function for D = {d_v , d_u } is separable:

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{L}_{vu}(v,u)=p(D| v,u) & = & p({d}_{v}| v,u)\,=\\ p({d}_{v}| v)p({d}_{u}| u) & = & {L}_{v}(v){L}_{u}(u);\end{array}\end{array}\end{eqnarray} \tag{ 7 }$$

2. v and i are a priori independent, which means that the prior P_vi (v, i) is separable:

$$\begin{eqnarray}&&{P}_{{vi}}(v,i)={P}_{v}(v){P}_{i}(i)\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&p(v| i)={P}_{v}(v);p(i| v)={P}_{i}(i).\end{eqnarray} \tag{ 9 }$$

Given these assumptions, the posterior distribution function (PDF) for $\cos i$ can be written as:

$$\begin{eqnarray}&&p(\cos i| D)\propto {P}_{{\rm{c}}{\rm{o}}{\rm{s}}{\rm{i}}}(\cos i)\int {L}_{v}(v){L}_{u}(v\sqrt{1-{\cos }^{2}i}){P}_{v}(v)dv.\end{eqnarray} \tag{ 10 }$$

where P_cosi(cosi) is uniform between 0 and 1 and the prior on the rotation rate P _v (v​) is uniform from 0 to break-up speed.

Converting this PDF in $\cos i$ to a PDF in i yields the distribution shown in Figure 4. We note that the posterior distribution for i_p is bimodal and symmetric around 90° simply because we do not know whether this spin axis vector (which has directionality) is pointed toward us or away from us. Therefore the mode and 68% highest probability density interval (HPDI) of i_p is 14° ± 4° for i_p < 90°, and 166° ± 4° for i_p > 90°.

Zoom In Zoom Out Reset image size

Our goal is to measure the true 3D angles between the three angular momentum vectors in this system–the companion spin axis, the companion orbit normal, and the disk normal. These angles are given by:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{op}}={\cos }^{-1}(\cos {i}_{p}\cos {i}_{o}+\sin {i}_{p}\sin {i}_{o}\cos ({{\rm{\Omega }}}_{o}-{{\rm{\Omega }}}_{p}))\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{dp}}={\cos }^{-1}(\cos {i}_{p}\cos {i}_{d}+\sin {i}_{p}\sin {i}_{d}\cos ({{\rm{\Omega }}}_{d}-{{\rm{\Omega }}}_{p}))\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{od}}={\cos }^{-1}(\cos {i}_{d}\cos {i}_{o}+\sin {i}_{d}\sin {i}_{o}\cos ({{\rm{\Omega }}}_{o}-{{\rm{\Omega }}}_{d}))\end{eqnarray} \tag{ 13 }$$

where Ψ_op is the true companion obliquity, Ψ_dp is the true spin-disk mutual inclination, and Ψ_od is the true orbit-disk mutual inclination. The PAs Ω_o , Ω_d , and Ω_p measure how the orbit, disk, and companion spin axis, respectively, are oriented on the sky plane. The nodal angle Ω_p is unknown.

The difference between the line-of-sight inclination of the companion spin axis i_p and that of the orbit normal i_o yields a lower limit on the true deprojected obliquity Ψ_op (Bowler et al. 2017):

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{op}}\gt | {i}_{p}-{i}_{o}| .\end{eqnarray} \tag{ 14 }$$

Similarly:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{dp}}\gt | {i}_{p}-{i}_{d}| \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{od}}\gt | {i}_{o}-{i}_{d}| .\end{eqnarray} \tag{ 16 }$$

Figure 5 shows the posteriors for ∣i_p − i_o ∣, ∣i_p − i_d ∣, and ∣i_o − i_d ∣. We also plot a random distribution in black for comparison, where i_p , i_o , and i_d are all drawn from uniform distributions in $\cos i$. We find that the 68% HPDI for ∣i_p − i_o ∣ lies between [32°, 119°]. For ∣i_o − i_d ∣ the 68% HPDI is [16°, 48°], and for ∣i_p − i_d ∣ we have the tightest 68% HPDI of [69°, 83°].

Zoom In Zoom Out Reset image size

Since these projected angles are all lower limits on the true 3D angles, we can say that our results tend to favor more "misaligned" orientations (defined here as mutual inclinations of 20°–180°). A schematic illustration of the line-of-sight architecture of the system is shown in Figure 6.

Zoom In Zoom Out Reset image size

We now calculate full probability distributions for all Ψ's using Equations (11)–(13). For Equations (11) and (12) we assume Ω_p is randomly drawn from a uniform distribution between 0 and 2π. For Equation (13), both Ω_o and Ω_d have been measured: Ω_o = ${99}_{-28}^{+26}$° or ${279}_{-29}^{+25}$°, and Ω_d = 1044 ± 03 or 284° ± 03 (Kalas et al. 2015; Lagrange et al. 2016; Nguyen et al. 2021). Figure 7 shows the resulting probability distributions for Ψ_op, Ψ_dp, and Ψ_od, along with a random mutual inclination distribution Ψ_random in black for reference.

Zoom In Zoom Out Reset image size

Not knowing Ω_p along with poor constraints on HD 106906 b's orbital elements (a consequence of the companion's wide orbital separation) leads to broad posterior distributions for the true deprojected angles Ψ_op and Ψ_od. By comparison, Ψ_pd is remarkably well constrained. All of these posteriors are bimodal, reflecting symmetries across 90°. For each Ψ posterior we calculate the mode and 68% HPDI for each half of the distribution below and above 90°. We find that the true companion obliquity Ψ_op is ${55}_{-16}^{+22}$° or ${125}_{-22}^{+16}$°, the true spin-disk angle Ψ_dp is ${84}_{-8}^{+6}$° or ${96}_{-6}^{+8}$°, and the true mutual inclination between the orbit and disk normals Ψ_od is ${39}_{-15}^{+20}$° or ${141}_{-20}^{+15}$° (Table 2). Each of these deprojected mutual inclinations deviates distinctly from the geometric prior–both Ψ_op and Ψ_od favor angles away from 90°, which is where the geometric prior peaks, and while Ψ_dp peaks around 90° it is a much tighter constraint than a random distribution.

Table 2. Measured Parameters

Parameter	Measured Value	Ref
${v}_{p}\sin {i}_{p}$	9.5 ± 0.2 km s−1	This work
Prot,p	4.1 ± 0.3 hr	Zhou et al. (2020)
ip	14° ± 4° or 166° ± 4°	This work
io	${56}_{-21}^{+12}$°	Nguyen et al. (2021)
id	85.3 ± 0.1 or 94.7 ± 0.1 deg	Kalas et al. (2015), Lagrange et al. (2016)
∣io − ip ∣	(32°, 119°)	This work
∣ip − id ∣	(69°, 83°)	This work
∣id − io ∣	(16°, 48°)	This work
Ψop	${55}_{-16}^{+22}$° or ${125}_{-22}^{+16}$°	This work
Ψdp	${84}_{-8}^{+6}$° or ${96}_{-6}^{+8}$°	This work
Ψod	${39}_{-15}^{+20}$° or ${141}_{-20}^{+15}$°	This work

Note. The three i inclinations presented here are all along our line of sight. The angle i_p is symmetric about 90° due to the fact that we do not know whether this spin angular momentum vector is pointing toward us or away from us. The angle i_d has two solutions because it is unclear whether the disk is rotating prograde or retrograde, so there are two combinations of i_d and Ω_d , where Ω_d is the PA of the ascending node, that are equally likely. The line-of-sight mutual inclinations ∣i_o − i_p ∣, ∣i_p − i_d ∣, and ∣i_d − i_o ∣ are lower limits on the true deprojected angles Ψ_op, Ψ_dp, and Ψ_od. Here we quote the mode and 68% highest probability density intervals for these angles.

Download table as:  ASCIITypeset image

We quantify the probability that each Ψ distribution yields an "aligned" state, which we define as Ψ ∈ (0°, 20°), or a "misaligned" state where Ψ ∈ (20°, 180°). The Bayesian odds ratio is p(m∣D)/p(a∣D), where p(m∣D) is the probability of misaligned state m given data D, and p(a∣D) is the probability of an aligned state a. The former probability is the integral of the posterior distribution p(Ψ∣D) (the colored histograms in Figure 7) from 20°–180°, while the latter probability is obtained by integrating the posterior distribution from 0°–20°.

We compute these integrals over the posterior distributions for all three mutual inclinations Ψ_op, Ψ_dp, and Ψ_od, as well as for the geometric prior distribution which is uniform in $\cos {{\rm{\Psi }}}_{\mathrm{random}}$ (Figure 7, black curve). In the absence of any data, the geometric prior yields an odds ratio favoring misalignment to alignment at 32:1 (2.2σ). The odds ratios for the true companion obliquity Ψ_op and the orbit-disk mutual inclination Ψ_od are only slightly larger, 39:1 (2.3σ) and 40:1 (2.3σ), respectively. Nevertheless, the fact that the shapes of both posteriors differ significantly from that of the prior indicates that the odds favoring misalignment for both Ψ_op and Ψ_od are not entirely driven by the prior. We can say there is tentative evidence that both of these mutual inclinations are large. By comparison, the true spin-disk angle Ψ_dp has an odds ratio that favors misalignment much more strongly than does our random prior–97081:1 (4.4σ). We can say with confidence that this system contains an edge-on disk and a planet spinning nearly pole on; even allowing for sky-plane degeneracy, our measurements indicate that the companion spin axis and disk normal are strongly misaligned.