Evidence for Primordial Alignment: Insights from Stellar Obliquity Measurements for Compact Sub-Saturn Systems

Despite decades of effort, the mechanisms by which the spin axis of a star and the orbital axes of its planets become misaligned remain elusive. Particularly, it is of great interest whether the large spin-orbit misalignments observed are driven primarily by high-eccentricity migration -- expected to have occurred for short-period, isolated planets -- or reflect a more universal process that operates across systems with a variety of present-day architectures. Compact multi-planet systems offer a unique opportunity to differentiate between these competing hypotheses, as their tightly-packed configurations preclude violent dynamical histories, including high-eccentricity migration, allowing them to trace the primordial disk plane. In this context, we report measurements of the sky-projected stellar obliquity ($\lambda$) via the Rossiter-McLaughlin effect for two sub-Saturns in multiple-transiting systems: TOI-5126 b ($\lambda=1\pm 48 ^\circ$) and TOI-5398 b ($\lambda=-8.1^{+5.3 \circ}_{-6.3}$). Both are spin-orbit aligned, joining a fast-growing group of just three other compact sub-Saturn systems, all of which exhibit spin-orbit alignment. In aggregate with archival data, our results strongly suggest that sub-Saturn systems are primordially aligned and become misaligned largely in the post-disk phase, as appears to be the case increasingly for other exoplanet populations.


INTRODUCTION
Stellar obliquity-i.e., the angle between the spin axis of the host star and the net orbit normal axis of the planets-is a powerful probe of a planetary system's formation history.While the precise mechanisms driving spin-orbit misalignment are not yet clear, it is important to distinguish whether misalignment is largely a product of high-eccentricity migration, which should be confined to short-period isolated planets 1 (e.g., hot Jupiters; see Dawson & Johnson 2018 and references therein), or some universal process that operates across planetary systems with a variety of architectures, e.g., magnetic warping (Lai et al. 2011;Romanova et al. 2013), interactions with stellar com-panions (Lubow & Ogilvie 2000;Batygin 2012), stellar flybys (Hao et al. 2013), or stellar gravity waves (Rogers et al. 2012(Rogers et al. , 2013)).
Compact multi-planet systems may be critical to differentiate between these two scenarios, as their tight orbital configurations are inconsistent with high-eccentricity migration induced by planet-planet scattering (Rasio & Ford 1996;Chatterjee et al. 2008), Lidov-Kozai cycling (Wu & Murray 2003;Fabrycky & Tremaine 2007;Naoz 2016), or secular interactions (Naoz et al. 2011;Wu & Lithwick 2011;Petrovich 2015) -all of which tend to disrupt the orbits of nearby planetary neighbors (Mustill et al. 2015).As such, if higheccentricity migration is the primary driver of misalignment, compact multi-planet systems will retain their primordially aligned configurations.On the other hand, if universal misalignment processes dominate, compact systems may be misaligned at a rate comparable to isolated systems.
The great majority of spin-orbit angle measurements have been made for close-in Jupiters, particularly hot Jupiters (see recent review by Albrecht et al. 2022), which rarely host nearby companions (Steffen et al. 2012;Huang et al. 2016;Hord et al. 2021;Wu et al. 2023).Therefore, it is useful to expand the census of stellar obliquities to other types of exoplanets.Sub-Saturns, broadly defined as giant planets with masses ranging between that of Neptune (∼17 M⊕) and Saturn (∼95 M⊕), represent a particularly interesting population for such studies.The core accretion model implies that they are 'failed gas giants' (Pollack et al. 1996), suggesting they have analogous origins to Jupiters (Dong et al. 2018;Lee 2019;Hallatt & Lee 2022).However, sub-Saturns are less amenable to detection than Jupiters and are relatively understudied as an exoplanet population, so few constraints on their formation histories exist.As of this writing, only about two dozen sub-Saturns have had their spin-orbit angle measured, demonstrating that these planets can take on a range of values from near-perfect alignment to almost completely retrograde.Yet, only three of these measurements have been made for sub-Saturns in compact multi-planet systems: Kepler-9 b (Wang et al. 2018), AU Mic b (Hirano et al. 2020a), and WASP-148 b (Wang et al. 2022).While this sample is limited in number, all three systems show evidence of spin-orbit alignment and thus provide tentative support for primordial alignment.
The Transiting Exoplanets Survey Satellite (TESS) mission (Ricker et al. 2015), now well into its Extended Mission phase, has produced a growing sample of sub-Saturns around bright stars.Of particular interest to this work, TESS is responsible for the discovery of at least two sub-Saturns in compact multi-planet systems: TOI-5126 b (Fairnington et al. 2023), andTOI-5398 b (Mantovan et al. 2022;Mantovan et al. 2024a).Given the brightness of their host stars and advances in Extreme Precision Radial Velocity (EPRV), these systems are amenable to stellar obliquity measurements.
In this study, we present Rossiter-McLaughlin (RM) effect (Holt 1893;Rossiter 1924;McLaughlin 1924;Queloz et al. 2000) measurements for both confirmed compact TESS sub-Saturn systems, TOI-5126 and TOI-5398, using the highprecision NEID spectrograph (Schwab et al. 2016) on the WIYN 3.5 m telescope at Kitt Peak, AZ.TOI-5126, a relatively bright (V = 10.1 mag) late-type F dwarf (T eff = 6297 K), hosts a hot sub-Saturn with a 5.5-day period and an outer warm sub-Neptune with a 17.9-day period (Fairnington et al. 2023).TOI-5398, a bright (V = 10.1 mag) early-type G dwarf (T eff = 6039 K), harbors a warm sub-Saturn with a 10.6-day period and an inner hot sub-Neptune with a 4.8day period (Mantovan et al. 2022;Mantovan et al. 2024a).Our reported spin-orbit angles for these two systems represent a significant contribution to the limited sample of sub-Saturns in compact multi-planet systems with such measurements.This work is additionally part of the 11th published outcome of the Stellar Obliquities in Long-period Exoplanet Systems (SOLES) survey, as detailed in a series of publications (Rice et al. 2021;Wang et al. 2022;Rice et al. 2022;Rice et al. 2023a;Hixenbaugh et al. 2023;Dong et al. 2023;Wright et al. 2023;Rice et al. 2023b;Lubin et al. 2023;Hu et al. 2024).
This paper is organized as follows.In Section 2, we outline our radial velocity observations.In Section 3, we describe our determination of stellar parameters.In Section 4, we detail our modeling of the RM signals and, subsequently, the stellar obliquities.In Section 5, we quantify the significance of our results and discuss their implications for the origins of spin-orbit misalignment and the formation histories of sub-Saturns.

