© The Authors 2024

1 Introduction

Various types of exotic planetary systems have been discovered over the last 30 years, including hot Jupiters, planets on extremely eccentric orbits, planets orbiting binary stars, planets orbiting white dwarfs, neutron stars, and others. A broad range of architectures has also been revealed by measurements of the host star’s obliquity, namely, the angle between the star’s spin axis and a planet’s orbital axis (or at least its sky projection). Thanks to many research groups pursuing these measurements, usually by means of the Rossiter–McLaughlin (RM) effect, we now have a sample of nearly 200 measurements including prograde, polar, and retrograde orbits (see, e.g., Albrecht et al. 2022, for a review).

Here, we present observations of 19 transiting exoplanets aiming to detect the RM effect and measure the sky-projected stellar obliquities. Together with the sample drawn from the literature, we also seek to investigate four issues regarding the distribution of obliquities.

(i) For most systems in the sample, the projected obliquity (λ) is known but not the obliquity itself (ψ) because of missing information about the inclination of the stellar rotation axis with respect to the line of sight (i_★). Using the subset of 57 systems for which ψ has been measured, Albrecht et al. (2021) found evidence for a preponderance of perpendicular planets, namely, a peak in the obliquity distribution near 90°. Following up on this result, Siegel et al. (2023) and Dong & Foreman-Mackey (2023) applied more advanced statistical techniques that allowed the entire sample to be used, and did not find evidence for a peak. With an enlarged sample, we wanted to revisit this issue and see whether there is a particular category of planets where such a peak exists.

(ii) Three theories have been proposed to explain the existence of hot Jupiters: in situ formation, disk-driven migration, and high-eccentricity or tidally-driven migration (see Dawson & Johnson 2018, for a review). The latter process would not only increase the orbit’s eccentricity, but might also raise its inclination relative to the star’s equatorial plane. We might therefore expect high eccentricities and obliquities to be statistically associated, and evidence for such an association has been reported (Rice et al. 2022). After noticing some errors in this earlier study, we decided to revisit the issue with an enlarged sample.

(iii) Stars with T_eff ≲ 6250 K that host hot Jupiters tend to have especially low obliquities. Indeed, the most precise such measurements (σ_λ < 2°) show an obliquity dispersion less than one degree – smaller than the 6.2° obliquity of the Sun relative to the Solar System’s invariable plane. Highlighting this result, Albrecht et al. (2022) interpreted the very low obliquities of hot Jupiters around cool stars as evidence for tidal obliquity damping. If this is true, by continuing to perform precise measurements, we might gain a better understanding of the evolution of hot Jupiters and rates of tidal dissipation. With the enlarged sample, we seek to see if this trend still holds and include the numerous (albeit less precise) obliquity measurements in the determination of the dispersion.

(iv) The first five measurements of λ for stars with compact multiple-transiting planetary systems (“multis”) were all consistent with low obliquities (Albrecht et al. 2013), following the blueprint of the Solar System. Since then, several exceptions have been discovered (Hjorth et al. 2021; Huber et al. 2013; Chaplin et al. 2013). The misalignments in these cases are suspected of being caused by the effects of an outer massive companion. We wanted to assess the obliquity distribution of the multis, and see whether all of the misaligned multis have outer companions.

The paper is structured as follows. Section 2 describes the general characteristics of the new data. Section 3 presents our general approach to analyzing the data. Section 4 discusses the issues and questions posed above and Section 5 summarizes our conclusions. The particulars for each of the 19 systems are discussed in Appendix A.

2 Observations

To measure the projected obliquity of each host star, we performed high-resolution optical spectroscopy over a time range (typically 4–8 hours) spanning a planetary transit and tried to detect the RM effect. Table 1 gives the names of the 19 systems, the names of the telescopes and instruments employed, and some key observational characteristics. The telescopes and instruments are further described in Table 2. Tables 3 and 4 give some of the basic parameters of the host stars and their planets. All of the apparent radial velocity (RV) measurements are available on Zenodo (see Data availability). These spectroscopic data were supplemented by the best available transit light curves, as described below.

2.1 Spectroscopy

Seven systems were observed with the Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2021) that can be fed by any of the four Very Large Telescopes at paranal Observatory, Chile. Another seven systems were observed with the FIber-fed Echelle Spectrograph (FIES; Frandsen & Lindberg 1999; Telting et al. 2014) mounted on the Nordic Optical Telescope (NOT; Djupvik & Andersen 2010) located on the Roque de los Muchachos, La Palma, Spain. For five systems, we used the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N; Cosentino et al. 2014) mounted on the Telescopio Nazionale Galileo (TNG), also on the Roque de los Muchachos. For one system, we used the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) mounted on the 10-m Keck-1 telescope on Mauna Kea, Hawai’i, USA.

The data from ESPRESSO and HARPS-N were reduced and radial velocities (RVs) were obtained using the standard data reduction software (DRS; Lovis & Pepe 2007; Pepe et al. 2021; Dumusque et al. 2021). The FIES data for all but two systems were reduced using FIEStool¹, and RVs were extracted using FIESpipe² following an approach similar to that described by Zechmeister et al. (2018). The FIES observations of WASP-136 and WASP-186 were reduced and RVs extracted following the procedure described by Gandolfi et al. (2015). To account for any instrumental RV drift of the FIES spectrograph, we observed the spectrum of a Thorium-Argon (ThAr) lamp in between each science exposure. The HIRES spectra were obtained with the iodine cell in the light path and the RV extraction followed the standard HIRES forward modeling pipeline of the California Planet Search (CPS; Howard et al. 2010) using the methodology of Butler et al. (1996).

For some systems, in addition to the apparent RVs we modeled the shape of the Cross Correlation Function (CCF) of the absorption lines. For the ESPRESSO and HARPS-N data, we used the CCFs that are automatically produced by the standard reduction pipelines. For the FIES data, we derived CCFs using a template spectrum based on the most suitable model atmosphere from the library of Castelli & Kurucz (2003), as implemented in FIESpipe.

2.2 Photometry

Light curves of all 19 of the systems are available from the databases of either the K2 mission (Howell et al. 2014) or the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015). For the K2 data, we used the light curves from the EVEREST pipeline (Luger et al. 2016, 2018). The TESS data were downloaded and extracted using the lightkurve package (Lightkurve Collaboration 2018), and the RegressionCorrector was applied to reduce the systematic variations due to scattered light. Our parametric model for each system was based on jointly fitting the photometric and spectroscopic data. In some cases, we supplemented the space-based photometry with ground-based photometry to refine the transit ephemerides.

3 Data analysis

3.1 Photometric data

The photometric data provide tight constraints on the orbital period (P), the time of conjunction of a transit chosen to be the reference epoch (T₀), the planet-to-star radius ratio (R_p/R_★), the cosine of the orbital inclination (cos i), and the ratio of the semi-major axis and stellar radius (a/R_★). Uniform prior probability distributions were adopted for these parameters.

To model the stellar limb darkening function, we used a standard quadratic law and allowed the sum of the two coefficients (the center-to-limb intensity ratio) to be an adjustable parameter, while holding the difference fixed at the value obtained from the tables of Claret et al. (2013); Claret (2018) as appropriate for the photometric bandpass and the star’s effective temperature (T_eff), surface gravity (log 𝑔), and metallicity ([Fe/H]). See Table 5 and Table A.1 for details. The sum of coefficients was subjected to a Gaussian prior with a width of 0.1 and a central value from the aforementioned tables.

In cases for which a nonzero orbital eccentricity has been reported (HD 118203 b, K2-261 b, K2-287 b, WASP-148 b, and WASP-186 b), we applied Gaussian priors on e and the argument of periastron (ω) using the values given in Table 3. For the multis, we applied Gaussian priors on the parameters of all planets apart from the transiting planet that was the target of our observations. The main effect of these planets on our model is to change the slope of the RV time series on the night of the transit that was observed spectroscopically. When two planets are transiting simultaneously, the shape of the transit light curve is also affected. Simultaneous transits occasionally occurred (for TOI-1130 and TOI-451A) within the long intervals spanned by the photometric data, but did not occur during any of our 19 spectroscopic observations.

We modeled the light curves using the batman package (Kreidberg 2015), which is based on the equations of Mandel & Agol (2002). In addition to accounting for the transit and white noise in the measurements, we included a Gaussian process (GP) using celerite (Foreman-Mackey et al. 2017). We generally employed a Matérn-3/2 kernel, which is characterized by two hyperparameters: the amplitude (A) and the timescale (τ). For the time series with 30 minute sampling, evenly spaced model light curves were created and integrated over to mimic a cadence of 2 min.

3.2 Stellar rotation

The analysis of the RM effect is aided by knowledge of the stellar rotation velocity or period. We tried to use the available spectroscopic and photometric data to measure these quantities independently of the transit data, as described below.

3.2.1 Projected rotation velocity

The amplitude of the RM effect scales with v sin i_★, the product of the stellar rotation speed and the sine of the line-of-sight inclination of the stellar spin axis. Therefore, knowing v sin i_★ helps us plan observations and analyze the results. When the RM effect is detected with a high signal-to-noise ratio (S/N), v sin i_★ can be determined precisely from the RM data. Indeed, for stars with a high obliquity, the effects of differential rotation might be detectable (Gaudi & Winn 2007; Cegla et al. 2016).

In more typical cases in which the RM effect is detected with a modest S/N, it is useful to determine v sin i_★ by modeling the spectral line broadening outside of transits, and thereby gain the option of applying a prior constraint on v sin i_★ when modeling the RM effect. For this purpose, we used the broadening function (BF) method of Rucinski (1999) as implemented in FIESpipe, which is based on an appropriate model atmosphere from Castelli & Kurucz (2003). The BFs were created from high S/N orders from out-of-transit spectra, where we visually inspected the result to ensure that the BF created from a given order was well-behaved. This typically meant selecting central orders typically with wavelengths ~5000–6000 Å. Orders containing telluric lines were excluded. For a given epoch we stacked the BFs from the individual orders to create the highest possible S/N BF for that epoch. We applied this method not only to the FIES spectra, but also the ESPRESSO and HARPS-N spectra. To extract v sin i_★ we fitted the theoretical rotational BF from Kaluzny et al. (2006, Eq. (2)). The model includes the rotational broadening parameters along with a parameter representing non-rotational broadening (in most cases dominated by the finite instrumental resolution).

Table 6 gives the results. The tabulated v sin i_★ is the median of the best-fit results of analyzing spectra from different epochs, and the tabulated uncertainty is the standard deviation between those epochs, which in most cases are unrealistically small. This is especially true for observations where we only have few out-of-transit observations. For systems where we only had two out-of-transit spectra, we used the first ingress as a third spectrum. Our measurements are generally in agreement with previously published results. The most significant deviations are seen for the most slowly-rotating stars, where the details of extracting line profiles and modeling non-rotational broadening are most important. When modeling the RM effect, we imposed Gaussian priors based on our v sin i_★ determinations, with a minimum width of σ = 0.5 km s⁻¹ to account for possible systematic errors.

To test whether our choice of the width of the uncertainty interval in the v sin i_★ prior would significantly affect our results for the projected obliquity, we selected four systems in our sample for which we expected the result of lambda to be most sensitive to v sin i_★ (low impact parameter, incomplete transit coverage). We then performed the same analysis described in Appendix A. However, with an uncertainty in v sin i_★ of 1.0 km s⁻¹. In all four cases, we find that the highest probability value for λ changes only by a small fraction of its uncertainty interval, relative to the 0.5 km s⁻¹ case.

3.2.2 Photometric rotation period

We attempted to measure the rotation period (P_rot) by seeking periodicities in the K2 and TESS light curves. The combination of the rotation period, v sin i_★, and R_★ can be used to constrain the stellar inclination (Masuda & Winn 2020), which in turn can be used in combination with i and λ to determine the “true” or three-dimensional obliquity.

We used the method based on the auto-correlation function (ACF) described by McQuillan et al. (2014). After removing the data affected by transits, we searched for peaks in the ACF as a function of lag. When multiple peaks were seen, we fitted a linear function between peak number and lag, the slope of which is an estimate of the rotation period (e.g., Hjorth et al. 2021).

Table 6 gives the results. Empty entries mark the cases for which no period could be determined. For two systems, TOI-813 and TOI-892, peaks were seen in the ACF (corresponding to periods of 10.1 ± 0.6 d and 2.76 ± 0.11 d, respectively), but because the statistical significance was weak, the results do not appear in Table 6. For WASP-50, a series of strong peaks seemed to imply a period of 7.2 ± 1.8 d, which contrasts with the value of 16.3 ± 0.5 d reported by Gillon et al. (2011). They also noted that the true rotation period could be 32.6 days if there are similar spot patterns on opposite stellar hemispheres. We folded the light curve using trial periods of 16.3 d or 32.6 d, but neither choice seemed compelling. Evidently, the variability of WASP-50 is complex, which is why Table 6 contains no entry for WASP-50.

We note that the estimated rotation periods of WASP-136 and WASP-186, 5.4 ± 0.5 d and 4.8 ± 1.9 d, are both close to the orbital period of the transiting planets. These systems might have been driven into spin-orbit synchronization by tidal interactions between the planet and star (e.g., Albrecht et al. 2012; Brown 2014; Penev et al. 2018).

3.3 Rossiter–McLaughlin effect

The RM effect is the distortion of a star’s absorption lines that arises from the combination of stellar rotation and a transiting planet. During a transit, the planet prevents the light from a portion of the photosphere from contributing to the disk-integrated spectral line. The missing spectral component causes the velocity profile of the spectral lines to have a slightly reduced flux at velocities centered on the rotational Doppler shift of the star at the “subplanetary” point, namely, the point on the star directly behind the center of the planet. Observations of the RM effect can therefore be used to determine the time series of subplanetary velocities, or equivalently, the trajectory of the transiting planet relative to the stellar rotation axis.

Neglecting differential rotation, the RM effect is mainly a function of λ, v sin i_★, R_p/R_★, and the transit impact parameter b ≡ d/R_★ cos i, where d is the star-planet distance at the time of conjunction. For precise modeling it is also necessary to take into account limb darkening and non-rotational spectral line broadening. To model non-rotational broadening due to the star, we used the Gray (2005) model with parameters for the velocity spread generated by microturbulence (ξ_★) and macroturbulence (ζ_★). We included a separate term for instrumental broadening, another source of non-rotational broadening. We neglected the effects of differential rotation and the convective blueshift, which for solar-type stars are not expected to produce detectable effects given the quality of our data (and indeed our modeling did not uncover any evidence for those effects, with one possible exception; see Appendix A.12.1). We did not model the effects of starspots, flares, or pulsations, since there was no evidence for such phenomena in our data.

The RM effect has been analyzed as a distortion of individual spectral lines (e.g., Albrecht et al. 2007), as distortion of the spectral cross-correlation function (e.g., Cegla et al. 2016), or as an “anomalous RV” obtained by an ordinary radial-velocity extraction code to a spectrum affected by the RM effect (e.g., Queloz et al. 2000). We used a combination of these approaches, depending on the system, as described in Appendix A.