OBSERVATIONS
The RM effect of TOI-5126 and TOI-5398 were observed using the WIYN/NEID spectrograph in the High Resolution (HR) mode (R ∼ 110, 000, Halverson et al. 2016;Schwab et al. 2016) on February 28, 2023and April 16, 2023, respectively.NEID is a fiber-fed (Kanodia et al. 2018), actively environmentally stabilized spectrograph (Stefansson et al. 2016;Robertson et al. 2019) with a wavelength coverage of 380 nm to 930 nm.For TOI-5126 b, we obtain 20 radial velocity (RV) measurements in HR mode with 1000second exposures from 02:20-07:48 UT.In total, 14 of these RVs cover the full transit while the remaining six provide out-of-transit baseline; one occurring 0.3 hr pre-ingress and the other five spanning nearly 1.5 hr post-egress.These observations occurred under atmospheric conditions with a seeing range of 1.2 ′′ -2.8 ′′ (median 1.75 ′′ ) and an airmass range of z = 1.03-2.31.At a wavelength of 5500Å, the NEID spectrograph achieved a signal-to-noise ratio (SNR) of 56 pixel −1 .For TOI-5398 b, we obtain 14 RV measurements in HR mode with 1250-second exposures from 02:55-09:04 UT, eight of which were in-transit, sampling nearly the entire transit (excluding ingress).The conditions featured a seeing range of 0.7 ′′ -1.5 ′′ (median 1.0 ′′ ) and an airmass range of z = 1.08-1.71.At the same wavelength, the spectrograph's SNR was 60 pixel −1 .We obtain an additional 14 NEID RVs covering TOI-5398 b's orbit interspersed between April 8, 2022 and December 27, 2023, four with 240-second exposure times and nine with 300-second exposure times, and we utilize these measurements in subsequent global fitting for this system.
The NEID spectra were analyzed using version 1.3.0 of the NEID Data Reduction Pipeline (NEID-DRP) 2 , and ra-dial velocities were derived via the cross-correlation function (CCF) method within NEID-DRP.We extracted the resulting barycentric-corrected RVs (denoted as CCFRVMOD within the NEID-DRP documentation), from the NExScI NEID Archive 3 .The RV data gathered from NEID for this study are illustrated in Figure 1 (note that we subtract the Keplerian baseline to more clearly demonstrate the RM effect) and are available through the Data Behind the Figure program (see figure caption for link).

Synthetic Spectral Fit By iSpec
All out-of-transit NEID spectra for TOI-5126 and TOI-5398 are corrected for RV shifts, co-added, and subsequently utilized to determine the respective stellar atmospheric parameters, including stellar effective temperature (T eff ), surface gravity (log g * ), metallicity ([Fe/H] ), and projected rotational velocity (v sin i * ).We use the synthetic spectral fitting technique provided by the Python package iSpec (Blanco-Cuaresma et al. 2014;Blanco-Cuaresma 2019) to measure these parameters.
We employ the SPECTRUM radiative transfer code (Gray & Corbally 1994), the MARCS atmosphere model (Gustafsson et al. 2008), and the sixth version of the GES atomic line list (Heiter et al. 2021), all integrated within iSpec, to create a synthetic model for all the combined out-of-transit NEID spectra (SNR of 112 for TOI-5126; SNR of 137 for TOI-5398).In our fitting process, we treat micro-turbulent velocities as a variable, allowing us to accurately represent the small-scale turbulent motions in the stellar atmosphere.Conversely, we determine macro-turbulent velocities using an empirical relationship (Doyle et al. 2014) that leverages established correlations with various stellar attributes.We select specific spectral regions to streamline the fitting process, focusing on the wings of the Hα, Hβ, and Mg I triplet lines, which are indicative of T eff and log g * , as well as Fe I and Fe II lines, which are crucial for constraining [Fe/H] and v sin i * .We then utilize the Levenberg-Marquardt nonlinear least-squares fitting algorithm (Markwardt 2009) to refine our spectroscopic parameters by iteratively minimizing the χ 2 value between the synthetic and observed spectra.The final spectroscopic parameters are detailed in Table 1.

SED+MIST Fit By EXOFASTv2
To ascertain additional stellar parameters, such as stellar mass (M * ) and radius (R * ), we utilize the MESA Isochrones & Stellar Tracks (MIST) model (Choi et al. 2016;Dotter 2016) in combination with a spectral energy distribution (SED) fit.Photometry for the SED fit was compiled from 3 https://neid.ipac.caltech.edu/various catalogs, including 2MASS (Cutri et al. 2003), WISE (Cutri et al. 2021), TESS (Ricker et al. 2015), and Gaia DR3 (Gaia Collaboration et al. 2023).We apply Gaussian priors based on our synthetic spectral fit to T eff and [Fe/H] , along with the parallax from Gaia DR3, and an upper limit for the V -band extinction from Schlafly & Finkbeiner (2011).To accommodate the approximate 2.4% systematic uncertainty floor in T eff , as suggested by Tayar et al. (2022), we increase the uncertainties of T eff to 150 K for both TOI-5126 and TOI-5398.
As mentioned in Fairnington et al. (2023) and Mantovan et al. (2024a) respectively, both TOI-5126 and TOI-5398 exhibit strong stellar rotation signals, providing additional stringent constraints on their stellar ages.We adopt Equation 7from Barnes (2007) to calculate the gyrochronological ages for these stars based on the stellar rotation periods derived from our work (see Section 4).The resulting gyrochronological ages for TOI-5126 and TOI-5398 are 0.13±0.02Gyr and 0.42 +0.07 −0.05 Gyr, respectively.Subsequently, we adopt these ages and their 3σ uncertainties as priors for the stellar ages in our SED fit.
We perform the SED fitting using the Differential Evolution Markov Chain Monte Carlo (DEMCMC) technique, integrated within EXOFASTv2 (Eastman 2017; Eastman et al. 2019), from which we obtain uncertainty estimates.The DEMCMC procedure was considered converged when the Gelman-Rubin diagnostic ( R; Gelman & Rubin 1992) fell below 1.01 and the count of independent draws surpassed 1000.Our final adopted stellar parameters are listed in Table 1; we find all derived values are consistent (within < 2σ) with results from discovery papers (see Fairnington et al. 2023for TOI-5126, Mantovan et al. 2024a for TOI-5398).