We fitted the spectral CCFs with a parameterized model that incorporates the RM effect, using a code first described by Knudstrup & Albrecht (2022) and Knudstrup et al. (2023). In this model, the stellar disk is discretized (with a radius of 100 pixels) and a spectrum is assigned to each pixel according to the model parameters for limb darkening, stellar rotation, and non-rotational broadening. A disk-integrated spectral line is created by summing the spectra of all pixels not concealed by the planet. The disk-integrated line profile is then compared to the observed spectral CCF.

An alternative is to fit the time series of subplanetary velocities derived from the residuals between the out-of-transit CCF and the in-transit CCFs. The subplanetary velocity is calculated as $$v_p(t)=v\sin i_{\star }x(t),$$(1)

where x(t) is the distance from the projected rotation axis divided by the stellar radius. Once the orbital parameters are chosen, x(t) can be calculated. It is also instructive to calculate x at the times of ingress (x₁) and egress (x₂) (Albrecht et al. 2011): $$\begin{array}{c}{x}_1=\sqrt{1-b^2}\cos \lambda -b\sin \lambda ,\\ x_2=\sqrt{1-b^2}\cos \lambda +b\sin \lambda .\end{array}$$(2)

Therefore, the very nearly linear trajectory of the planet can be calculated as a function of time for given values of t, T₀, the transit duration (T₄₁), b, and λ.

The simplest and most common approach is to fit a time series of anomalous RVs. In the model, the anomalous RV is computed as a function of vp by taking into account the size of the planet’s shadow, the limb darkening function, and the response of the radial-velocity extraction code to the RM distortion; for this purpose we used the code described by Hirano et al. (2011).

The different ways to sense the RM effect are illustrated in Fig. 1 for the case of WASP-50. For all 19 systems in our sample, we fitted the time series of anomalous RVs (the “RV-RM” method). For some systems, we also analyzed the distortions of the CCF or fit the time series of subplanetary velocities. Multiple methods were used to test for consistency in the results, and because for some systems it seemed plausible that the more sophisticated models would give more information about λ. In general, we expect line-profile modeling to be advantageous when (R_p /R_★)²v sin i_★ is larger than the non-rotational broadening (Albrecht et al. 2022).

We determined the best-fit parameters and their uncertainties via Markov Chain Monte Carlo (MCMC) sampling of the posterior probability distribution in parameter space. We used the emcee package (Foreman-Mackey et al. 2013). To ensure convergence we both inspected the chains visually and calculated the rank normalized R-hat³. As mentioned above, we always fit the transit light curve alongside the spectroscopic data.

When fitting the RVs and the line distortion directly, we used Gaussian priors for the macroturbulence and microturbulence parameters, and we adopted a fixed value for the instrumental broadening parameter of each spectrograph (see Table 2). The mean values for our Gaussian priors were obtained from the relationships presented by Doyle et al. (2014) for stars with T_eff ∊ [5200 K, 6400 K] and log 𝑔 ∊[4.0, 4.6], otherwise, Bruntt et al. (2010) for stars with T_eff ∊[5000 K, 6500 K] and log 𝑔 > 4.0. The width of the prior in all cases was 1 km s⁻¹.

We note that the properties of WASP-136, TOI-813, and HD 118203 are on the borderline of the range of properties for which the aforementioned relationships involving turbulent broadening are applicable; we applied them anyways. On the other hand, the properties of WASP-172 and LTT 1445A are far outside the applicable ranges. For WASP-172 (T_eff = 6900 ± 150 K), we expect a high macroturbulent velocity. We chose 6.0 km s⁻¹ and increased the width of the prior to 4.0 km s⁻¹ to reflect our uncertainty. For the microturbulent velocity, we adopted ξ_★ = 1.5 ± 1.5 km s^–1. In the other end of the spectrum, we expect the turbulent velocities of the M dwarf LTT 1445A to be smaller. We chose ζ = 1.0 ± 1.0 km s^–1 and ξ_★ = 1.0 ± 1.0 km s^–1. All of the priors are summarized in Table.

For each system, we performed two RV-RM analyses, one in which we adopted a uniform prior on v sin i_★, and one in which we applied a Gaussian prior on v sin i_★ based on the observed line broadening of the out-of-transit spectrum (Table 6) and a width of 0.5 km s^–1. When analyzing the line profiles or Doppler shadow, we always applied the Gaussian prior on v sin i_★.

When fitting the time series of subplanetary velocities, we did not simultaneously fit the light curve. Instead, we applied Gaussian priors on P, T₀, R_p/R_*, i, and a/R* based on an external fit to the light curve, and adopted uniform priors on λ and v sin i_★.

In some cases, to boost the S/N, we stacked the CCFs from different exposures in a way similar to Johnson et al. (2014), but with a code based more directly on the work of Hjorth et al. (2019). In these cases, we extracted the Doppler shadow as shown in Fig. 1, and then transformed the velocities by correcting for the slope in Eq. (2). The white line/slope in Fig. 1 is in this way used to bring each pixel directly beneath it in the Doppler shadow to be located at 0 km s^–1. We then collapsed the slope corrected shadow by summing along the ordinate (time stamps) for a given geometry (λ,v sin i_★,b). When the model parameters match reality, we expect to see a sharp peak at 0 km s^–1 in the summed spectrum, and we found the height of the peak by fitting a Gaussian to it.

Hjorth et al. (2019) did this for a dense grid in the parameter space of b, v sin i_★, and λ. For our purposes, since b was always tightly constrained by the light curve, we neglected the uncertainty in b and simply adopted the best-fit value from the RV-RM method. We then searched the two-dimensional parameter space of v sin i_★ and λ for the strongest peak; in practice this was done by fitting a two-dimensional Gaussian function. We did not run an MCMC when stacking the CCFs, nor did we fit the light curve simultaneously. The quoted uncertainties are taken as the widths (σ) of the fitted two-dimensional Gaussian and as such do not reflect the actual precision. The resulting values therefore only serves as an indication to whether the results in (λ,v sin i_★) from the RV-RM runs are in (qualitative) agreement with that seen in the CCFs.

Fig. 1 Rossiter–McLaughlin effect displayed in four ways, for the illustrative case of WASP-50. Top right: schematic of the planet transiting the stellar disk. The wedge represents the best-fit λ and its uncertainty. The color scheme is based on the star’s rotational Doppler shift; the same color scheme is used in the top left and lower left panels to convey the location of the planet on the star. Top left: spectral lines. The thick black curve is the average out-of-transit cross-correlation function (CCF) and the superimposed white curve is a CCF observed near mid-transit. Middle left: spectral line distortions. Shown are the results of subtracting five of the in-transit CCFs from the mean out-of-transit CCF. Solid curves show the data. Dashed curves are the best-fit Lorentzian profiles, with vertical lines marking their centers. The color coding is done according to the apparent velocity anomaly going from redshifted (in red) to blueshifted (in blue). Bottom left: Doppler shadow. Each row shows the spectral line distortion at a particular time. The white contours are amplitude levels of the distortion. Each small circle marks the best-fit subplanetary velocity (the Doppler shift of the portion of the star directly behind the planet’s center), with colors that convey the corresponding anomalous radial velocity shown in the lower right panel. Dark vertical lines mark the measured v sin i★, and the diagonal white line is the best-fit trajectory. Bottom right: anomalous radial velocities. The dark shaded area is the range of times when the planet’s silhouette is completely contained within the stellar disk.

3.4 Results

Table 7 gives our best estimates of λ and v sin i_★ for each system analyzed as described in Appendix A. For cases in which the rotation period could be measured, the table also includes the stellar inclination, i_★, calculated as per Masuda & Winn (2020), and the three-dimensional obliquity. Other key parameters are available on Zenodo (see Data availability).

4 Discussion

We return to the four questions posed in the introduction about the distribution of obliquities of stars with transiting planets. We combined the λ measurements presented in this paper (and listed in Table 7) with those measured previously. Specifically, we used the catalog assembled by Albrecht et al. (2022) with the help of the TEPcat catalog (Southworth 2011), and several other measurements that have appeared in the more recent literature (up until February 26th, 2024). Key parameters for these additional systems, along with the systems presented in this study are available on Zenodo (see Data availability).

Fig. 2 shows plots of λ and ψ as a function of T_eff, a/R_* and M_p/M_*.

As mentioned in their respective sections (Appendices A.7 and A.8) we excluded our results for LTT 1445Ab and TOI-451Ab in the ensemble analyses below. Furthermore as this study builds upon the sample presented in Albrecht et al. (2022) systems excluded in there are also excluded here (based on the criteria outlined in their Appendix A).

4.1 A possible preference for polar planets

Albrecht et al. (2021) investigated a sample of 57 systems for which ψ has been determined by combining measurements of the projected obliquity (λ) and the stellar inclination angle (i_★). Of the 19 values of ψ that are inconsistent with 0°, 18 are between 80° and 125°. They hypothesized that the misaligned systems show a preference for approximately polar orbits. Siegel et al. (2023) and Dong & Foreman-Mackey (2023) followed up by replicating the preceding work and analyzing the larger sample of systems, for which only λ was measured but not i_★. The larger sample did not show a significant preference for polar orbits, suggesting that the sample analyzed by Albrecht et al. (2021) was somehow biased toward polar orbits.

We revisited the issue with a sample that is about 40% larger. The current tally of systems with λ measurements is 205, of which ψ is known for 87. In the analysis that follows, when considering the ψ measurements, we decided to omit the 6 systems for which the gravity-darkening technique was employed, because parameter degeneracies make it difficult to measure ψ when it is near 0° or 180°. When considering only the λ measurements, we retained the six gravity darkening measurements.

First, we repeated the statistical tests performed by Albrecht et al. (2021) using the enlarged set of ψ measurements. We selected the 27 systems for which cos ψ < 0.75 and computed two statistics that quantify clustering: the dispersion around cos ψ = 0 and the standard deviation relative to the mean. We then used Monte Carlo simulations to find the probability of obtaining clustering statistics at least as large as those that were observed if the obliquities were drawn from an isotropic distribution (i.e., a uniform distribution in cos ψ). A subtlety of these calculations is that for given values of the observables | cos i| (from the transit impact parameter), | sin i_★| (from rotational broadening) and λ (from the RM effect), there are two closely-spaced solutions for ψ: $$\cos \psi =\left|\sin i_{\star }\sin i\right|\cos \lambda \pm \left|\cos i_{\star }\cos i\right|.$$(3)

Following Albrecht et al. (2021), our Monte Carlo simulations take this discrete degeneracy into account by randomly choosing one of the two solutions in each iteration. The results were p = 9.4 × 10^–4 for the dispersion and p = 3.8 × 10^–4 for the standard deviation. Applying the same test to the sample of 14 systems that were available to Albrecht et al. (2021) and were not based on gravity darkening gives p = 1.4 × 10^–2 and 2.1 × 10^–3, respectively. Thus, the inclusion of 13 new data points has reduced the p-values by factors of 15 and 5.5, allowing a firmer rejection of the null hypothesis.

Next, we used the Dong & Foreman-Mackey (DFM) code⁴ provided by Dong & Foreman-Mackey (2023) to fit a model to the obliquity data consisting of the sum of two beta distributions, representing the aligned and misaligned populations.⁵ Each beta distribution is characterized by three hyperparameters: a weight, w, a mean, µ, and an inverse variance, κ. We adopted the same priors on the hyperparameters as outlined in Dong & Foreman-Mackey (2023). the results for the hyperparameters when the DFM code is applied to different samples of the data are available on Zenodo (see Data availability), and Fig. 3 shows the posterior probability density for cos ψ.

The left panel of Fig. 3 shows the results of analyzing the complete sample of projected obliquity measurements. There is a strong peak at cos ψ = 1 representing the aligned systems, and a much broader distribution representing the misaligned systems. The distribution falls off toward retrograde systems (cos ψ < 0.75) but otherwise there is little to no evidence for a concentration near cos ψ = 0, in agreement with the findings by Dong & Foreman-Mackey (2023) and Siegel et al. (2023).

The middle panel of Fig. 3 shows the results of analyzing only those systems for which constraints on i_★ are available, but without using the i_★ information to perform the inference. The inferred probability density of cos ψ does not look too different from the case described above.

The right panel of Fig. 3 shows the results of analyzing only those systems for which constraints on are available and using the i_★ information to perform the inference. This is the type of sample that led Albrecht et al. (2021) to propose that there is a peak near cos ψ = 0 and, indeed, we find a similar result.

The preceding results suggest that the obliquity distribution of the systems for which measurements of both λ and have been reported is different from that of the systems for which λ was measured and no constraints are available on i_★. This is true even though we have excluded inclination measurements obtained via gravity darkening, thereby eliminating the bias discussed by Siegel et al. (2023).

It is probably relevant that the λ & sample is enriched in the types of systems that are known to show frequent misalignments: sub-Saturns (0.025 M_J < M_p < 0.2 M_J), and hot Jupiters around hot stars (a/R_★ < 7, M_p > 0.3 M_J, T_eff ≥ 6250 K). Stellar inclination measurements in those systems almost all come from the v sin i_★ method - the combination of a photometric rotation period, rotational Doppler broadening of the spectral absorption lines, and the stellar radius. In the entire sample, there are 17 sub-Saturns, of which inclination measurements are available for 11, of which 5 have cοs ψ < 0.75. The entire sample also contains 41 hot Jupiters around hot stars, of which 15 have inclination measurements. After excluding the gravity-darkening measurements, we are left with 10 hot Jupiters, all of which have cos ψ < 0.75.

By combining the 5 misaligned hot Neptunes and the 10 misaligned hot Jupiters, we constructed a sample of 15 systems and subjected it to the same Monte Carlo experiments described above to test the null hypothesis that they are drawn from an isotropic distribution. For the dispersion-around-zero and standard deviation, we found p = 3.1 × 10^–3 and p = 3.5 × 10^–3, respectively. We also applied the DFM code to a sample consisting only of sub-Saturns, and a sample consisting only of hot Jupiters around hot stars, with results that are given on Zenodo (see Data availability) and displayed in Fig. 4, in the same format as in Fig. 3. For the sub-Saturns, a central peak appears to be favored by the data, although the uncertainty in the distribution is large due to the small number of systems (17). A central peak is even more pronounced in the sample of hot Jupiters around hot stars (41 systems including the measurements employing the gravity darkening method). For the subset of systems with inclination measurements from this hot Jupiter sample, since all of the systems are misaligned, we modeled the distribution with a single Beta distribution rather than using two components. There is a clear peak at cos ψ ≈ –0.2 or ψ ≈ 100°. Indeed, all of these ten systems are consistent with ψ = 90° with varying degrees of uncertainty, and a simple kernel-density estimate of the distribution shows a peak at cos ψ = 0. Both distributions are shown in Fig. 4.

To investigate the importance of the priors assumed for the hyperparameters in the "DFM" code for these subpopulations, we tried adopting less informative priors on the variance, κ or rather log κ, hyperparameter. Here, we tried a uniform prior instead of the Gaussian prior used in Dong & Foreman-Mackey (2023, i.e., 𝒰(–4, 10) instead of 𝒩(0, 3)). The resulting hyperparameters are available on Zenodo (see Data availability) and distributions are displayed in Fig. 5. As the variance (κ) is now poorly constrained the distributions do look somewhat different, and the peaks for the misaligned populations, especially for the hot Jupiter sample, are not as sharp. At present there are probably too few systems or too little information to properly infer the distributions for these two subpopulations, without assuming some value for the variance to describe the morphology. When applying the information on we have for the subset of ten systems, the resulting distribution looks very similar to that in Fig. 4.

Although the results are too prior-dependent to be sure, these results lead us to hypothesize that the preponderance of polar planets is a phenomenon specific to hot Jupiters around hot stars and possibly also sub-Saturns. The hot stars, in particular, are more amenable to the determination of via the v sin i_★ method, given their higher rotation velocities, possibly explaining their over-representation in the λ & i_★ sample.

If there truly is a preponderance of systems with nearly polar orbits, it would seem to be a clue about the history of these systems. For example, for the sub-Saturn population, Petrovich et al. (2020) proposed that nearly polar orbits of Neptune-mass planets can result from a secular resonance between the nodal precession frequencies induced by the protoplanetary disk and by an outer companion. We might test this theory by searching for distant massive companions to the polar-orbiting planets. In any case, Petrovich et al. (2020) noted that their theory is not applicable to hot Jupiters.

Tidal evolution might lead to a pile-up of polar orbits of hot Jupiters around hot stars. Lai (2012) demonstrated that orbital orientations can linger near 90° when damping is from the dissipation of inertial waves in the convective envelope of a star. In general, this theory requires the planet’s orbital angular momentum to be larger than the star’s spin angular momentum to avoid producing more systems near 180° than are observed (e.g., Li & Winn 2016).

Fig. 2 Obliquities and projected obliquities as a function of Teff (left), a/R* (center), and Mp/M* (right). Triangles indicate a measurement whose value is outside of the plotted range. In all panels, the color of each data point conveys Teff using the color scale shown in the lower left panels. The contours are KDEs illustrating the density of measurements in a given parameter space.

Fig. 3 Hierarchical Bayesian inference of the obliquity distribution using the code of Dong & Foreman-Mackey (2023). The left panel is for all systems for which λ was measured, the middle panel is for systems for which both λ and i★ were measured but without using the i★ information, and the right panel is for systems for which both λ and i★ were measured and making use of the i★ information.

Fig. 4 Hierarchical Bayesian inference of the obliquity distribution for subsets of the sample that are suspected of being especially prone to misalignment. The left panel is for all sub-Saturns with λ measurements, without using any i★ information even when it is available. The middle panel is for all hot Jupiters around hot stars with λ measurements, without using any i★ information even when it is available. The right panel is for hot Jupiters around hot stars for which both λ and i★ are known and utilized. In the right panel, the black upward tick marks on the horizontal axis are the individual values of cos ψ, and the horizontal bands surrounding the tick marks convey the uncertainties. The solid gray curve is a KDE estimated from these measurements.

Fig. 5 Hierarchical Bayesian inference of the obliquity distribution for the same subsets as in Fig. 4, but here a uniform prior, 𝒰(–4, 10), was applied to log κ in the "DFM" code.

4.2 Exploring a possible obliquity-eccentricity correlation

Explaining the existence of giant planets on orbits well within the "ice line" of the protoplanetary disk is one of the oldest unsolved problems in exoplanetary science. Disk-driven migration is one possibility, and would generically yield low eccentricities and obliquities (e.g., Lin et al. 1996; Dawson & Johnson 2018). In this scenario, the observed cases of high eccentricity and obliquity arise from other processes. An alternative idea is to combine eccentricity-raising dynamical interactions - such as planet-planet scattering or Kozai-Lidov cycles - and tidal dissipation (e.g., Fabrycky & Tremaine 2007; Chatterjee et al. 2008; Nagasawa et al. 2008). If the eccentricity reaches a high enough value, or more to the point, if the periastron distance reaches a low enough value, the damping of tidal oscillations in either the planet or star can convert orbital energy into heat and thereby shrink the semi-major axis.

In many of the proposed scenarios for raising of the eccentricity, the orbital inclination relative to the initial plane is also increased, which could lead to a high stellar obliquity. Therefore, if the dynamical scenarios are correct, one might expect to see a correlation between e and λ for the affected planets. An important caveat is that tidal damping tends to reduce both the eccentricity and stellar obliquity, but on timescales that might differ by orders of magnitude. This is not only because eccentricity and obliquity tides excite different types of perturbations, but also because eccentricity damping can take place inside either body while significant obliquity damping requires the energy to be dissipated inside the star (Ogilvie 2014). Thus, tides might confound the interpretation of any correlation. For example, if eccentricity damping is much faster than obliquity damping, then any initial correlation might be diminished on astronomically short timescales. Or, tidal dissipation might increase the correlation by boosting the number of low-e, low-/l systems even if there were little or no correlation before tidal evolution.

Regarding the latter point, a/R_★ and λ are known to be positively correlated for cool stars (Albrecht et al. 2022). This has been interpreted as evidence for tidal dissipation, since dissipation rates should decline sharply with a/R_★. Furthermore, a/R_★ is known to be positively correlated with e, presumably for the same reason (although the causal variable might be a/R_p in that case). Taken together, these two correlations might cause e and λ to be correlated without invoking any new explanation.

Nevertheless, it is well worth testing for a correlation between e and λ in the current sample and Rice et al. (2022) investigated both of these parameters. Specifically they found at the Kraft break, the cumulative sum of λ for systems with eccentric orbits differs by 6.5σ from the λ distribution of systems with circular orbits. For systems with giant planets they obtain an 8.7σ result. Their procedure was as follows, as we understand it. They selected systems with periastron distances ≤ 0.1 AU for which a value of λ was reported in the TEPCat obliquities table, looked up e in the NASA Exoplanet Archive (ΝΕΑ) Confirmed Planets list, and separated the systems into two samples: an eccentric sample (e ≥ 0.1), and a circular sample (e = 0). They omitted borderline systems with 0 < e < 0.1. They arranged the planets in order of the star’s effective temperature, giving each planet a position i in a line of Ν planets. They examined the cumulative sums of |λ| for the eccentric sample, with the jth cumulative sum defined as the sum of all |λ| values for i ≤ j. Since the eccentric sample contained fewer systems and had a different distribution of stellar types than the circular sample, Rice et al. (2022) used a Monte Carlo procedure to make "matching" circular samples for comparison. They created 5000 matching circular samples by drawing randomly (with replacement) from the circular sample the same number of stars with T_eff < 6100 Κ and the same number of stars with T_eff > 6100 Κ as in the eccentric sample. They then compared the Nth cumulative sums of the eccentric sample and the ensemble of matching circular samples, and found them to differ by 6.5σ, where σ was determined from the spread amongst the 5000 trials. For convenience we will denote the "number-of-sigma" by the symbol Δ.

They also tried some variations on the samples. When performing the same steps but with the samples restricted to giant planets, defined as M_p > 0.3 M_J, they found Δ = 8.7. To test the notion that a/R_★ is a "hidden variable" responsible for the e/λ correlation, they created samples that were not divided by eccentricity, but rather on whether a/R_★. is greater or smaller than 12. In this case, they found Δ = 4.9, and since this is smaller than 6.5 they concluded that the variation in a/R_★. cannot explain the entire effect they observed between e and λ.

With the enlarged sample, we decided to revisit these intriguing results. After implementing the method of Rice et al. (2022) to the best of our ability, we found Δ = 2.1 for the whole sample shown in Fig. 6 and Δ = 4.1 for planets with M_p > 0.3 M_J. To check on our code, we tried reproducing their results using the systems known at that time and were unsuccessful. Further investigation showed that the main reason for the discrepancy can be traced to 6 systems in the eccentric sample that should have been in the circular sample. These 6 systems (with details in the footnotes) are HAT-P-23 b and HAT-P-32 b⁶, KELT-6 b⁷; and HATS-14 b, HATS-70 b, & Kepler-63 b⁸. Taken together, these errors caused σ to appear larger than in reality.

We tried some variations on the samples. When including systems with 0 < e < 0.1 in the circular group, instead of omitting them, we found Δ = 1.8. When focusing on giant planets (M_p > 0.3 M_J), we obtained Δ = 4.6. We note that there are only 10 known eccentric giant planets around cool stars, making the results sensitive to an error for an individual system. When we tried dividing the sample according to whether a/R_★ is smaller or greater than 12, instead of by eccentricity, we found Δ = 8.1 as shown in Fig. 7. Thus, in our case, the σ statistic is larger when testing the effect of a/R_★ than when testing for the effect of eccentricity.

Fig. 6 Obliquity distribution for stars hosting planets on eccentric and circular orbits. Left: the cumulative sum of |λ| for eccentric planets shown in red, sorted by Teff. The black line is the average of 5000 randomly sampled sets from the circular population with the gray shaded area denoting the 1σ and 2σ intervals. The vertical blue line denotes the Kraft break (Teff = 6100 K). Right: histogram showing the distribution of the circular population and the value of the eccentric population at the Kraft break. The Gaussian is a fit to the histogram. The circular and eccentric samples only differ by 2.1 σ, meaning the eccentric population does not appear to be significantly more misaligned. Adapted from Rice et al. (2022).

Fig. 7 Obliquity distribution for stars hosting planets on close-in and wide orbits. Same as Fig. 6, but for systems with planets on close-in or wide orbits with the division given for a/R★ = 12. The systems with planets on wider orbits appear to be significantly more misaligned compared to those with planets closer in.

4.3 Tidal realignment

The main evidence that tidal obliquity damping occurs on astro-physically relevant timescales, despite many researchers’ early expectations (see, e.g., Queloz et al. 2000; Winn et al. 2005), is that the types of systems where one might expect tidal dissipation rates to be fastest - with massive, close-orbiting planets and host stars with thick convective envelopes - do indeed have an obliquity distribution that is concentrated around 0°. And conversely, systems for which misalignment is proposed to be part of the formation process but for which tidal dissipation rates are expected to be much slower - somewhat wider-orbiting giant planets, and stars with radiative envelopes - have a broad obliquity distribution. For some details on these findings, see Schlaufman (2010); Winn et al. (2010), and Hébrard et al. (2011). To demonstrate this point more explicitly, Albrecht et al. (2012) constructed a crude "tidal dissipation figure of merit" for each system, meant to be a proxy for the tidal dissipation rate, and showed that the obliquity distribution does indeed narrow for the systems where the figure of merit predicts more rapid dissipation.

However, the precise mechanisms of tidal dissipation are not well understood and are probably quite complex, depending on forcing frequency, stellar structure, rotation rate, etc. Another problem with tidal damping is that the expected endpoint of tidal evolution for most of the known close-orbiting giant planets is tidal orbital decay and destruction (Brown et al. 2011). Thus, an observed planet around a star that experienced tidal damping must somehow have avoided orbital decay before the star was aligned to the observed degree.

Several theories exist that might prevent or at least delay planet destruction for long enough to be consistent with the data. The orbital angular momentum of some lucky planets may exceed the host’s spin angular momentum, causing the star to align before its orbital decay (see, e.g., Hansen 2012; Valsecchi & Rasio 2014; Dawson 2014). If tidal damping is primarily driven by the dissipation of inertial waves in a star’s convective zone, then alignment could be faster than orbital decay (Lai 2012; Lin & Ogilvie 2017; Damiani & Mathis 2018). The damping of gravity waves in radiative cores might similarly cause alignment before tidal in-spiral (Zanazzi et al. 2024). Another idea is that if the efficiency of tidal dissipation drops sharply with increasing forcing frequency (shorter orbital periods), then alignment can be arranged to occur before tidal evolution slows down dramatically, allowing the orbit to survive (Barker & Ogilvie 2010; Penev et al. 2018; Anderson et al. 2021; Ma & Fuller 2021). Alternatively, and more speculatively, tides from the planet might only affect a star’s outer convective zone, allowing the planet to donate a smaller amount of its orbital angular momentum to spin up the star and save itself (Winn et al. 2010). With all of this theoretical uncertainty, we considered whether there could be a relatively model-independent way to seek corroborating evidence for tidal obliquity damping.

The Sun’s obliquity is 7.155° relative to the ecliptic (Beck & Giles 2005) and 6.2° relative to the Solar System’s invariable plane (see, e.g., Souami & Souchay 2012; Gomes et al. 2017). If we assume this angle to be representative of the degree to which stars and their planetary systems are initially aligned - admittedly, a big assumption - then we would expect some systems that have experienced strong tidal damping to have obliquities much smaller than 6°. This can be tested with high-precision RM observations of the types of systems where strong tidal damping is expected and comparisons to systems where it is not expected.

Focusing on only the most precise measurements (σ_λ < 2°) in cool hot Jupiter systems, Albrecht et al. (2022, Section 3.1.8) found the dispersion of λ measurements to be 0.91°. Given an average measurement uncertainty of 0.82°, this should probably be regarded as an upper limit on the true dispersion of a few degrees. There are drawbacks to this simple approach. First, one must decide what constitutes a "precise" measurement. Somewhat arbitrarily, Albrecht et al. (2022) chose a threshold precision of 2°. Second, it ignores the information contained in the much more numerous and less precise measurements.

We tried overcoming these problems with a hierarchical Bayesian model (HBM), following up on the work by Siegel et al. (2023), who carried out similar tests to those described below. (To eliminate any suspense, we obtained similar results, too.) The hyper-parameter of interest is the dispersion of the underlying λ distribution around 0°, independent of measurement uncertainties. To identify systems for which rapid tidal dissipation is expected, we required T_eff < 6100 K, a/R* < 10, and M_p > 0.3 M_J. We also required σ_λ < 50°. This led to a sample of 60 systems, only one of which stood out: TOI-858Bb (Hagelberg et al. 2023), which is located in a wide binary system and has a retrograde orbit. Since tidal dissipation is evidently not the only important obliquity-altering process in this system (perhaps due to the stellar companion) we decided to omit it from consideration, leaving a sample of 59 systems. Fig. 8 shows the key characteristics of those 59 systems (as well as some systems with higher values of a/R_★). The data points with a precision in λ better than 2° are highlighted; the measurement of WASP-50 of –2.9 ± 1.2° presented in this paper is included in this sample.