OBLIQUITY MODELING
We determine the sky-projected spin-orbit angle λ for TOI-5126 b and TOI-5398 b by jointly fitting transit photometry from their respective discovery papers (including twominute cadence TESS light curves), and all in-transit and out-of-transit RVs from our NEID observations (see Section 2), using a modified version of the allesfitter package (Günther & Daylan 2021).The modified version of allesfitter incorporates the Hirano et al. (2011) RM model implemented in tracit (Hjorth et al. 2021;Knudstrup & Albrecht 2022).We exclude the transit data for the non-sub-Saturn planet in both systems but model their RVs (TOI-5126 c at P = 17.9 d and TOI-5398 c at P = 4.77 d).
For TOI-5126 b, we adopt the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP) light curves from the TESS Science Processing Operations Center (SPOC), modeling eight full transits and one partial transit from  Sector 45, four full transits and one partial transit from Sector 46, and four full transits from Sector 48.Within allesfitter, we simultaneously fit the TESS and NEID data together with all photometric data utilized in TOI-5126 b's discovery paper (Fairnington et al. 2023), including observations from the CHaracterising ExOplanets Satellite (CHEOPS, Benz et al. 2021) as well as ground-based 1.0 m telescopes at McDonald Observatory in Fort Davis, TX and at the Cerro Tololo Inter-American Observatory (CTIO) in Cerro Tololo, Chile, which are both part of the Las Cumbres Observatory Global Network (LCOGT; Brown et al. 2013).Due to their large scatter relative to the expected planetary RV semi-amplitudes, we do not make use of the out-of-transit RVs presented in Fairnington et al. (2023), which were obtained from the Tillinghast Reflector Echelle Spectrograph (TRES) at the Fred Lawrence Whipple Observatory (FLWO) and the CHIRON spectrograph on the SMARTS telescope at CTIO.
For TOI-5398 b, we also consider the PDCSAP TESS light curves (we did not find substantive evidence for the overcorrections in the PDCSAP data reported by Mantovan et al. 2024a), modeling two full transits from TESS Sector 48, as well as photometry from the LCOGT 1.0 m telescopes at McDonald Observatory and Teide Observatory and from the 1.2 m KeplerCam at FLWO on Mount Hopkins, Arizona (all sourced from Mantovan et al. 2024a).We additionally incorporate the out-of-transit RVs from Mantovan et al. (2024a), which were obtained from the High Accuracy Radial velocity Planet Searcher (HARPS-N, Cosentino et al. 2012) at the Italian Telescopio Nazionale Galileo (TNG).Recently, Mantovan et al. (2024b) published an independent RM effect measurement of TOI-5398 b; therefore, we model the RM effect for this target both with and without the inclusion of their published in-transit RVs (holding all other factors equal), which were also obtained from the HARPS-N spectrograph4 .To account for strong stellar activity in the RVs, we apply a rotational Gaussian Process (GP) regression kernel formulated by Foreman-Mackey et al. (2017), where P rot represents the period of stellar rotation, L denotes the coherence timescale, τ signifies the interval between two successive data points, B is a scaling hyperparameter that adjusts the amplitude, and C serves as the balance parameter for the periodic and nonperiodic components of the GP kernel.
For TOI-5398, we note that our adopted GP model is most similar to "Case 1" from Mantovan et al. (2024a), though the results from their multi-dimensional "Case 3" were ultimately adopted.
For both TOI-5126 b and TOI-5398 b, we adopt priors but utilize the results reported in the literature (Fairnington et al. 2023 andMantovan et al. 2024a, respectively) as our initial guesses for R p /R * , (R * + R p )/a, cos(i), T 0 , P , and K (see Table 1 for parameter definitions).We fixed eccentricity (e) and the argument of periastron (ω) to 0 due to the nearzero eccentricity reported by the discovery papers of TOI-5126 b and TOI-5398 b.We consider transformed linear and quadratic limb-darkening coefficients q 1 and q 2 , in addition to physical linear and quadratic limb-darkening coefficients u 1 and u 25 , for both NEID and HARPS-N, initializing q 1 and q 2 at 0.5 for each instrument.As an initial guess for v sin i * , for which we apply Gaussian priors, we adopt the values derived from our synthetic spectral fit.Finally, we initialize λ at 0 • for both systems.All parameters (except e and ω) are allowed to vary during the fitting process.To account for the short-term overnight instrumental systematics and stellar variability in the RM fit, we employ a quadratic baseline for TOI-5126 and a constant baseline for TOI-5398, though we note that all tested baselines (constant, linear, quadratic, and cubic) produced consistent λ values (within 1σ).
For the fits of both systems, we sample the posterior distributions for all modeled parameters using an affine-invariant Markov Chain Monte Carlo (MCMC) grounded in emcee (Foreman-Mackey et al. 2013a) with 100 walkers.Best-fit parameters were obtained after 200,000 accepted steps per walker were reached (note that each ran 10,000 burner steps); our final MCMC results and associated 1σ uncertainties are summarized in Table 1.Additionally, all Markov Chains reached > 50× their autocorrelation lengths, indicating convergence.
We find good agreement (< 2σ difference) with all stellar, planetary, and derived parameters reported in the discovery papers of both systems (Fairnington et al. 2023 for TOI-5126 b and Mantovan et al. 2024a for TOI-5398 b), except T 0;b (2.5σ discrepant) for TOI-5126, which may be due to tentative transit timing variations (see Fairnington et al. 2023).Critically, our RM fits show that both sub-Saturns are consistent with spin-orbit alignment.For TOI-5126 b, we find a best-fit projected obliquity of λ = 1 ± 48 • .For TOI-5398 b, we compute a best-fit λ = −8.1 +5.3 −6.3 • using in-transit RVs from both NEID (presented in this work) and HARPS-N (presented in Mantovan et al. 2024b), while we derive a marginally aligned value of λ = −24 +14 −13 • using our NEID data alone (note that both values are within 2σ of the λ = 3.0 +6.8 −4.2 reported by Mantovan et al. 2024b).Figure 1 reveals the best-fit joint models and associated residuals for both systems, while Table 1 contains all fitted and derived parameters.
It is useful to additionally measure the 3D stellar obliquity ψ, which yields a system's true spin-orbit configuration.
To estimate ψ, we first derive the stellar inclination i * and stellar equatorial velocity v eq for both TOI-5126 and TOI-5398 by applying the Bayesian inference method presented by Masuda & Winn (2020) and Hjorth et al. (2021) to our values for stellar radius R * , stellar rotation period P rot , and cos i * .We utilize the aforementioned TESS PDCSAP light curves to compute the stellar rotation periods for each target, employing the autocorrelation function (ACF) implemented in SpinSpotter (Holcomb et al. 2022).Correspondingly, we find P rot = 4.05 ± 0.21 days for TOI-5126 and P rot = 7.443±0.041days for TOI-5398.However, due to the effects of latitudinal differential rotation (Epstein & Pinsonneault 2014;Aigrain et al. 2015) and the systematic uncertainty floor (of ≈ 4.2%) on the derivation of stellar radii (Tayar et al. 2022), we adopt a 10% uncertainty floor on our derived stellar rotation periods, yielding P rot = 4.05 ± 0.41 days and P rot = 7.443 ± 0.74 days for TOI-5126 and TOI-5398, respectively.These estimates are consistent with those reported by their respective discovery papers, which both applied the Lomb-Scargle periodogram to photometric data, rather than the ACF, in order to ascertain their adopted rotation periods (Fairnington et al. 2023 find P rot = 4.602 +0.071 −0.067 days for TOI-5126, Mantovan et al. 2024a find P rot = 7.34 ± 0.05 days for TOI-5398).From our Bayesian analysis, we find i * = 90 ± 15 • and v eq = 14.36 ± 1.02 km/s for TOI-5126 and i * = 91 ± 17 • and v eq = 8.09 ± 0.87 km/s for TOI-5398, and compute the true 3D stellar obliquities as (see Eqn. 9 of Fabrycky & Winn 2009): where i is the inclination angle of the planet.For TOI-5126 b, we find a true obliquity of ψ = 37.1 +20.3 −33.6 • .Similarly, we derive ψ = 16.4 +7.6 −10.3 • for TOI-5398 b (in good agreement with the ψ = 13.2 ± 8.2 • derived by Mantovan et al. 2024b), demonstrating that both systems are consistent with alignment.