The result of the HBM was a dispersion in λ of 1.4 ± 0.7°, depicted in Fig. 8 as a Gaussian function with a width of 1.4°. For prograde systems the expected endpoint of the alignment process is λ = 0°, which is why our HBM did not include a hyperparameter for the mean of the obliquity distribution. However, as an experiment to check for surprises or systematic errors (such as the neglect of the convective blueshift) we repeated the analysis after allowing for the mean to be a second hyperparam-eter. The results were a dispersion of 1 .5 ± 0.8° and a mean of –0.2 ± 0.7°, giving no clear evidence for surprises or systematic errors.

In the light of these results, our measurement of WASP-50 (–2.9 ± 1.2°) is more interesting than we originally anticipated. According to the rough criteria described above, the tidal dissipation timescale for this system should be faster than in most of the other systems in the sample. Gillon et al. (2011) reported an age of 7.0 ± 3.5 Gyr based on comparing the usual photometric and spectroscopic observables to stellar-evolutionary models. However, the estimated age based on the chromospheric activity and rotational period was only 0.8 ± 0.4 Gyr (which could also be a result of tidal spin-up, e.g., Tejada Arevalo et al. 2021). If the old age is correct, the system is more likely to have had enough time to achieve good alignment than if the young age is correct. As discussed in Section 3.2.2, Gillon et al. (2011) reported a rotation period for WASP-50 of 16.3 ± 0.5 d or potentially twice this value. Using these values ψ comes out to ${34}_{-9}^{+12°}$ and ${5.0}_{-3.0}^{+1.9\circ }$, respectively. The value of the rotation period has dramatic consequences for the inferred orientation of the system; the shorter period we identified in the TESS light curve implies an even more misaligned system, which seems worth checking on, given that so many other similar systems are well-aligned. Future TESS observations might be able to resolve the issue.

Considering planets on wider orbits, with 10 ≤ a/R_★ < 20, the HBM yielded a dispersion of 8 ± 2°. The spread in λ is therefore seen to increase with orbital separation (as also discussed in Section 4.2). When restricted to planets with 20 ≤ a/R_★ < 30 (a sample of only 7 systems), the HBM yielded a dispersion of 4 ± 3°, which is consistent with the 10 ≤ a/R_★ < 20 results. It would be interesting to obtain more measurements for planets with large separations to bolster these results and possibly even investigate the separation-dependence of tidal dissipation rates.

Fig. 8 Precise projected obliquity measurements. Here we highlight the most precise (σλ < 2°) measurements of λ for systems with cool stellar hosts (Teff < 6100 K). WASP-50 is highlighted with a black center. The color-coding is done according to Teff and the marker sizes scale with the tidal alignment timescale (as given by Zahn 1977). The gray shaded area shows the Gaussian distribution inferred from cool hosts (Teff < 6100 K) with close-in (a/R★ < 10), massive (Mp < 0.3 MJ) planets. The horizontal lines demarcate the Solar System obliquity of 7.155° (Beck & Giles 2005). The light gray error bars show the less precise measurements. Adapted from Albrecht et al. (2022, Figure 10).

4.4 Compact multi-transiting systems, with and without outer companions

At a time when only five λ measurements had been made for stars with more than one known transiting planets, Albrecht et al. (2013) noted that they were all consistent with good alignment. The number is now 22, including four systems we have added to the tally in this work. For WASP-148, λ had already been measured by Wang et al. (2022), and our result is consistent with theirs but more precise. For TOI-1130, we found $\lambda =4_{-6}^{+5\circ }$, while for two other systems, LTT 1445A and TOI-451A, our formal results were that the orbits are prograde, but given the low S/N of our measurements we do not include them in the HBM presented at the tail end of this section. We do include an additional 4 systems for which there are good constraints on the stellar inclination even though λ is not known. Our discussion is therefore based on the properties of 24 systems.

Figure 9 is a family portrait of these systems. The Solar System is shown at the top for scale and reference. There are 19 systems for which the data are consistent with zero obliquity; there are 5 systems with large obliquities relative to at least one transiting planet: HD 3167, Kepler-56, HIP 41378, Kepler-129, and K2-290. HD 3167 has a very unusual architecture in which the planets appear to be on perpendicular orbits (Bourrier et al. 2021). This case is illustrated in Fig. 9 as the outer planet floating over Kepler-129. HIP 41378 is a five-planet system in which Grouffal et al. (2022) measured the projected obliquity with respect to planet d to be ${57}_{-18}^{+26}°$. However, due to the long duration (~12.5 hr and P~278 d), a full transit was not observed, preventing us from having high confidence in the result.

In Fig. 9, Kepler-30 is grouped with the aligned systems, but there are conflicting indications about the obliquity. By using two different methods Sanchis-Ojeda et al. (2012) found λ = −1 ± 10° and λ = 4 ± 10°, indicating a well-aligned system. However, Morgan et al. (2023) found $i_{\star }={43}_{-9}^{+15}°$, indicating a large misalignment, based on the combination of the rotation period, v sin i_⋆, and R_⋆. However, this result depends critically on the value v sin i_⋆ = 2.0 ± 0.2 km s⁻¹ reported by Fabrycky et al. (2012), and the measurement uncertainty of 0.2 km s⁻¹ is likely to have been underestimated. Systematic errors of at least 0.5 km s⁻¹ are typical, due to effects such as turbulent and instrumental broadening. A value of 3.0 ± 0.5 km s⁻¹ would allow for consistency with alignment (i_⋆ ≈ 90°).

In contrast to hot Jupiter systems, where misalignments might have occurred in the tumultuous process of high-eccentricity migration after disk dispersal, the misalignments in multiplanet systems might have occurred while the disk was present. Gravitational interactions between a planet-forming disk, a binary companion, and a spinning host star can be enhanced when the disk looses mass, magnifying the host star’s tilt through a secular resonance (e.g. Batygin & Adams 2013). For a binary companion to generate spin-orbit misalignments, the stellar spin must be gravitationally decoupled from the planet forming in the disk, otherwise the star’s spin axis and planet’s orbit normal precess in-unison (Eq. (80) of Zanazzi & Lai 2018). The spin is decoupled from the orbits of LTT 1445Ab, and K2-290Ab,c, but not for TOI-451Ab, when torqued by their binary companions.

Similarly, sweeping secular resonances from a dispersing disk can enhance gravitational interactions between forming planets, magnifying their inclinations (e.g. Ward 1981). However, for planetary companions to generate misalignments, they must also supply the system with a large angular momentum deficit to misalign an inner planet. Based on the work by Petrovich et al. (2020); Zanazzi & Chiang (2024), we assume a planetary companion (mass m₂, semi-major axis a₂) could secularly misalign an inner planet (mass m₁, semi-major axis a₁) if they are well separated (a₂/a₁ > 2), and they have a large angular momenta ratio $(\left(m_2{a}_2^{1/2}\right)/\left(m_1{a}_1^{1/2}\right)>10$, see e.g. Eq. (22) of Zanazzi & Chiang 2024). The planet pairs Kepler-56 c,d and Kepler-129 c,d meet this criteria, while Kepler-25 c,d, WASP 47 b,c, and TOI-1130 c,d do not, although we note that the companion TOI-1130 d is as of yet poorly constrained (Korth et al. 2023). LTT 144A might be aligned when a misalignment is predicted by Zanazzi & Lai (2018), while the other systems with distant, massive companions appear consistent with spin-orbit misalignments being generated (or not) while the gas disk is present (Fig. 9).

Secular interactions from companions exciting an inner planet’s inclination during the disk-hosting phase is not the only mechanism which can explain misalignments in multiplanet systems. For instance, Gratia & Fabrycky (2017) postulate a dynamical instability between three outer planets could have misaligned Kepler-56 b and c, while Best & Petrovich (2022) find secular chaos induced by the two binary companions of K2-290A is sufficient to misalign K2-290Ab and c, with both mechanisms taking place after disk dispersal. Additional stellar obliquity measurements in multiplanet systems with companions (e.g. LTT 1445A), and further constraints on the presence or absence of companions in inclined multiplanet systems (e.g. HD 3167, HIP 41378), could constrain which mechanisms are responsible for generating misalignments.

For the multiplanet systems which have obliquities consistent with alignment (again barring TOI-451A and LTT 1445A), we applied the HBM framework described in the previous section to infer the mean µ and dispersion σ of the intrinsic distribution of obliquities. The results were µ = 0.0 ± 1.7° and σ = 2.0 ± 1.5°. When we removed Kepler-30 c from the sample, out of the concern raised above, we found µ = 0.0 ± 1.8° and σ = 2.1 ± 1.6°. Thus, the obliquity dispersion of the host stars of multis is similar to that of cool stars with hot Jupiters. However, tidal obliquity damping is expected to be much weaker for the multis than for the hot Jupiters, because of the planets’ lower masses and relatively wider orbits. The low obliquities of the multis might therefore be primordial, namely, due to the flatness of the protoplanetary disk and its strong coupling to the accreting young star. From this point of view, the solar obliquity of 6° may be unusually high – and if the high solar obliquity is simply a peculiarity of the Solar System, then the evidence for tidal obliquity damping presented in the previous section would be undermined. Finally, we note that an additional aligned system could potentially be the HD 148193 system ($\lambda ={14}_{-8}^{+7}°$, Appendix A.2), if the inner transiting candidate is confirmed as a bona fide planet at some point.

Fig. 9 Obliquities of stars with multiple transiting planets, as constrained by either the RM effect or asteroseismology. Each horizontal line shows the a/R⋆ values of the planets in a given system, with the solar system on top (with gaps to allow Jupiter and Neptune to be shown). For the planets, the size of each circle conveys the planet’s radius relative to Jupiter. For the stars, symbol size conveys the radius relative to the Sun, with a color conveying the effective temperature. Gray circles are planets for which the RM effect has been measured – except for β Pic b where spectro-interferometry was used. The wedge indicates the projected obliquity for a given system, based on the most precise measurement for any planet in the system, with the exception of HD 3167 for which two planets appear to have very different inclinations. The asteroseismic measurements constrain i⋆ and not λ; in those cases the wedge should therefore be read as aligned (Kepler-50 and Kepler-65) or misaligned (Kepler-56, Kepler-129). Planetary (tan circle) or stellar companions (red star, M-dwarfs) to a given system are shown to the right, where √(÷) denotes that the companion could (not) have influenced the obliquity of the inner planetary system. The asterisks for LTT 1445A and TOI-451A denote that these measurements should be taken with a grain of salt. Adapted from Wang et al. (2022, Figure 3).

5 Conclusions

We presented new observations of the RM effect for 19 stars with transiting planets. We combined these measurements with results from the literature to arrive at a sample of 205 systems for which λ has been determined with reasonable confidence. With this sample, we revisited our four primary questions, as detailed below:

• (i)

  The possible preponderance of perpendicular planets was first raised as a possibility by Albrecht et al. (2021) and then questioned by Siegel et al. (2023) and Dong & Foreman-Mackey (2023). Our updated sample does not settle the matter. It remains the case that a “polar peak” is seen in the obliquity distribution of the sample of planets for which λ and i_⋆ have been measured, but no such peak is seen clearly when analyzing the larger sample of systems for which λ has been measured but not necessarily i_⋆. We showed that the two samples might differ from each other because certain kinds of systems (i.e., hot stars with hot Jupiters) are more amenable to i_⋆ measurements than others and also more likely to have near-polar orbits. The subsample of sub-Saturns around cooler stars also seems to have a preference for polar orbits, although the sample needs to be increased to be sure. Thus, the resolution of this issue still seems likely to be astrophysically interesting;

• (ii)

  With respect to the plausibility of a correlation between orbital eccentricity and stellar obliquity for close-orbiting giant planets, we concluded that the correlation is weaker than previously reported by Rice et al. (2022). Furthermore, we argued that the correlation exists because both eccentricity and obliquity are correlated with a third variable, orbital separation. The causal mechanism is likely to be tidal dissipation, which lowers both eccentricity and obliquity at short orbital separations;

• (iii)

  We considered the obliquities of stars suspected of having been subject to tidal obliquity damping and how low they might be. Measurements of λ of cool stars with close-orbiting giant planets (T_eff < 6100 K, a/R_⋆ < 10, M_p > 0.3 M_J) have an intrinsic dispersion of 1.4 ± 0.7° in λ, according to our hierarchical Bayesian model (HBM). Such fine alignment can be considered as supporting evidence for tidal obliquity damping, if the solar obliquity of 6° is taken as representative of a typical un-damped “primordial” obliquity. More precise measurements of λ in systems hosting massive planets around cool stars might even help to quantify the dependence of tidal dissipation rates on orbital separation;

• (iv)

  Finally, we considered how well aligned the central stars of compact systems of multiple transiting planets are. Except for a few notable exceptions, the λ measurements in compact multiplanet systems are consistent with low obliquities. For this sample, we found an intrinsic dispersion in λ of 2.1 ± 1.6°. In this sense, the multis are better aligned than the Solar System, indicating that compact multis are born unusually flat or that the Sun is unusually tilted, complicating the conclusions drawn above about tidal obliquity damping. The few examples of misaligned stars with compact systems of multiple transiting planets might be explained by interactions between the protoplanetary disk and a stellar companion or by interactions with outer giant planets after the epoch of planet formation.

While the more than 200 measurements of the projected obliquity allow us to investigate different trends, it is clear that the sample as a whole is very diverse. Conclusions based on the whole sample might therefore not apply to all types of systems. As the community continues to amass measurements of the (projected) obliquity in more systems, it will be possible to make more robust inferences about the underlying distributions and trends when considering the more homogeneous subpopulations.

Data availability

Phase-folded light curves for all systems analyzed in Appendix A are available on Zenodo (https://zenodo.org/records/13847422), in particular figures similar to Fig. A.8 for KELT-4A and XO-7 are available. Additional tables containing MCMC posteriors, all RVs presented, and extensions to Tables A1 and A2 of Albrecht et al. (2022) can also be found on Zenodo (https://zenodo.org/records/13897109).

Radial velocities are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/690/A379

Acknowledgements

The authors thank the anonymous referee for helpful and insightful comments that greatly improved the quality of the manuscript. The authors would like to thank R. I. Dawson, J. Siegel, J. Dong, D. Foreman-Mackey, M. Rice, and S. Wang for helpful comments and discussions. Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant agreement no.: DNRF106). E.K. and S.H.A. acknowledge the support from the Danish Council for Independent Research through grant No.2032-00230B. Parts of the numerical results presented in this work were obtained at the Centre for Scientific Computing, Aarhus, https://phys.au.dk/forskning/faciliteter/cscaa/. This research has made use of data obtained from or tools provided by the portal https://exoplanet.eu/home/ of The Extrasolar Planets Encyclopaedia. This research has made use of NASA’s Astrophysics Data System. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exo-planet Exploration Program. This work is based on observations made with the Very Large Telescope (VLT) at Cerro Paranal, Chile, the Telescopio Nazionale Galileo (TNG) at Roque de los Muchachos, La Palma, Spain, and the Nordic Optical Telescope (NOT) also at Roque de los Muchachos. NOT is owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofísica de Canarias. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. The K2 light curves were extracted through the EVEREST pipeline (Luger et al. 2016, 2018) hosted on the Mikulski Archive for Space Telescopes (MAST). The TESS light curves were extracted using Lightkurve (Lightkurve Collaboration 2018). EsoReflex (Freudling et al. 2013) was used to extract CCFs from the ESPRESSO data. This work made use of the following Python packages; NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), corner (Foreman-Mackey 2016), PyMC (Wiecki et al. 2023), and Astropy (Astropy Collaboration 2013, 2018, 2022).

Appendix A Analysis of the individual systems

Below, we describe the particulars of each system, our analysis, and the results. The details for each system are also summarized in Tables 1, 3, and 4.

A.1 HD 118203

HD 118203 b was discovered by da Silva et al. (2006) through RV monitoring, and was later found to be a transiting planet by Pepper et al. (2020) using TESS data. TESS data are available for HD 118203 from Sectors 15, 16, 22, and 49 with a 2 minute cadence. The planet is a hot Jupiter on an eccentric orbit around a slightly evolved G-type star (Pepper et al. 2020). The TESS data centered on the mid-transit time are available on Zenodo (see Data availability).

We observed two transits of HD 118203 b with FIES. The first transit began on March 24, 2020 at around UT 21:00 and ended at 04:15. The observations were interrupted twice due to high humidity. The exposure time was 900 s and the time between exposures was 1171 s, with the time in between spent on ThAr calibration. Two exposures on this night were heavily affected by weather and the resulting RVs deviate significantly from the rest (see the gray points in Fig. A.1). These two data points were omitted from the analysis.

Since the first transit had been interrupted by the weather, we organized observations of a second transit on the night starting on June 18, 2022. Unfortunately, the first part of the transit was missed because the telescope was being used to observe a target-of-opportunity. Our transit observations started at around UT 22:30 and continued until around 03:00. As the FIES spectrograph was refurbished in the summer of 2021 resulting in an improved throughput, we opted for a shorter exposure time of 600 s for the second transit, which resulted in a sampling time of around 768 s. We note that the refurbishment introduced an offset in RV between observations. As such we allowed the systemic velocity parameter to differ between the two nights: γ_FIES applies to the observations taken before summer 2021, and γ_FIES+ applies to subsequent observations.

As mentioned in Section 3 we performed two MCMC analyses, one with a uniform prior on v sin i_⋆ and one with a Gaussian prior. For this system, the results differed substantially. The uniform prior in v sin i_⋆ yielded $\lambda =-{74}_{-9}^{+6}°$ and the Gaussian prior yielded $-{23}_{-38}^{+25}°$. The data obtained from the first transit observations are obviously very noisy. This paired with the second transit observations only covering the second half of the transit meant that the amplitude of the signal was allowed to vary significantly, which made v sin i_⋆ wander off to large values when applying applying a uniform prior. In this case v sin i_⋆ ended up being ${26}_{-23}^{+13}$ km s⁻¹, which is both poorly constrained and at face value in disagreement with the value we derived and the value found in the literature. Therefore, our best estimate is $-{23}_{-38}^{+25}°$ based on the analysis in which v sin i_⋆ was subjected to a Gaussian prior.

A.2 HD 148193

HD 148193 is also known as TOI-1836. The planet HD 148193 b (TOI-1836.01) was confirmed and characterized by Chontos et al. (2024). It is a warm Saturn around a subgiant star in a system with a companion star. A second transiting planet candidate was announced later (TOI-1836.02), which is potentially a super-Earth-sized planet on a 1.7-day orbit. However, the signal-to-noise ratio of the observed transits is modest and the planet remains unconfirmed.

TESS data are available with 30-min cadence from Sectors 16 and 25; with 2-min cadence in Sectors 23, 24, 49, 50, 51, and 52; and with 20-sec cadence in Sector 56. The light curve is available on Zenodo (see Data availability).

We used HARPS-N at the TNG to observe a transit of HD 148193 b. We started observations around UT 21:30 and continued until UT 05:45 on the night starting on June 15, 2021. The exposure time was 1140 s. After accounting for overheads, the sampling time was 1168 s. The transit observations are shown in Fig. A.2.

The combined fit to the light curve and RVs gave $\lambda ={29}_{-13}^{+12}°$. However, we obtained only one pre-ingress RV and two post-egress RVs. With such sparse out-of-transit data, any mismatch between the out-of-transit RV slope observed on that night and the calculated slope based on the planet’s orbital parameters can lead to systematic errors. This problem, combined with the reasons discussed in Section 3.3 and the large projected rotation speed of this system (~ 6 km s⁻¹), led us to analyze the line distortions and Doppler shadow. We found $\lambda ={14}_{-8}^{+7}°$, with results shown in Fig. A.3. We adopted the value from the Doppler-shadow analysis as our best estimate for λ.

A.3 K2-261

K2-261 was observed in Campaign 14 of the K2 mission with 30-minute cadence. The planet was confirmed and characterized independently by Johnson et al. (2018), who named the planet K2-261 b, and Brahm et al. (2019), who named it K2-161 b. Here, we adopt the name K2-261 b, which is also the name chosen by Ikwut-Ukwa et al. (2020). The planet is a warm Saturn-sized planet on an eccentric orbit around a slightly evolved G-type star (Brahm et al. 2019). K2-261 was observed by TESS in Sectors 9, 35, 45, and 46. TESS data are available with 2-min cadence in Sectors 9,45, and 46, and with 20-sec cadence in Sector 35. The light curve is available on Zenodo (see Data availability).

For K2-261 b we used HARPS-N at the TNG to observe a transit. The data are shown in Fig. A.4. We carried out the observations on the night starting on March 24, 2022, from around UT 20:15 until 03:30. We used an exposure time of 600 s, resulting in a sampling time of around 625 s.

We found a value for v sin i_⋆ for K2-261 that is somewhat smaller than the value of 3.70 ± 0.17 km s⁻¹ reported by Johnson et al. (2018), both from our BF analysis (1.6 km s⁻¹) and from our RM observations. The amplitude of the signal is not as strong as would be expected if v sin i_⋆ were really 3.70 ± 0.17 km s⁻¹. A uniform prior on v sin i_⋆ resulted in v sin $i_{\star }={1.0}_{-0.6}^{+0.5}$ km s⁻¹ and $\lambda ={29}_{-56}^{+35}°$, and applying a Gaussian prior yielded a value for the projected obliquity of ${32}_{-33}^{+48}°$. We adopt the latter as our best estimate for λ.

A.4 K2-287

K2-287 was observed in Campaign 15 of the K2 mission. Jordán et al. (2019) discovered and characterized the system, finding the host to be a G-dwarf and the planet to be a warm Saturn-sized planet on an eccentric orbit.

We observed a transit of K2-287 b using ESPRESSO. Observations were carried out from UT 02:06 to 06:16 on the night July 6, 2021. The exposure time was 540 s, which including overhead comes out to a sampling time of roughly 600 s. The RVs derived from the spectra are shown in the left half of in Fig. A.5.

Fig. A.1 RM effect observed in HD 118203. Left: Time series and residuals of the FIES+ RVs (orange) and FIES RVs (green). The black curve is the best-fit model and the gray shading shows the 1 and 2σ confidence intervals created by randomly drawing values from the MCMC chains for b, v sin i⋆, Rp/R⋆, and λ. The dashed curve is the best-fit model assuming λ = 0°. Gray circles are the data that were obtained during poor atmospheric conditions and were excluded from the analysis. Right: Correlations between λ, v sin i⋆, and b. The red lines illustrate the prior that was applied (either uniform or Gaussian).

Fig. A.2 RV-RM effect for HD 148193. Same as for Fig. A.1, but with RVs from HARPS-N.

Given the 30 min time averaging of the K2 photometric time series, and the rather long orbital period of the planet, the K2 light curve is relatively poorly sampled. Furthermore, this star has evaded TESS observations to this point. As of the time of writing, there are no scheduled TESS observations that would cover this star.

This situation, along with the lack of post-egress spectroscopoic observations, prompted us to search for a suitable archival ground-based light curve covering the egress. The one we found was obtained at the Adams Observatory at Austin College on the night starting June 30, 2018⁹. More recently, the system was observed with the CHaracterising ExOPlanet Satellite (CHEOPS; Benz et al. 2021). Borsato et al. (2021) observed three separate transits with a 1-min cadence, and the data from all three provide good coverage of the entiure transit. In our analysis, we included the de-trended light curves presented by Borsato et al. (2021). The K2, CHEOPS, and ground-based light curves are available on Zenodo (see Data availability).

The two RV-RM analyses yielded similar results for λ, with ${21}_{-11}^{+14}°$ for the case of a uniform prior on v sin i_⋆ and ${23}_{-13}^{+12}°$ for the case of a Gaussian prior on v sin i_⋆. We also investigated the in-transit subplanetary velocities, which are shown in Fig. A.6 along with the best-fitting slope. In addition to the errors obtained when extracting the location of the subplanetary velocities, we add two different jitter terms in quadrature; one for the epochs taken when the planet is completely within the stellar disk (dark gray area in Fig. A.6), and another when the planet is at the limb (light gray). This is because there generally (see, e.g., WASP-50 in Section A.12) seems to be some additional scatter in the subplanetary velocities at the limb, which the errors from fitting a Gaussian to the residual peak do not fully capture.

Fig. A.3 Doppler shadow for HD 148193. Left: The observed Doppler shadow created by subtracting the in-transit CCFs from an average out-of-transit CCF. Middle: Model Doppler shadow. Right: Residuals from subtracting the model from the observed shadow. The summed residuals are given in the adjacent panels. The color bar demarcate the contrast. Green horizontal lines denote points of first (solid), second (dashed), third (dashed), and fourth (solid) contact. The solid, orange, vertical lines denote ±v sin i⋆.

Fig. A.4 RV-RM effect for K2-261. Same as for Fig. A.1, but with RVs from HARPS-N.

From this we found values of v sin i_⋆ = 1.16 ± 0.19 km s⁻¹ and $\lambda ={37}_{-9}^{+8}°$, which are consistent with the values derived from the RV-RM effect. Somewhat arbitrarily, we adopt $\lambda ={23}_{-13}^{+12}°$ as our best estimate; this derives from the RV-RM fit with a Gaussian prior on v sin i_⋆. The posterior distributions are shown in Fig. A.5.

A.5 KELTS-3

Pepper et al. (2013) reported the discovery of the KELT-3 system, which consists of a hot Jupiter orbiting a late F star. TESS data with 2-min cadence are available from Sectors 21 and 48. The light curve is available on Zenodo (see Data availability).

We used FΓES to observe a transit on the night starting April 8, 2022, with observations starting around UT 20:50 and continuing until UT 01:30 April 9, 2022. The exposure time was 900 s, which including overhead resulted in a cadence of 1080 s. The RM-RV effect is shown on the left side of Fig. A.7.

As was the case for K2-261, the value for KELT-3’s projected rotation velocity drawn from the literature (10.2 ± 0.5 km s⁻¹; Pepper et al. 2013) is higher than what we obtained from the BF method (7.6 km s⁻¹). The observed amplitude of the RM effect suggests that the true v sin is even lower. Our RV-RM analysis gave 6.4 ± 0.4 km s⁻¹ when using a uniform prior in v sin i_⋆, and λ = −5 ± 4°. When we applied a Gaussian prior on v sin i_⋆, we obtained −6 ± 5°.

Given the discrepancy between the different measurements of v sin i_⋆, we decided to look at the stacked CCFs resulting from our FIES spectra. In Fig. A.8 we show the contours of the peak values obtained from stacking the CCFS on the 2D grid in the space of v sin i_⋆ and λ, compared to our posterior from the RV-RM fit. From this analysis we found v sin i_⋆ = 5.96 ± 0.03 km s⁻¹ and λ = 1.06 ± 0.16°, implying that the system is well-aligned, and that v sin i_⋆ is closer to 6 km s⁻¹ than the value of 10.2 ± 0.5 km s⁻¹ reported by Pepper et al. (2013). Although evidently the uncertainties in some or all of the v sin i_⋆ determinations have been underestimated, we believe that the stacked-CCF method gives the most reliable result for v sin i_⋆.

Fig. A.5 RV-RM effect for K2-287. Same as for Fig. A.l, but with RVs from ESPRESSO.

Fig. A.6 Subplanetary velocities for K2-287. Subplanetary velocities for K2-287 with color-coding as described in Fig. 1. The dark shaded area corresponds to the time the planet is completely within the stellar disk and the light shaded area denote ingress and egress. The error-bars include a jitter term, where two separate terms have been added in quadrature; one for the epochs for the full transit and one for the epochs when the planet is partially outside of the limb of the star. The horizontal black lines in the top panel mark v sin i⋆ from the fit.

This lends further credence to the RV-RM fit with a uniform prior on v sin i_⋆, and as our observations of KELT-3 have a high S/N and cover the transit very well. We therefore adopt the result from this RV-RM fit, λ = −5 ± 4°, as our best estimate for this svstem.

A.6 KELT-4A

KELT-4 is a hierarchical triple star system discovered by Eastman et al. (2016). The primary star, KELT-4A, is an F star orbited by a hot Jupiter. The close binary system KELT-4BC is engaged in a wide orbit with KELT-4A. TESS data are available from Sector 48 with 2-min cadence. The light curve is available on Zenodo (see Data availability).the TESS cameras cannot resolve KELT-4A from KELT-4BC, we included a dilution factor in our analysis of δM = 3.5 ± 0.1 (estimated from Table 3 in Eastman et al. 2016) to correct for the light from KELT-4BC.

On the night starting January 27, 2021, we observed a transit of KELT-4 Ab, starting the observations UT 01:40 and ending them UT 06:45 January 28, 2021. The exposure time was 900 s, meaning the effective sampling was a spectrum every 1090 s.

The observations are shown in Fig. A.9 with our results from a run invoking a Gaussian prior on v sin i_⋆. From this we got $\lambda ={80}_{-22}^{+25}°$. Applying a uniform prior on v sin i_⋆ we found $\lambda ={77}_{-23}^{+25}°$. Evidently, applying a Gaussian or a uniform prior on v sin i_⋆ makes little difference, despite the somewhat noisy data, most likely because of the rather large impact parameter and large misalignment.

As for KELT-3, we investigated the stacked CCFs for this system. The results coincided with the result from our RV-RM fit, and when fitting a 2D Gaussian to the grid of v sin i_⋆ and λ we obtained values of v sin i_⋆ = 6.0 ± 0.4 km s⁻¹ and λ = 88.9 ± 0.5°. The results are consistent with the polar configuration implied by the RV-RM fit. Our best estimate for λ for this system, ${80}_{-22}^{+25}°$ is based on the RV-RM fit with a Gaussian prior on v sin i_⋆.

A.7 LTT 1445A