DISCUSSION
While sub-Saturns broadly appear to boast a wide range of stellar obliquities, the few previously confirmed to be in compact multi-planet systems are exclusively aligned.Our RM measurements for two of such recently-confirmed sub-Saturns, TOI-5126 b and TOI-5398 b, continue this trend, as both are consistent with alignment.We show the updated stellar obliquity distribution for the sub-Saturn population as a function of stellar effective temperature in Figure 2. To construct this sub-Saturn sample, we combine the catalogs of Albrecht et al. (2022) and TEPCat 6 (Southworth 2011) accessed on March 7, 2024.We consider only RM effect measurements, which exhibit near-uniform sensitivity to the full range of possible λ values (compared to nonuniform methods such as starspot tracking or gravity darkening, Siegel et al. 2023;Dong & Foreman-Mackey 2023) and constitute the majority of obliquity measurements (Triaud 2018;Albrecht et al. 2022), prioritizing those selected in Albrecht et al. ( 2022) (we otherwise adopt the "preferred" measurements from TEPCat).Finally, we apply a mass cut of 17 M ⊕ < M pl < 95 M ⊕ , and select only single-star systems by cross-matching our sample with the multi-star catalog of Rice et al. (2024), which is based on the most recent data from Gaia DR3.This yields 22 sub-Saturn systems with λ measurements to comprise our statistical sample.We further distinguish between sub-Saturns that are isolated and those in compact multi-planet systems, which we define as having a nearby companion with a period ratio of < 4. All sub-Saturns considered in this work are presented in Table 2 (note that we supplement any missing parameters using the NASA Exoplanet Archive7 and values adopted by their respective λ reference papers).Our RM measurements allow for more reliable statistical analyses to be performed on the relationship between compact configurations and spin-orbit alignment for sub-Saturns (Section 5.1) and thus more robust constraints on the origins of spin-orbit misalignment as well as the formation histories of the sub-Saturn population (Sections 5.2 and 5.3).