LTT 1445 is a triple star system consisting of three M-dwarfs. The primary star hosts an Earth-sized planet on a 5-day orbit, LTT l445Ab, which was discovered by Winters et al. (2019). A second transiting planet, LTT 1445 Ac, was reported by Winters et al. (2022), and an RV signal from a third planet was detected by Lavie et al. (2023). The orbit of LTT l445Ab is sandwiched between the orbits of the other two planets.

LTT 1445A was observed by TESS in Sectors 4 and 31, and the data are available with 20-sec and 2-min cadence, respectively. The light curves show a photometric quasiperiodicity that probably arises from rotation of either star В or star С (Winters et al. 2019). To account for this rotational modulation, our model for the light curve including a GP with a kernel consisting of two stochastically driven damped harmonic oscillators. As the amplitude of the rotational modulation seemed to change between the two sectors (Winters et al. 2022), we did not require the hyper-parameters to be the same in both sectors. The light curve is available on Zenodo (see Data availability).

Fig. A.7 KELT-3. Same as for Fig. A.1. Here also with RVs from FIES.

Fig. A.8 Peak of the stacked CCFs of KELT-3. The peak of the stacked CCFs of KELT-3 created by assuming the best-fitting value of b from the RV-RM fit. The gray contours show a peak around (v sin i⋆,λ)=(6.0 km s−1, 1.0°). The red contours are the same as created from the posterior shown to the right in Fig. A.7.

We observed a spectroscopic transit of LTT 1445 Ab employing HARPS-N. The observations were carried out on the night starting September 5th, 2020. We used an exposure time of 900 s, yielding a sampling time of 920 s. The first exposure was not useful because the telescope was erroneously aimed at LTT 1445BC. The remainder of the observations are shown in Fig. A.10.

As can be seen in Fig. A.10 no RM effect is apparent in the RV data. Formally, we found $\lambda ={22}_{-83}^{+98}°$, implying a gentle preference for a prograde orbit. However, as described by Albrecht et al. (2011), such low S/N results are probably even more uncertain than the formal uncertainty suggests.

Usually, if the RV-RM effect has a detectable amplitude, then the scatter of the RVs inside the transit should be larger than the scatter of the RVs outside the transit. For LTT 1445A, the in-transit scatter is σ = 5.9 m s⁻¹ after subtracting the best-fit orbital solution. The out-of-transit scatter is σ = 7.1 m s⁻¹ (dropping to σ = 6.0 m s⁻¹ if the single pre-transit datum is excluded). The in-transit residuals relative to the best-fit RM model have a scatter of σ = 5.5 m s⁻¹. In light of these results, we do not consider the RM effect to have been detected at all. Better results would probably require a bigger telescope.

A.8 TOI-451A

TOI-451 is a young multiplanet system in the Pisces-Eridanus stream consisting of the G-dwarf planet host, TOI-451A, and a wide-orbiting (likely) M-dwarf binary (Newton et al. 2021). The planetary system around TOI-451A consists of three planets with radii 1.9, 3.1, and 4.1 R_⊕ and periods of 1.9, 9.2, and 16 days. TESS observed the system in Sector 4, 5, and 31 with a cadence of 2 min. The TESS observations of TOI-451A are available on Zenodo (see Data availability).

We observed a transit of TOI-45lAb, the innermost and smallest planet, using ESPRESSO. The observations were carried out on the night starting on September 24, 2022, from around UT 04:15 until 08:30 September 25, 2022. The exposure time was set to 720 s, resulting in a spectrum every 755 s.

As reported by Newton et al. (2021), the light curve of TOI-451 displays a rather large amplitude rotational modulation. We therefore adopted the same kernel as for LTT 1445 (described in the previous section) instead of the Matérn-3/2 kernel. We included all three planets in our transit modeling. We applied Gaussian priors with values drawn from Newton et al. (2021) for the transit parameters to the two, outer planets (c and d) for which we did not observe a spectroscopic transit. For planet b we followed our usual approach.

The RV-RM analysis with a uniform prior on v sin i_⋆ gave $\lambda =-{27}_{-74}^{+64}°$ and the analysis with a Gaussian prior gave $-{23}_{-40}^{+37}°$. All we can conclude is that the orbit is prograde.

We investigated the RV scatter in the same manner we did for LTT 1445. We obtained σ = 2.8 m s⁻¹ and σ = 2.4 m s⁻¹ for the in-transit and out-of-transit scatter, respectively, and when subtracting the RM effect from the in-transit data we found σ = 2.0 m s⁻¹. With these results we cannot be confident that the RM effect was detected. For better results, one might observe a second transit of TOI-45lAb or give up and choose to observe one of the larger planets (for which opportunities to observe transits are rarer).

Fig. A.9 KELT-4A. Same as for Fig. A.1. Here also with RVs from FIES.

Fig. A.10 LTT 1445A. Same as for Fig. A.1, but with RVs from HARPS-N.

We did not use any results from LTT 1445A or TOI-451A in the ensemble analyses described in Section 4.2, Section 4.1, and Section 4.3. However, the results are displayed in Section 4.4, which focuses on multiplanet systems.

A.9 TOI-813

Eisner et al. (2020) used TESS data to discover a Saturn-sized planet on an 84-day orbit around this subgiant star. TESS data are available with 30-min cadence in Sector 2; 2-min cadence in Sectors 5, 8, 11, 30, 33, and 39; and 20-sec cadence in Sectors 27 and 61. Due to the planet’s long orbital period, transits only occurred during Sectors 2,27, and 61. The light curve is available on Zenodo (see Data availability).

We observed a transit of TOI-813 b on the night January 20, 2023, using ESPRESSO. The exposure time was set to 800 s, yielding a sampling time of around 832 s. This transit was simultaneously observed with TESS. We observed the transit starting just before ingress and continued for for 6.5 hours until the target set. These observations cover roughly half of the 13.14 hr transit as shown in Fig. A.12.

In an attempt to cover as much of the transit as possible, the two last exposures were acquired through a large airmass (> 2.2), and the RVs based on those exposures deviate from the trend established by the preceding RVs. Including the two high-airmass RVs in our fit resulted in $\lambda =-{60}_{-21}^{+19}°$ suggesting that the system is significantly misaligned. Omitting these two exposures yielded −32 ± 23°. In both cases described above, a Gaussian prior on v sin i_⋆ was applied. Excluding the last two exposures while running an MCMC with a uniform prior in v sin i_⋆, we obtained $-{53}_{-25}^{+22}°\text{ for }\lambda$.

We also analyzed the Doppler shadow for this system and found a value of $-{53}_{-28}^{+25}°$, which is more in-line with the results obtained using a uniform prior on v sin i_⋆ or including the two data points at high airmass. However, the lack of observations outside of the transit left us without a spectral CCF that is free of the RM effect.

Fig. A.11 RV-RM effect for TOI-451. Same as for Fig. A.1, but with RVs from ESPRESSO.

Fig. A.12 RV-RM effect for TOI-813. Same as for Fig. A.1, but with RVs from ESPRESSO.

Since only half of the transit was covered, we found it necessary to apply a Gaussian prior on v sin i_⋆ to obtain reliable results. We also decided not to trust the two high-airmass RVs. We adopt the resulting value of −32 ± 23° as our best estimate of λ for TOI-813. The RM effect has a high enough amplitude to be detectable using smaller telescopes and less precise spectrographs. Another observation covering the second half of the transit is warranted, and would allow us to judge more definitively if the system is aligned or not.

A.10 TOI-892

TOI-892 b is a warm Jupiter on a circular orbit around an F-type dwarf (Brahm et al. 2020). The system was first observed by TESS in Sector 6 with 30 min cadence, and then in Sector 33 with a 2 min cadence. The light curve is available on Zenodo (see Data availability).

We used HARPS-N to observe a transit of TOI-892 b. We started observations around UT 19:15 and continued until 03:15 on the night starting on January 1, 2022. The exposure time was set to 960 s, which including overhead comes out to a sampling time of around 987 s. The observations were interrupted due to high humidity.

The results with and without the Gaussian prior on v sin are in agreement, giving $\lambda ={10}_{-35}^{+39}°$ and $9_{-23}^{+30}°$, respectively. Since about half of the transit was not observed, the run including a Gaussian prior on v sin i_⋆ yielded a better constraint on λ, as one would expect.

We analyzed the Doppler shadow and found $\lambda =-{12}_{-16}^{+18}°$, which differs from the RV-RM results more than one would expect for two analyses of the same data. Since the the shadow shown in Fig. A.14 looks convincing, we decided to adopt the value of $-{12}_{-16}^{+18}°$ as our best estimate of À for this system.

A.11 TOI-1130

TOI-1130 is а K dwarf with two transiting planets, an inner Neptune and an outer Jupiter, discovered by Huang et al. (2020). We used ESPRESSO to observe a transit of the larger planet, TOI-1130 c, on the night starting on June 1, 2021. The observations were carried out from around UT 01:30 to 06:00 on June 2, 2021. We opted for an exposure time of 660 s, yielding a sampling time of 701 s. TOI-1130 was observed by TESS in Sector 13 in a 30 minute cadence and in Sectors 27 and 67 in a 20 second cadence. The TESS light curve for planet с is available on Zenodo (see Data availability).

Fig. A.13 RV-RM effect for TOI-892. Same as for Fig. A.1, but with RVs from HARPS-N.

Fig. A.14 Doppler shadow for TOI-892. Same as in Fig. A.3, but for TOI-892.

Because the planets in TOI-1130 are in or near a 2:1 mean-motion resonance, the planets show relatively large transit timing variations (TTVs). The TTV amplitudes are on the order of a few hours for planet b and 15 min for planet с (Korth et al. 2023). Based on the comprehensive model of the system developed by Korth et al. (2023), the mid-transit time of planet с on June 2, 2021 should have occurred at BJD 2459367.6304 ± 0.0003 d (J. Korth, private communication). We used this result as a Gaussian prior on T₀ in our analysis.

The TTVs also affect the phase-folding process that is part of the construction of a high-S/N light curve. To deal with the TTVs, we identified the transits of both planet b and с in the TESS data, and fitted a model the TESS data alone, holding fixed all of the parameters except for the individual midtransit times. Then, in our RV-RM analysis, we fitted the TESS data jointly with the RV time series while applying Gaussian priors on the mid-transit times with a width of 0.005 d. We also applied a Gaussian prior on the stellar mean density of ρ_⋆ = 2.98 ± 0.18 g cm⁻³ (Korth et al. 2023). Most of the transits of planet b were irrelevant, but there was one instance of overlapping transits of planets b and c.

Our BF analysis gave a lower projected rotation velocity for TOI-1130 (1.3 km s⁻¹) than we found in the literature (4.0 ± 0.5 km s⁻¹, Huang et al. 2020). Our result is more in line with the upper limit of 3 km s⁻¹ reported by Korth et al. (2023). From the RV-RM analysis, we obtained v sin i_⋆ = 1.3 ± 0.5 km s⁻¹ and $\lambda =4_{-6}^{+5}°$ when applying a uniform prior on v sin i_⋆. With a Gaussian prior on v sin i_⋆, we obtained $\lambda =4_{-6}^{+5}°$, which we regard as the best estimate. We need to stress that the accuracy of the result does depend on the prior on T₀ that was obtained from the model of (Korth et al. 2023).

We investigated the subplanetary velocities for TOI-1130 to see if we could properly trace the deformation of the stellar lines for this grazing transit, where (unlike for K2-287) the planet only obscures part or the stellar limb, it proved more difficult to locate the distortion in the first and last few in-transit spectra, as shown in Fig. A.16. We also found it necessary to adopt a Gaussian prior on v sin i_⋆.

Fig. A.15 RV-RM effect for TOI-1130. Same as for Fig. A.1, but with RVs from ESPRESSO.

Fig. A.16 Subplanetary velocities for TOI-1130. In the top panel the distortion from subtracting the in-transit CCF from the out-oſ-transit are shown as colored (wiggly) lines. Lines of corresponding colors with a black outline are (Gaussian) fits to the distortion in order to extract the position of the distortion denoted by the dashed vertical line. The bottom panel shows the location of the distortion plotted around the midtransit time. The squares denote epochs for which the decrease in flux from the transiting planet is less than 0.75%. The solid white (black dashed) line shows a model for an (anti-)aligned configuration.

When fitting the subplanetary velocities using all of the data, the posterior for λ had a peak centered on 0°, and a secondary peak near λ = 180°. We found the secondary peak to be a consequence of including the noisiest data points from the first and last few in-transit spectra. Out of concern over systematic errors in our model for limb-grazing transits, we repeated the fit after omitting the in-transit data that were obtained when the loss of light was smaller than 0.75% (the square data points in Fig. A.16). The result was $\lambda =6_{-9}^{+10}°$. Although excluding data without a firm justification is undesirable, we judge the evidence as indicating that any misalignment is smaller than about 20°.

A.12 WASP-50

WASP-50 is a G9 dwarf hosting a hot Jupiter on a 1.9 d orbit discovered by Gillon et al. (2011). TESS data are available from observations in Sector 4 and 31 with 2 min and 20 sec cadence, respectively. The TESS transit light curve is available on Zenodo (see Data availability).

On the night starting on August 24, 2022, we observed a transit of WASP-50 b with ESPRESSO. The observations started at around UT 05:10 and continued until UT 09:40 on August 25, 2022. An exposure time of 240 s ensured a cadence of 270 s.

The time series is well-sampled and the S/N is very high, exactly the case where we do not expect the results to be sensitive to the prior on v sin i_⋆. Indeed, for this system we found −2.9 ± 1.2° when applying a uniform prior and $-{2.9}_{-1.3}^{+1.2}°$ for a Gaussian prior.

The high-quality data prompted us to investigate the RM effect further. We examined the subplanetary velocitiesm, which were shown in Fig. 1, but in Fig. A.18 we display them in a slightly different manner (similar to the format of Fig. A.6 for K2-287). As was the case for K2-287, we included two jitter terms, one for when the planet’s silhouette is on or near the limb, and another when it is completely within the stellar disk. The result of our analysis using the subplanetary velocities was v sin i_⋆ = 2.77 ± 0.14 km s⁻¹ and λ = −1.8 ± 1.0°. The agreement with the RV-RM result is reassuring, but for our best estimate we chose the result −2.9 ± 1.2° from the RV-RM analysis.

A.12.1 Convective blueshift

The term “convective blueshift” refers to the net blueshift of the disk-integrated stellar spectrum caused by the imbalance in brightness between different components of convective cells. The upwelling gas is hotter, brighter, and blueshifted, while the sinking gas is cooler, fainter, and redshifted. Shporer & Brown (2011) showed that this effect contributes to the anomalous RVs with amplitudes of a few m s⁻¹, and variability on a timescale of several minutes. Within our sample, WASP-50 is the only system for which the quality of the data is high enough to contemplate detecting and modeling the convective blueshift.

Fig. A.17 RM-RV effect for WASP-50. Same as for Fig. A.l, but with RVs from ESPRESSO.

Fig. A.18 Subplanetary velocities for WASP-50. Subplanetary velocities for WASP-50 with color-coding as described in Fig. 1. The dark shaded area correspond to the time the planet is completely within the stellar disk and the light shaded area denote ingress and egress. The errorbars include a jitter term, where two separate terms have been added in quadrature; one for the epochs for the full transit and one for the epochs when the planet is partially outside of the limb of the star. The horizontal black lines in the top panel mark v sin i⋆ from the fit.

We modeled the convective blueshift as described by Shporer & Brown (2011). In the model, the stellar disk is pixelated, and each pixel is assigned an intensity according to a limb-darkening law and a net convective velocity (V_CB) directed toward the center of the star. The radial component of each blueshifted surface element is therefore V_CB cos θ where θ is the angle between the normal to the stellar surface and our line of sight. The value of V_CB is expected to range from about −200 m s⁻¹ for K-type stars to −1000 m s⁻¹ for F-type stars (Dravins & Nordlund 1990a,b; Dravins 1999). The solar value is about -300 m s⁻¹ (Dravins 1987). Since WASP-50 is similar to the Sun, with T_eff = 5400 ± 100 K, a reasonable expectation is V_CB ≈ −200 m s⁻¹.

Fig. A.19 Convective blueshift in WASP-50. Same as the left panel of Fig. A.17, but also showing a model as a red solid line that includes the effect of convective blueshift (VCB = −300 m s−1).

Since the anomalous RV caused by the convective blueshift is smallest near the limb, and WASP-50 b has a high impact parameter (b ≈ 0.7), the overall effect of the convective blueshift is expected to be rather weak, as shown in Fig. A.19. Nevertheless, we tested whether the inclusion of the convective blueshift in our RV-RM model leads to a substantial improvement in the quality of the fit to the data,. We created a grid of λ values between −45° and +45° in steps of 0.25°. For each choice of λ, we computed the Bayesian information criterion (BIC; Schwarz 1978) for RV-RM models with and without the convective blueshift. The standard model has a minimum BIC for λ = −2.9° (as expected given our result of −2.9 ± 1.2°), whereas the model with the convective blueshift has a minimum BIC for λ = −3.9°. The difference in the optimized BIC values is 4.5, favoring the model without the convective blueshift. We calculated the difference in BIC (ABIC) between the two models at each λ. ABIC changes sign around λ = −4.9° meaning that the net effect of including convective blueshift results in slightly lower (more misaligned) values for λ.

Even though we did not find that including the convective blueshift significantly improves the fit, we decided to repeat our RV-RM analysis after including the convective blueshift. For convenience, we did not jointly fit the TESS light curve; we simply held b, R_p /R_⋆ , P, and T₀ fixed at the best-fit values from our previous run above. We applied a Gaussian prior on V_CB centered on −400 m s⁻¹ with a width of 100 m s⁻¹. From this run we obtained λ = −3.5 ± 1.3° again suggesting a slightly more misaligned system. A smaller absolute value for V_CB (probably more in-line with the expected value for WASP-50) would make the difference in λ from the two runs even smaller as such we do not expect convective blueshift to significantly affect the results. Thus, when attempting to measure λ with a precision at the level of one degree, the results depend on how much one trusts the model for the convective blueshift and its assumed amplitude.

A.13 WASP-59

The WASP-59 system was discovered by Hébrard et al. (2013). The system consists of a warm Jupiter, WASP-59 b, orbiting a K-dwarf. We observed a transit of WASP-59 b with HIRES. The observations were carried out on November 23, 2016, starting at UT 05:19 to 09:46 with an exposure time of 787 s. Including overhead the sampling time came out to 835 s. The series of 20 observations spanned a transit.

TESS observed the system in Sector 56 with a cadence of with a cadence of 2 min. Given the long time interval that elapsed between our HIRES observations and the TESS observations, one might be concerned about the extra uncertainty in the transit ephemerides. Fortunately, we also arranged for groundbased photometric observations of a transit a few weeks after our spectroscopic observation. We used KeplerCam, which is mounted on the 1.2 m telescope at the Fred L. Whipple Observatory on Mount Hopkins, Arizona (Szentgyorgyi et al. 2005). We observed with an r′ filter on the night starting that began on the night of December 8, 2016. The light curves are available on Zenodo (see Data availability).

Because the HIRES spectra were imprinted with iodine absorption lines, we did not derive a value for v sin i_⋆ with the BF method. Instead, we adopted the value of 2.3 ± 1.2 km s⁻¹ reported by Hébrard et al. (2013). When using this result as a Gaussian prior, we obtained $\lambda ={16}_{-16}^{+17}°$ and v sin i_★ $={0.61}_{-0.17}^{+0.16}$ km s⁻¹, clearly suggesting that the Gaussian prior is not particularly informative. Indeed assuming a uniform prior on v sin i_★ we got v $\sin i_{\star }={0.58}_{-0.16}^{+0.18}$ km s⁻¹ and $\lambda ={16}_{-16}^{+18}°$, which we adopt as a out final value for λ. The resulting RM curve and posteriors are shown in Fig. A.20.

A.14 WASP-136

WASP-136 b is an inflated hot Jupiter orbiting an F9 subgiant star. The system was discovered by Lam et al. (2017). TESS data are available with 2-min cadence in Sector 29 and 20-sec cadence in Sector 41. The light curve is available on Zenodo (see Data availability).

We observed a transit with FIES on the night starting August 30, 2021, with observations running from UT 23:00 to UT 05:50 August 31, 2021. We opted for an exposure time of 1140 s, yielding a sampling time of 1320 s.

Both of the RV-RM fits resulted in a significant misalignment for the WASP-136 system. We obtained $\lambda =-{57}_{-14}^{+9}°$ when applying a uniform prior on v sin i_★, and $-{38}_{-18}^{+16}°$ when applying a Gaussian prior. We adopt the latter as our best estimate of the projected obliquity.

A.15 WASP-148

The WASP-148 system was discovered by Hébrard et al. (2020). The system consists of two giant planets orbiting a late-G dwarf with strong TTVs. The inner planet, WASP-148 b, is a transiting hot Jupiter with a period of 8.8 days. The star’s projected obliquity with respect to the orbit of planet b was measured by Wang et al. (2022) to be $\lambda =-8_{-10}^{+9}°$.

We observed a transit of WASP-148 b using HARPS-N on the night starting March 21, 2021, with observations starting around UT 01:45 to 06:15 March 22, 2021. The exposure time was set to 900 s, resulting in a sampling time of 920 s. Because the strong TTVs lead to extra uncertainty in the transit ephemerides, we also arranged for simultaneous photometric observations of the transit with the Multicolor Simultaneous Camera for studying Atmospheres of Transiting exoplanets (MuSCAT-2; Narita et al. 2015) mounted on the Telescopio Carlos Sánches, Tenerife, Spain. MuSCAT-2 allows for simultaneous observations in four filters; we used a griz filter set. The light curves, after detrending, are available on Zenodo (see Data availability).When modeling these light curves, we did not use GP regression. The limb-darkening coefficients tabulated by Claret et al. (2013) for these filters and for the appropriate stellar parameters are reproduced in Table A.1.

At first, we fitted the RM effect for this system jointly with the TESS and MuSCAT-2 light curves, allowing the mid-transit time of each transit to be a freely adjustable parameter (as we did for TOI-1130, Section A.11). We also tried fitting only the MuSCAT-2 light curves and not the TESS light curve. In both cases, we found λ = 16 ± 40°, which is consistent with, but less precise than, the measurement reported by Wang et al. (2022).

Almenara et al. (2022) performed a “photodynamical” analysis of the system, in which the light curve is fitted with a model in which all the mutual gravitational interactions between the star and both planets are taken into account. They used the data presented by Hébrard et al. (2020) along with additional groundbased data and TESS data. Given their comprehensive modeling effort, we decided to use the parameters from Almenara et al. (2022) as priors, and fitted only the HARPS-N RV measurements and the MuSCAT-2 light curve. The free parameters relevant to the light curve were the mid-transit time, a/R_⋆ , cos i, R_p /R_⋆ , and the sum-of-limb-darkening coefficients for each filter. The period was held fixed. In addition, we included a Gaussian prior on the projected obliquity using the value of $\lambda =-8_{-10}^{+9}°$ as to incorporate all the information available for the system. Again, we performed two MCMCs obtaining a value of −2 ± 8° for both the run with a uniform and the one with a Gaussian prior on v sin i_⋆, respectively.

As this is a dynamically interesting system, for which a large collection of transit times over a long interval can help to constrain the parameters of the planetary system, we measured the mid-transit time based on a fit solely to the MuSCAT-2 light curves. The results are given in Table A.1.

A.16 WASP-172

WASP-172 is an F star with a bloated hot Jupiter discovered by Hellier et al. (2019). Data from TESS are available with 30- min cadence from Sector 11 and 20-sec cadence from Sectors 37 and 38. The TESS light curve is available on Zenodo (see Data availability).

On the night starting June 1, 2022, we observed a transit using ESPRESSO with observations being carried out from UT 23:00 to UT 06:40. We opted for an exposure time of 900 s, resulting in a sampling time of 940 s. The RVs are shown in Fig. A.23, where three data points are marked with gray points as they were excluded from the analysis, since they were taken at high airmass.

Fig. A.20 RM-RV effect for WASP-59. Same as for Fig. A.1, but with RVs from HIRES.

Fig. A.21 RV-RM effect for WASP-136. Same as for Fig. A.1. Here also with RVs from FIES.

By fitting the RV-RM effect and the TESS light curve, we found λ = 118 ± 6° and 115 ± 6° for λ when applying a uniform prior or a Gaussian prior on v sin i_⋆ , respectively. Given the high quality of the data, and the relatively rapid rotation of the star $\left(v\sin i_{\star }={13.0}_{-0.5}^{+0.4} \text{km} {\text{s}}^{-1}\right)$, we decided to analyze the Doppler shadow. The results are displayed in Fig. A.24. The Doppler- shadow analysis gave $\lambda ={121}_{-13}^{+11}°$, consistent with the RV-RM results. We adopted 115 ± 6° as our best estimate for λ, based on the RV-RM fit with a Gaussian prior on v sin i_⋆.

The ESPRESSO transit spectroscopy presented here was also used by Seidel et al. (2023) in an investigation of the planet’s atmosphere. They reported the detection of sodium and hydrogen absorption, as well as a tentative detection of iron absorption. In their study, a value of λ = 121 ± 13° was adopted to correct for the RM effect, which was taken from our Doppler-shadow analysis prior to the completion of this paper.

A.17 WASP-173A

WASP-173Ab is a hot Jupiter discovered by Hellier et al. (2019). The host star is of spectral type G3V and has a wide-orbiting stellar companion. The planet was also dubbed KELT-22Ab by Labadie-Bartz et al. (2019), who reported an independent discovery of the transits using data from the KELT data. Here, we use the name WASP-173Ab.

TESS observed the system twice. Data with 2-min cadence are available from Sector 2, and with 20-sec cadence from Sector 29. The TESS light curve is available on Zenodo (see Data availability).We observed a transit on the night starting July 23, 2022, with observations beginning at UT 05:40 until UT 10:00. The exposure time was set to 555 s, yielding a sampling time of 590 s. The RVs are shown to the left in Fig. A.25.

Our analysis followed similar steps as for WASP-172 (see the previous section). We performed the usual RV-RM analyses and also an analysis of the Doppler shadow, for which the results are shown in Fig. A.26. The RV-RM results were $\lambda ={16}_{-31}^{+38}°$ when applying a uniform prior on v sin i_★ and ${11}_{-20}^{+32}°$ for a Gaussian prior, while the shadow run resulted in a value for λ of ${25}_{-11}^{+12}°$.

Fig. A.22 WASP-148. Same as for Fig. A.1, but with RVs from HARPS-N.

Fig. A.23 WASP-172. Same as for Fig. A.1, but with RVs from ESPRESSO. The gray data points were excluded in the fit due to high airmass.

The v sin i_⋆ derived from the Doppler-shadow analysis was ${7.2}_{-0.4}^{+0.4}\text{\hspace{0.17em}km} {\text{s}}^{-1}$, which is a bit higher than the results obtained from the RV-RM anlayses $\left({6.6}_{-0.8}^{+0.6} \text{km} {\text{s}}^{-1}\right)$ and a Gaussian prior $\left({6.5}_{-0.4}^{+0.3} \text{km} {\text{s}}^{-1}\right)$. The RV-RM-based values are closer to the value derived from the BF of 6.6 km s⁻¹. As our best estimate of λ, we therefore adopt the value of ${11}_{-20}^{+32}°$ from our run applying a Gaussian prior on v sin i_⋆.

A.18 WASP-186

Schanche et al. (2020) discovered the WASP-186 b, a massive hot Jupiter with an eccentric orbit around a mid-F star. As noted by Schanche et al. (2020), WASP-186 was observed by TESS in Sector 17 and produced data with a 30-min cadence. In addition, TESS observed this system in Sectors 42 and 57 and the data are available with 2-min and 20-sec cadence, respectively. The light curves are available on Zenodo (see Data availability).

Using FIES, we observed a transit on the night October 11, 2021, with observations starting UT 22:30 and continuing until UT 02:50 October 12, 2021. The exposure time was set to 900 s, and the time between mid-exposures was 1080 s. The transit observations are shown in Fig. A.27.

Fig. A.24 Doppler shadow for WASP-172. Same as in Fig. A.3, but for WASP-172.

Fig. A.25 RV-RM effect for WASP-173A. Same as for Fig. A.1, but with RVs from ESPRESSO.

Our two RV-RM fits gave $\lambda =-{10}_{-32}^{+29}°\text{ and }\lambda =-{10}_{-16}^{+18}°$ for the run with a uniform and the one with a Gaussian prior on v sin i_⋆ , respectively. The values are in agreement and the prior on v sin i_⋆ leads to a more precise result for λ. Therefore, we chose $\lambda =-{10}_{-16}^{+18}°$ as our best estimate for WASP-186.

A.19 XO-7

XO-7 b was discovered by Crouzet et al. (2020) as part of the XO project (McCullough et al. 2005). The transiting planet is a hot Jupiter. There is also a massive companion on an orbit with a period of at least 2 years around the G2V host star (Crouzet et al. 2020). TESS data are available with 30-min cadence in Sectors 18, 19, and 20; with 2 min cadence in Sectors 40, 47, 53, 53; and with 20-sec cadence in Sectors 59 and 60. The TESS light curves are available on Zenodo (see Data availability).

We observed a transit of XO-7 b with FIES on the night starting on August 27, 2022, with observations starting around UT 00:15 and continuing until 04:00 August 28, 2022. We opted for an exposure time of 840 s, resulting in a sampling time of 1000 s. The RVs from the transit night are shown in Fig. A.28. The egress was not completely covered, and no observations were made after the egress.