Sub-Saturns in Compact Multi-planet Systems are Significantly Aligned
It is clear from Figure 2 that spin-orbit alignment is much more common for sub-Saturns in compact multi-planet systems than for isolated sub-Saturns, while no clear relation exists for the host star's effective temperature (as seen in previous works, e.g., Albrecht et al. 2022).Our two stellar obliquity measurements increase the sample of sub-Saturns in compact configurations to five, enabling the first statistical tests of this trend.We investigate the statistical significance of preferential alignment in compact systems by comparing the number of misaligned systems found in the compact multi-planet sub-Saturn sample (zero) with that found in a series of random draws from the isolated sub-Saturn sample.We perform this statistical analysis using the projected 2D spin-orbit angles λ, rather than the true 3D spin-orbit angles ψ, for two reasons: i) considering solely ψ is statistically prohibitive, as only 14 of the 22 systems in our sample have ψ measurements, and ii) nearly all 22 systems in our sample orbit cooler stars (T eff < 6250 K), which tend to have edge-on inclinations such that ψ ≈ |λ| (see Louden et al. 2021; note that the average difference between ψ and |λ| is just 13.2 • for the 14 systems with measurements for both angles).
Of the 22 sub-Saturns with λ measurements considered in this work, five are in compact multi-planet systems, while the other 17 are isolated.Of the five in compact configurations, none are misaligned, while 9/17 (53%) isolated sub-Saturns are misaligned.We define misalignment as |λ| > 10 • at the 1σ level (i.e., |λ| − σ λ > 10 • ) and λ > 0 • at the 2σ level (i.e., |λ| − 2σ λ > 0 • ), where we adopt the 1σ uncertainties (σ λ ) reported for each λ measurement (see Table 2).In order to perform a fair comparison between the compact multiplanet sample and the isolated sample, we randomly select and classify the alignment of five isolated sub-Saturns.We perform 100,000 iterations of these random draws, counting the number of misaligned systems for each, and fit a Gaussian function over the resulting distribution of N systems per count.The upper-right sub-panel of Figure 2 reveals the result, illustrating that sub-Saturns in compact multi-planet systems are more often aligned than isolated sub-Saturns to a significance level of 2.6σ.
We verified that this significance is robust to the precise misalignment criteria adopted by re-running the above procedure with the following alternative cuts: λ > 15 • and λ > 0 • at the 2σ level (2.6σ significance), (ii) λ > 20 • and λ > 0 • at the 2σ level (2.6σ significance), (iii) λ > 15 • (2.3σ significance), and (iv) λ > 20 • (2.3σ significance).Assuming our original alignment criteria, we additionally vary the Neptune/sub-Saturn mass boundary and re-run our procedure: (i) 15 M ⊕ (3.0σ significance), and (ii) 0.1 M J (1.8σ significance).We note that while the significance appears to degrade with an increasing mass bound, this trend is likely not astrophysical, as the significance is directly proportional to the number of sub-Saturns in compact systems considered, which drops from six assuming a 15 M ⊕ bound to only three assuming a 0.1 M J (≈ 32 M ⊕ ) bound.Finally, we note that the young (∼ 20 Myr), aligned, compact sub-Saturn system AU Mic still hosts a dusty debris disk (Plavchan et al. 2020;Hirano et al. 2020a), meaning that its spin-orbit angle and compactness has likely not yet had the opportunity to be affected by various pathways to trigger post-disk eccentric migration (Wu et al. 2023).Therefore, while this system provides strong support for primordial alignment, it may bias our statistical test; we re-run our original procedure excluding AU Mic (such that there are only four compact sub-Saturn systems) and find that the significance is maintained, but reduced to 2.2σ.
As the spin-orbit angles of compact multi-planet systems should remain largely unaltered following the dispersal of the disk, the low stellar obliquities observed for sub-Saturns and other types of exoplanets in compact systems suggest that most planetary systems are initially formed spin-orbit aligned.Consequently, the dominant mechanism driving misalignment in single-star systems is likely not a universal process that operates indiscriminately across different types of systems, but instead is inherent to those with certain postdisk dynamical histories.More broadly, support for misalignment being acquired well into the post-disk phase rather than near the onset of system formation comes from evidence that i) both young stellar systems (≲ 100 Myr old) and planets still embedded within their debris disks are aligned (Kraus et al. 2020;Plavchan et al. 2020;Hirano et al. 2020a;Martioli et al. 2020;Palle et al. 2020;Albrecht et al. 2022;Johnson et al. 2022), and ii) hot Jupiter systems that are misaligned tend to be older (Hamer & Schlaufman 2022).In combination with the large spin-orbit angles commonly observed for misaligned isolated hot Jupiters and sub-Saturns, these facts are generally consistent with violent eccentric migration as the dominant driver of misalignment (Dawson & Johnson 2018;Wu et al. 2023); identifying the most relevant migration pathway(s), however, is beyond the scope of this work.The sub-panel in the top right displays the resultant distributions of the number of misaligned sub-Saturns found in compact multiplanet systems (dashed line) and the number of misaligned isolated sub-Saturns after performing 100,000 random draws (bars), with statistical significance (2.6σ) indicated.In both panels, orange corresponds to sub-Saturns in compact multi-planet systems while blue corresponds to isolated sub-Saturns.