The RV-RM results for λ are $-4_{-23}^{+27}°\text{ and }{0}_{-14}^{+20}°$ for the runs with a uniform and a Gaussian prior on v sin i_⋆, respectively, with the posterior distribution for the latter case shown in Fig. A.28. Despite the incomplete coverage of the transit, the high impact parameter of b ~ 0.7 enabled decent constraints to be placed upon λ.

The posterior for λ, shown on the right side of Fig. A.28, has a tail on the left side extending toward −90°, allowing for the possibility that system is strongly misaligned. This tail is probably caused by the lack of egress/post-egress observations and also shows up in the 2σ contours in the left panel of Fig. A.28. We wondered if this unusual feature is related to the lack of egress and post-egress observations. We investigated the issue by examining the stacked spectral CCFs. The stacked CCFs showed contours that largely agree with the posteriors from our RV-RM. However, a peak is seen at around (v sin i_★ , λ) = (4 km s⁻¹ , −60°) that is more consistent with the negative tail of the posterior from the RV-RM analysis, favoring a more misaligned system. A second and more complete transit observation seems warranted. For now, we adopted $0_{-14}^{+20}°$ as our best estimate, based on the run with a Gaussian prior on v sin i_⋆.

Fig. A.26 Doppler shadow for WASP-173A. Same as in Fig. A.3, but for WASP-173A.

Fig. A.27 RV-RM effect for WASP-186. Same as for Fig. A.1. Here also with RVs from FIES.

A.20 WASP-26

WASP-26 b is a hot Jupiter orbiting an early G-type star, discovered by Smalley et al. (2010). Albrecht et al. (2012) analyzed the RV-RM effect based on data obtained with the Planet Finder Spectrograph (PFS; Crane et al. 2010) mounted at the 6.5 m Magellan II telescope. The RV time series is shown in the left half of Fig. A.29. They reported $\lambda =-{34}_{-26}^{+36}°$.

This case is unlike all of other cases presented in this paper in that we did not obtain any new data. Rather, we wanted to take advantage of TESS data that did not exist at the time of the prior study, and correct an error that was made in the prior study by Albrecht et al. (2012). The error was the adoption of an incorrect prior constraint on the transit duration (T₁₂, in the terminology of Albrecht et al. 2012), which can be seen in their Table 3. This error resulted in an incorrect value of b that affected the results for λ.

Rather than placing a prior on T₂₁, we followed our usual procedure of fitting the TESS light curve jointly with the RVs. We used TESS 20-sec data from Sector 29. The light curve is available on Zenodo (see Data availability). As in Albrecht et al. (2012), we applied a Gaussian prior on v sin i_★ of 2.4 ± 1.5 km s⁻¹, and as usual, we applied Gaussian priors on the macro- and microturbulent velocities of 4.0 ± 1.0 km s⁻¹ and 1.1 ± 1.0 km s⁻¹, respectively. We also applied a Gaussian prior on T₀ taken from Anderson et al. (2011), as did Albrecht et al. (2012). The reason for this choice is that about 10 years elapsed between the spectroscopic transit observations and the TESS observations; small systematic uncertainties in the TESS-based ephemeris might propagate into large uncertainties in the relevant transit time. We found a projected obliquity of $-{16}_{-10}^{+14}°$. Fig. A.29 shows the RV-RM curve and posteriors.

Fig. A.28 RV-RM effect for XO-7. Same as for Fig. A.1. Here also with RVs from FIES.

Fig. A.29 RV-RM effect for WASP-26. Same as for Fig. A.1, but with RVs from PFS.

References

All Tables

All Figures

Fig. 1 Rossiter–McLaughlin effect displayed in four ways, for the illustrative case of WASP-50. Top right: schematic of the planet transiting the stellar disk. The wedge represents the best-fit λ and its uncertainty. The color scheme is based on the star’s rotational Doppler shift; the same color scheme is used in the top left and lower left panels to convey the location of the planet on the star. Top left: spectral lines. The thick black curve is the average out-of-transit cross-correlation function (CCF) and the superimposed white curve is a CCF observed near mid-transit. Middle left: spectral line distortions. Shown are the results of subtracting five of the in-transit CCFs from the mean out-of-transit CCF. Solid curves show the data. Dashed curves are the best-fit Lorentzian profiles, with vertical lines marking their centers. The color coding is done according to the apparent velocity anomaly going from redshifted (in red) to blueshifted (in blue). Bottom left: Doppler shadow. Each row shows the spectral line distortion at a particular time. The white contours are amplitude levels of the distortion. Each small circle marks the best-fit subplanetary velocity (the Doppler shift of the portion of the star directly behind the planet’s center), with colors that convey the corresponding anomalous radial velocity shown in the lower right panel. Dark vertical lines mark the measured v sin i★, and the diagonal white line is the best-fit trajectory. Bottom right: anomalous radial velocities. The dark shaded area is the range of times when the planet’s silhouette is completely contained within the stellar disk.

In the text

	Fig. 2 Obliquities and projected obliquities as a function of Teff (left), a/R* (center), and Mp/M* (right). Triangles indicate a measurement whose value is outside of the plotted range. In all panels, the color of each data point conveys Teff using the color scale shown in the lower left panels. The contours are KDEs illustrating the density of measurements in a given parameter space.
In the text

	Fig. 3 Hierarchical Bayesian inference of the obliquity distribution using the code of Dong & Foreman-Mackey (2023). The left panel is for all systems for which λ was measured, the middle panel is for systems for which both λ and i★ were measured but without using the i★ information, and the right panel is for systems for which both λ and i★ were measured and making use of the i★ information.
In the text

Fig. 4 Hierarchical Bayesian inference of the obliquity distribution for subsets of the sample that are suspected of being especially prone to misalignment. The left panel is for all sub-Saturns with λ measurements, without using any i★ information even when it is available. The middle panel is for all hot Jupiters around hot stars with λ measurements, without using any i★ information even when it is available. The right panel is for hot Jupiters around hot stars for which both λ and i★ are known and utilized. In the right panel, the black upward tick marks on the horizontal axis are the individual values of cos ψ, and the horizontal bands surrounding the tick marks convey the uncertainties. The solid gray curve is a KDE estimated from these measurements.

In the text

	Fig. 5 Hierarchical Bayesian inference of the obliquity distribution for the same subsets as in Fig. 4, but here a uniform prior, 𝒰(–4, 10), was applied to log κ in the "DFM" code.
In the text

Fig. 6 Obliquity distribution for stars hosting planets on eccentric and circular orbits. Left: the cumulative sum of |λ| for eccentric planets shown in red, sorted by Teff. The black line is the average of 5000 randomly sampled sets from the circular population with the gray shaded area denoting the 1σ and 2σ intervals. The vertical blue line denotes the Kraft break (Teff = 6100 K). Right: histogram showing the distribution of the circular population and the value of the eccentric population at the Kraft break. The Gaussian is a fit to the histogram. The circular and eccentric samples only differ by 2.1 σ, meaning the eccentric population does not appear to be significantly more misaligned. Adapted from Rice et al. (2022).

In the text

	Fig. 7 Obliquity distribution for stars hosting planets on close-in and wide orbits. Same as Fig. 6, but for systems with planets on close-in or wide orbits with the division given for a/R★ = 12. The systems with planets on wider orbits appear to be significantly more misaligned compared to those with planets closer in.
In the text

Fig. 8 Precise projected obliquity measurements. Here we highlight the most precise (σλ < 2°) measurements of λ for systems with cool stellar hosts (Teff < 6100 K). WASP-50 is highlighted with a black center. The color-coding is done according to Teff and the marker sizes scale with the tidal alignment timescale (as given by Zahn 1977). The gray shaded area shows the Gaussian distribution inferred from cool hosts (Teff < 6100 K) with close-in (a/R★ < 10), massive (Mp < 0.3 MJ) planets. The horizontal lines demarcate the Solar System obliquity of 7.155° (Beck & Giles 2005). The light gray error bars show the less precise measurements. Adapted from Albrecht et al. (2022, Figure 10).

In the text

Fig. 9 Obliquities of stars with multiple transiting planets, as constrained by either the RM effect or asteroseismology. Each horizontal line shows the a/R⋆ values of the planets in a given system, with the solar system on top (with gaps to allow Jupiter and Neptune to be shown). For the planets, the size of each circle conveys the planet’s radius relative to Jupiter. For the stars, symbol size conveys the radius relative to the Sun, with a color conveying the effective temperature. Gray circles are planets for which the RM effect has been measured – except for β Pic b where spectro-interferometry was used. The wedge indicates the projected obliquity for a given system, based on the most precise measurement for any planet in the system, with the exception of HD 3167 for which two planets appear to have very different inclinations. The asteroseismic measurements constrain i⋆ and not λ; in those cases the wedge should therefore be read as aligned (Kepler-50 and Kepler-65) or misaligned (Kepler-56, Kepler-129). Planetary (tan circle) or stellar companions (red star, M-dwarfs) to a given system are shown to the right, where √(÷) denotes that the companion could (not) have influenced the obliquity of the inner planetary system. The asterisks for LTT 1445A and TOI-451A denote that these measurements should be taken with a grain of salt. Adapted from Wang et al. (2022, Figure 3).

In the text

Fig. A.1 RM effect observed in HD 118203. Left: Time series and residuals of the FIES+ RVs (orange) and FIES RVs (green). The black curve is the best-fit model and the gray shading shows the 1 and 2σ confidence intervals created by randomly drawing values from the MCMC chains for b, v sin i⋆, Rp/R⋆, and λ. The dashed curve is the best-fit model assuming λ = 0°. Gray circles are the data that were obtained during poor atmospheric conditions and were excluded from the analysis. Right: Correlations between λ, v sin i⋆, and b. The red lines illustrate the prior that was applied (either uniform or Gaussian).

In the text

	Fig. A.2 RV-RM effect for HD 148193. Same as for Fig. A.1, but with RVs from HARPS-N.
In the text

Fig. A.3 Doppler shadow for HD 148193. Left: The observed Doppler shadow created by subtracting the in-transit CCFs from an average out-of-transit CCF. Middle: Model Doppler shadow. Right: Residuals from subtracting the model from the observed shadow. The summed residuals are given in the adjacent panels. The color bar demarcate the contrast. Green horizontal lines denote points of first (solid), second (dashed), third (dashed), and fourth (solid) contact. The solid, orange, vertical lines denote ±v sin i⋆.

In the text

	Fig. A.4 RV-RM effect for K2-261. Same as for Fig. A.1, but with RVs from HARPS-N.
In the text

	Fig. A.5 RV-RM effect for K2-287. Same as for Fig. A.l, but with RVs from ESPRESSO.
In the text

Fig. A.6 Subplanetary velocities for K2-287. Subplanetary velocities for K2-287 with color-coding as described in Fig. 1. The dark shaded area corresponds to the time the planet is completely within the stellar disk and the light shaded area denote ingress and egress. The error-bars include a jitter term, where two separate terms have been added in quadrature; one for the epochs for the full transit and one for the epochs when the planet is partially outside of the limb of the star. The horizontal black lines in the top panel mark v sin i⋆ from the fit.

In the text

	Fig. A.7 KELT-3. Same as for Fig. A.1. Here also with RVs from FIES.
In the text

	Fig. A.8 Peak of the stacked CCFs of KELT-3. The peak of the stacked CCFs of KELT-3 created by assuming the best-fitting value of b from the RV-RM fit. The gray contours show a peak around (v sin i⋆,λ)=(6.0 km s−1, 1.0°). The red contours are the same as created from the posterior shown to the right in Fig. A.7.
In the text

	Fig. A.9 KELT-4A. Same as for Fig. A.1. Here also with RVs from FIES.
In the text

	Fig. A.10 LTT 1445A. Same as for Fig. A.1, but with RVs from HARPS-N.
In the text

	Fig. A.11 RV-RM effect for TOI-451. Same as for Fig. A.1, but with RVs from ESPRESSO.
In the text

	Fig. A.12 RV-RM effect for TOI-813. Same as for Fig. A.1, but with RVs from ESPRESSO.
In the text

	Fig. A.13 RV-RM effect for TOI-892. Same as for Fig. A.1, but with RVs from HARPS-N.
In the text

	Fig. A.14 Doppler shadow for TOI-892. Same as in Fig. A.3, but for TOI-892.
In the text

	Fig. A.15 RV-RM effect for TOI-1130. Same as for Fig. A.1, but with RVs from ESPRESSO.
In the text

Fig. A.16 Subplanetary velocities for TOI-1130. In the top panel the distortion from subtracting the in-transit CCF from the out-oſ-transit are shown as colored (wiggly) lines. Lines of corresponding colors with a black outline are (Gaussian) fits to the distortion in order to extract the position of the distortion denoted by the dashed vertical line. The bottom panel shows the location of the distortion plotted around the midtransit time. The squares denote epochs for which the decrease in flux from the transiting planet is less than 0.75%. The solid white (black dashed) line shows a model for an (anti-)aligned configuration.

In the text

	Fig. A.17 RM-RV effect for WASP-50. Same as for Fig. A.l, but with RVs from ESPRESSO.
In the text

Fig. A.18 Subplanetary velocities for WASP-50. Subplanetary velocities for WASP-50 with color-coding as described in Fig. 1. The dark shaded area correspond to the time the planet is completely within the stellar disk and the light shaded area denote ingress and egress. The errorbars include a jitter term, where two separate terms have been added in quadrature; one for the epochs for the full transit and one for the epochs when the planet is partially outside of the limb of the star. The horizontal black lines in the top panel mark v sin i⋆ from the fit.

In the text

	Fig. A.19 Convective blueshift in WASP-50. Same as the left panel of Fig. A.17, but also showing a model as a red solid line that includes the effect of convective blueshift (VCB = −300 m s−1).
In the text

	Fig. A.20 RM-RV effect for WASP-59. Same as for Fig. A.1, but with RVs from HIRES.
In the text

	Fig. A.21 RV-RM effect for WASP-136. Same as for Fig. A.1. Here also with RVs from FIES.
In the text

	Fig. A.22 WASP-148. Same as for Fig. A.1, but with RVs from HARPS-N.
In the text

	Fig. A.23 WASP-172. Same as for Fig. A.1, but with RVs from ESPRESSO. The gray data points were excluded in the fit due to high airmass.
In the text

	Fig. A.24 Doppler shadow for WASP-172. Same as in Fig. A.3, but for WASP-172.
In the text

	Fig. A.25 RV-RM effect for WASP-173A. Same as for Fig. A.1, but with RVs from ESPRESSO.
In the text

	Fig. A.26 Doppler shadow for WASP-173A. Same as in Fig. A.3, but for WASP-173A.
In the text

	Fig. A.27 RV-RM effect for WASP-186. Same as for Fig. A.1. Here also with RVs from FIES.
In the text

	Fig. A.28 RV-RM effect for XO-7. Same as for Fig. A.1. Here also with RVs from FIES.
In the text

	Fig. A.29 RV-RM effect for WASP-26. Same as for Fig. A.1, but with RVs from PFS.
In the text