On the Temperature-Obliquity Relation
There is an interesting discrepancy between the obliquity distributions of isolated hot Jupiters and isolated hot sub-Saturns: hot Jupiters are commonly misaligned around hot stars (T eff ≳ 6250 K) but exclusively aligned around cool stars (i.e., the T eff − λ relation, see Winn et al. 2010a;Albrecht et al. 2012Albrecht et al. , 2022)), while hot sub-Saturns are commonly misaligned around cool stars (and no data exist for those orbiting hot stars, see Figure 2).Realignment via tidal dissipation offers a plausible explanation for this phenomenon; cool stars are susceptible to realignment by massive hot Jupiters, while lower-mass hot sub-Saturns are less likely to realign them on timescales shorter than the system's age (e.g., see Eqn. 2 of Albrecht et al. 2012).Hence, both hot Jupiters and hot sub-Saturns may indeed become misaligned through high-eccentricity migration, but only cool stars with hot Jupiters can be realigned.This tidal realignment mechanism, however, requires finely-tuned parameters in order to completely realign all hot Jupiters around cool stars without leaving any on polar or retrograde orbits (Rogers et al. 2013;Li & Winn 2016), which are not seen in observations.
Alternatively or in addition to this tidal realignment scenario, it is possible that sub-Saturns are subject to misalignment mechanisms that may not operate for Jupiters.For instance, Petrovich et al. (2020) illustrates that secular reso-nance between an outer gas giant companion and an inner sub-Saturn, both embedded in a decaying disk, can drive the sub-Saturn to a polar orbit.The level of obliquity excitation depends on the ratio of the angular momentum between the inner and outer orbits, with nearly polar orbits preferentially excited for inner sub-Saturn-mass planets over Jupiter-mass planets.This resonance sweeping offers a natural explanation for the excess of sub-Saturns observed to be on polar orbits, regardless of stellar effective temperature (Bourrier et al. 2018;Hixenbaugh et al. 2023), as well as a potential mechanism to stunt the runaway growth of sub-Saturns since they are assumed to be embedded within an actively dispersing disk.However, this mechanism predicts similar misalignment rates for compact multi-planet and isolated systems, holding all other factors equal.This would imply that isolated sub-Saturns, which are more often misaligned, host outer companions at a higher rate than sub-Saturns in compact systems, but no statistical studies on these relative rates have yet been completed.
To explain the T eff − λ relation and spin-orbit misalignment in single-star exoplanetary systems more broadly, Hixenbaugh et al. ( 2023) outline a unified model wherein planets gain misalignments if other planets of comparable (or higher) mass form within their disk.Specifically, hotter, more massive stars tend to host more massive protoplanetary disks (Williams & Cieza 2011;Andrews et al. 2013;Andrews 2020), which are more likely to form multiple Jupiter-mass planets (Johnson et al. 2010;Ghezzi et al. 2018;Yang et al. 2020).These Jupiters may interact with each other to produce spin-orbit misalignments after the disk dissipates.By contrast, the limited disk mass around cooler stars is unlikely to produce more than one Jupiter-mass planet, if any (Andrews et al. 2013;Ansdell et al. 2016;Pascucci et al. 2016;Dawson & Johnson 2018;Yang et al. 2020).As a result, these giants, which have no comparable-mass planetary perturbers, should largely retain their primordial alignments throughout the post-disk phase.Nevertheless, despite their less massive disks, cool stars may still form several sub-Saturn planets that undergo an array of dynamical interactions in the postdisk phase that are capable of exciting their stellar obliquities, analogously to multiple Jupiters around hot stars.As presentday compact systems may have avoided violent dynamical interactions, they preserve their primordial spin-orbit angle, allowing them to serve as a window into the early, undisturbed state of planetary systems.

Figure 1 .
Figure1.Keplerian signal-subtracted radial velocities (top panels) and RM model residuals (bottom panels) for TOI-5126 (left) and TOI-5398 (right).In the top panels, the best-fit RM models are shown with red dashed lines, and their uncertainties are indicated by grey shadows, while the planets' spin-orbit configurations are depicted in the bottom left portion of the plot.Note that the NEID RV measurement obtained ≈ 1 hour prior to ingress for TOI 5398 b was taken on April 10, 2022 via 300-second exposure as part of our out-of-transit follow-up for this system.All NEID data used in this work are available here.

Figure 2 .
Figure2.Sky-projected stellar obliquity λ versus stellar effective temperature T eff for all sub-Saturns considered in this work, with 1σ error bars shown.The sub-panel in the top right displays the resultant distributions of the number of misaligned sub-Saturns found in compact multiplanet systems (dashed line) and the number of misaligned isolated sub-Saturns after performing 100,000 random draws (bars), with statistical significance (2.6σ) indicated.In both panels, orange corresponds to sub-Saturns in compact multi-planet systems while blue corresponds to isolated sub-Saturns.

Table 2 .
Sub-Saturns Around Single Stars With a Stellar Obliquity Measurement