Rotation Periods, Inclinations, and Obliquities of Cool Stars Hosting Directly Imaged Substellar Companions: Spin-Orbit Misalignments are Common

The orientation between a star's spin axis and a planet's orbital plane provides valuable information about the system's formation and dynamical history. For non-transiting planets at wide separations, true stellar obliquities are challenging to measure, but lower limits on spin-orbit orientations can be determined from the difference between the inclination of the star's rotational axis and the companion's orbital plane ($\Delta i$). We present results of a uniform analysis of rotation periods, stellar inclinations, and obliquities of cool stars (SpT $\gtrsim$ F5) hosting directly imaged planets and brown dwarf companions. As part of this effort, we have acquired new $v \sin i_*$ values for 22 host stars with the high-resolution Tull spectrograph at the Harlan J. Smith telescope. Altogether our sample contains 62 host stars with rotation periods, most of which are newly measured using light curves from the Transiting Exoplanet Survey Satellite. Among these, 53 stars have inclinations determined from projected rotational and equatorial velocities, and 21 stars predominantly hosting brown dwarfs have constraints on $\Delta i$. Eleven of these (52$^{+10}_{-11}$% of the sample) are likely misaligned, while the remaining ten host stars are consistent with spin-orbit alignment. As an ensemble, the minimum obliquity distribution between 10-250 AU is more consistent with a mixture of isotropic and aligned systems than either extreme scenario alone--pointing to direct cloud collapse, formation within disks bearing primordial alignments and misalignments, or architectures processed by dynamical evolution. This contrasts with stars hosting directly imaged planets, which show a preference for low obliquities. These results reinforce an emerging distinction between the orbits of long-period brown dwarfs and giant planets in terms of their stellar obliquities and orbital eccentricities.


Introduction
Measurements of stellar obliquity-the orientation between a star's spin axis and a planet's orbital angular momentum vector-provide valuable clues about how planets form, migrate, and dynamically interact (Winn & Fabrycky 2015;Albrecht et al. 2022). In the absence of external perturbers, the axisymmetric collapse of a molecular cloud core is expected to result in mutual alignment between the stellar spin and the protoplanetary disk's axis of rotation, which would then be inherited by any planets that form in the disk.
Many processes can disrupt this alignment and produce a wide range of stellar and planetary obliquities during and after the phase of planet formation. These mechanisms can act on the star, the planet, or the protoplanetary disk and include the influence of passing stars (Heller 1993;Batygin et al. 2020); secular torques from wide stellar and substellar companions (Fabrycky & Tremaine 2007;Batygin 2012;Anderson et al. 2016;Lai et al. 2018); and gravitational scattering with other planets (Chatterjee et al. 2008;. Primordial misalignments of protoplanetary disks (both "intact" and "broken") appear to be common, indicating that some planets probably form with substantial spin-orbit angles (e.g., Bate et al. 2010;Huber et al. 2013;Marino et al. 2015;Ansdell et al. 2020;Epstein-Martin et al. 2022). Even the Sun is misaligned by about 6°with respect to the solar system's invariable plane, but the origin of this enigmatic offset continues to be debated (e.g., Adams 2010; Bailey et al. 2016;Spalding 2019).
Large transiting planets have provided a wealth of information about stellar obliquities, largely through individual and collective studies of Rossiter-McLaughlin measurements (McLaughlin 1924;Rossiter 1924;Queloz et al. 2000;Gaudi & Winn 2007). It is now clear that nonzero obliquities are common among field stars hosting short-period planets Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. (Fabrycky & Winn 2009;Schlaufman 2010;Muñoz & Perets 2018;Louden et al. 2021), and especially those with temperatures above the Kraft break at ≈6250 K (Winn et al. 2010;Schlaufman 2010). As an ensemble, however, most compact multiplanet systems tend to have lower stellar obliquities (Albrecht et al. 2013;Morton & Winn 2014;Winn et al. 2017). A complex picture is emerging in which several mechanisms are probably responsible for tilting and in some contexts realigning the orbits of close-in planets over many different timescales (e.g., Morton & Johnson 2011;Albrecht et al. 2012;Dawson 2014).
Less, however, is known about stellar obliquities for systems hosting long-period companions. At wide separations, stellar obliquities can be used to provide clues about the formation and orbital evolution of companions amenable to direct imaging (Bowler 2016). In particular, planets and brown dwarfs have been discovered at separations spanning tens to thousands of au and may originate from a combination of outward gravitational scattering (Boss 2006;Veras et al. 2009;Bailey & Fabrycky 2019), in situ formation in a disk via pebble accretion or disk instability (Durisen et al. 2007;Lambrechts & Johansen 2012), dynamical capture (Perets & Kouwenhoven 2012), or direct collapse from a fragmenting molecular cloud core (Bate 2012). There is growing evidence that most directly imaged planets within about one hundred astronomical unit are formed in a disk based on their eccentricity distributions , mass functions (Wagner et al. 2019), and demographic trends (Nielsen et al. 2019;Vigan et al. 2021). Spin-orbit angles (between the stellar spin axis and companion orbital plane) and spin-spin angles (between the spin axes of the host star and companion) can provide another unique perspective on the origin and orbital evolution of this population.
Measuring stellar obliquities is challenging at wide separations. The Rossiter-McLaughlin effect becomes less practical as a tool to study projected spin-orbit orientations because the geometric probability of a transit is lower, transit events are less frequent, and transit durations become longer; this explains why few systems with orbital periods beyond a few tens of days have had stellar obliquity measurements . Our strategy in this study is to use stellar rotation periods, projected rotational velocities, and stellar radii to infer the spin orientation of stars hosting imaged substellar companions. The goal is to uniformly establish stellar inclinations for the full sample of host stars, then determine stellar obliquities for the subset of these with measurements of the orbital inclination from orbit monitoring. As a byproduct of this analysis, we also present rotation periods for 62 host stars, which can be used to refine the ages of these systems using gyrochronology.
This paper is organized as follows. In Section 2 we summarize the geometry and overall framework to constrain inclinations and obliquities for host stars of imaged companions. A description of our sample, new high-resolution spectroscopy, and Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) light curves are detailed in Section 3. Section 4 includes our * v i sin measurements, periodicity analysis, stellar line-ofsight inclinations, and obliquity constraints. The interpretation of our results and a discussion of sample biases can be found in Section 5. Our conclusions are summarized in Section 6.

Stellar Obliquity Geometry for Directly Imaged Companions
The problem of observationally measuring the angle between the stellar spin axis and a companion's orbital plane, ψ, has a long history in the context of stellar binaries and the coplanarity of disks (e.g., Hale 1994;Watson et al. 2011). For clarity we briefly review the basic geometric setup of this problem as it applies to directly imaged companions. Figure 1 shows the relevant orbital elements of a visual binary companion (left diagram) and the orientation of the spin axis of its host star (middle diagram). The observer is looking down along the z-axis toward the x-y sky plane. The orbital plane of an imaged planet is characterized by its inclination i o with respect to the sky plane, which is equal to the angle between the z-axis and the orbital angular momentum vector Figure 1. Orbital and rotational geometry relevant for obliquity constraints of stars hosting imaged companions. The observer is looking down from the z-axis. The xy plane is the plane of the sky, and the star is centered on the origin. The x-axis points north and the y-axis points east. Left: the inclination of the companion's orbital plane, i o , intersects the sky along the line of nodes, which is oriented by Ω o from true north to the ascending node. In this configuration the orbital angular momentum vector ¾  L o points in the hemisphere facing the observer, with ¾  L o pointing toward the observer along the z-axis when i o = 0°. Middle: the star's orientation in space determines the inclination of its equatorial plane, i * , and the position angle of its line of nodes. The sense of the stellar spin sets its rotational angular momentum vector ¾  * L and the longitude of ascending node Ω * . Right: if only the inclinations i o and i * are measured, but not their true (or relative) orientations in the plane of the sky (Ω o and Ω * ), then each angular momentum vector can fall anywhere along nested cones opening toward or away from the observer. If Ω o is known (point A 1 , for example), but not Ω * , then the obliquity can be as small as |i o -i * | (at B 1 ), can reach i o + i * for prograde orbits (B 2 ), or can be as large as π − |i o − i * | (B 3 ) for retrograde orbits.
¾  L o , and the longitude of ascending node Ω o , defined from celestial north (here the x-axis) increasing eastward (toward the y-axis) to the line of nodes. i o spans 0-π and Ω o spans 0-2π. In general, orbit monitoring with relative astrometry results in an ambiguity between Ω o and Ω o + 180°, which can be distinguished with radial velocities (RVs) of the companion (e.g., Snellen & Brown 2018;Ruffio et al. 2019Ruffio et al. , 2021 or its host star (e.g., Crepp et al. 2012;Bowler et al. 2018).
The setup is similar for the spinning star: its equatorial plane is inclined by i * from the sky plane, equal to the angle between the observer's line of sight and the star's rotational angular momentum vector ¾  * L . The direction in which ¾  * L points is determined by its orientation on the sky, Ω * . With few exceptions for direct measurements of stellar oblateness (e.g., Belle 2012) and resolved rotation (Kraus et al. 2020), Ω * is generally unknown.
The relationship between i o , Ω o , i * , and Ω * can be derived by an observer-to-orbit coordinate system transformation (Fabrycky & Winn 2009) or via the spherical law of cosines (Glebocki & Stawikowski 1997): For transiting planets, i o is known (≈90°) and the projected spin-orbit orientation ΔΩ = Ω o − Ω * (often denoted as λ) can be measured with the Rossiter-McLaughlin effect (Holt 1893;McLaughlin 1924;Rossiter 1924). In this case, i * in general is unknown, and y cos ≈ ( ) l * i sin cos , so λ is a lower limit on the true obliquity.
For directly imaged planets, i o can be measured through orbit monitoring. In some cases the inclination of the host star can be determined through asteroseismology (e.g., Zwintz et al. 2019) or the projected rotational velocity technique, which relies on knowing * v i sin , the rotation period P rot , and the stellar radius R * (Shajn & Struve 1929). This latter approach is the focus of this study and is possible for spotted stars whose rotational modulations manifest as periodic light-curve variations.
There is an important distinction between simple coplanarity of the orbital and stellar equatorial planes, versus genuine alignment of the orbital and rotational angular momentum vectors between the companion and the host star-the primary angle of interest in this study. When considering only coplanarity (which could include spin-orbit "anti-alignment," or retrograde orbits) and both i o and i * are known (but not ΔΩ), the sum and absolute difference between the orbital and stellar inclinations provide boundaries on the true spin-orbit alignment: On the other hand, if one is concerned with angular momentum alignment in the system, and only i o and i * are known, there exists a broader range of possibilities for ψ because the measurement of i * does not contain information about the sense of the spin. Depending on the direction of rotation, the spin angular momentum vector may be pointing toward the observer or away from the observer into the plane of the sky. In this case the constraint on ψ becomes: Regardless, if i o and i * differ then ψ is nonzero and the system's orbital and rotational planes are misaligned by at least |i o − i * |. However, if i o and i * are identical then the system is consistent with being aligned but the true obliquity may be much larger. This is illustrated in the rightmost diagram in Figure 1. If only i o and i * are known, the projection of the orbital and rotational angular momentum vectors will trace out nested cones with a degeneracy between the cones opening toward and away from the observer. If i o is fixed at point A 1 , for example, and a stellar inclination i * is measured, the minimum value of ψ will occur when Ω * = Ω o at point B 1 . If Ω * = Ω o + 180°(B 2 ), ψ can reach i o + i * . But if ¾  * L points away from the observer (residing on the lower cone opening in the −z direction), ΔΩ = 180°(B 3 ), and ψ can reach a value as high as 180°− |i o − i * |. For instance, if i o = i * = 30°, the orbital and rotational planes can be perfectly coincident (ψ = 0°) or misaligned by up to 60°, and the spin-orbit angle can be misaligned by up to 120°. This study focuses on measurements of the minimum misalignment, Δi = |i o − i * |, for stars hosting directly imaged substellar companions.

Sample of Substellar Hosts
Our sample originates from a compilation of 177 stars hosting imaged substellar companions spanning a wide range of orbital separations. The list includes discoveries based on high-contrast imaging as well as wide common-proper motion companions identified from seeing-limited surveys. The compilation is assembled predominantly from the more focused lists in Deacon et al. (2014), Bowler (2016), and Bowler et al. (2020a), as well as additional systems identified more recently. The properties of the original parent sample are broad. Host star masses span the hydrogen burning limit up to supernova progenitors (Squicciarini et al. 2022); companions range in mass from 2 to about 75 M Jup and orbital separations from a few au to thousands of astronomical unit ( Figure 2).
From this list we identify stars with rotation periods and * v i sin measurements-either new in this work or previously published. When considering rotation periods, we restrict our analysis to spectral types of F5 or later to reduce confusion between rotation periods and stellar pulsations (e.g., ). 12 This is especially true of the intermediate-mass γ Dor pulsators, which span late-A to mid-F spectral types and have characteristic g-mode oscillation frequencies of a few hours to a few days (Kaye et al. 1999;Aerts 2021). Light curves and periodograms of stars in the remaining sample are further visually inspected for signs of pulsations.
Two exceptions are made to this spectral type cut. HIP 21152 is an F5 member of the Hyades with a recently discovered brown dwarf companion (Bonavita et al. 2022;Franson et al. 2023;Kuzuhara et al. 2022). Its light-curve periodogram shows significant peaks that are not integer harmonics of its strongest signal, and are therefore not likely to be caused by rotational modulation. We therefore remove this star from our sample. 12 Note that this is not an explicit cut on mass because spectral type generally evolves during the pre-and post-main-sequence phases with changing T eff . However, for mid-F stars the difference is modest. On the main sequence, F5 corresponds to an effective temperature of ≈6510 K and a mass of about 1.4 M e (Drilling & Landolt 2000;Pecaut & Mamajek 2013). Similarly, for the young (≈18 Myr) equal-mass F5 binary AK Sco, Czekala et al. (2015) found a total dynamical mass of 2.49 ± 0.10 M e , implying individual masses of about 1.25 M e . 51 Eri is a young F0 member of the β Pic moving group with an imaged planet at 11 au (Macintosh et al. 2015).  found the host to be a γ Dor pulsator, casting doubt on previous rotation periods determined from light curves. However, they also found a plausible core rotation rate of -+ 0.9 0.1 0.3 days, so we include this star in our sample but recognize that detailed asteroseismic analysis using a longer time series is needed to more reliably constrain the rotation properties of this massive star.
Finally, all light curves with low-amplitude modulations are carefully examined on a case-by-case basis for signs that the variations might be instrumental rather than astrophysical in nature. Optical artifacts from bright sources and scattering from earthshine can result in low-level changes that may mimic real signals. Our criteria for retaining light curves are that the signals must show clear and consistent periodic brightness changes that are evident from visual inspection, rotation periods must be shorter than half of the TESS sector baseline (13 days for a single sector), and amplitudes must be substantially higher than the expected sensitivity floor for the target star brightness.
After implementing these cuts, our final sample for this study comprises 62 host stars with rotation periods, 53 of which also have projected rotational velocities and stellar radius determinations. These measurements are detailed in the following subsections and are summarized in Tables 2 and Figure 2. A Gaia color-magnitude diagram of host stars in our final sample is shown in Figure 3.

High-resolution Spectroscopy with the Tull Spectrograph
High-resolution optical spectra for 22 targets from our broader sample were obtained with the Tull Coudé Spectrograph (Tull et al. 1995) at McDonald Observatory's 2.7 m Harlan J. Smith telescope between 2019 March and 2022 February. The 1 2 slit and E2 grating were used with the TS23 setup for all observations, which resulted in a resolving power of R = λ/δ λ ≈ 60,000. Fifty-six échelle orders were simultaneously captured with the TK3 Tektronix CCD from 3870-10500 Å in (largely nonoverlapping) segments ranging from 60-200 Å. On most nights we also observed several bright, slowly rotating RV and * v i sin standards from Chubak et al. (2012) and Soubiran et al. (2018) spanning F through M spectral types. A summary of the observations can be found in Table 1.
Two-dimensional traces for the science observations are defined using a spectral flat field calibration frame taken on the same night. After bias subtraction and correction for known bad pixels, each curved spectral order is resampled to a linear two-dimensional trace using subpixel interpolation. Night sky lines are subtracted from each order using dispersed sky regions near the ends of the projected slit. Orders are optimally extracted following the method of Horne (1986), and fifthorder polynomial fits to the extracted flat spectra are used to remove the blaze function.
A wavelength solution is derived for each of the 56 orders using a ThAr emission line spectrum obtained on the same night and with the same setup as the science observations. A ThAr line list comprising a total of 9145 emission lines with documented relative line intensities was assembled spanning 3785-10507 Å from Lovis & Pepe (2007; for λ < 6915 Å) and Murphy et al. (2007; for λ > 6915 Å). For each order, a synthetic emission line spectrum was generated at the resolving power of the Tull spectrograph. To automate the process of mapping pixel values to wavelengths, we iteratively solved for the coefficients of a third-order polynomial model by maximizing the cross-correlation function (CCF) between the Figure 2. Overview of the host star spectral types and rotation periods (left), projected rotational velocities (middle), and companion spectral types and separations (right). Systems with stellar obliquity constraints are plotted in dark blue. Companions with obliquity constraints in this study reside within ≈250 au where orbital motion is more readily detectable. . Gaia G BP -G RP color-magnitude diagram of stars hosting imaged substellar companions in this study. With the exception of 51 Eri, the bluest star in the figure, we focus on cool stars of F5 and later to avoid confusion between stellar rotation and pulsations in our light-curve analysis. Host stars are color-coded by rotation period and range from 0.08-44 days. Long rotation periods are either determined from multiple TESS sectors or are adopted from previous measurements in the literature. observed spectral order and a "true" emission spectrum using the AMOEBA algorithm (Nelder & Mead 1965). With reasonable initial estimates for the coefficients, this approach performed well, and the solution was visually verified for each order on each night.

TESS Light Curves
We used the lightkurve (Lightkurve Collaboration et al. 2018) package to download the 2 minute cadence Science Processing Operations Center (SPOC) Pre-search Data Conditioning Simple Aperture Photometry light curve Stumpe et al. 2012Stumpe et al. , 2014Jenkins et al. 2016) from the Mikulski Archive for Space Telescopes (MAST). 13 Two targets (GJ 3305 and SDSS J130432.93+090713.7) were analyzed using their 30 minute full-frame images (FFIs), which were reduced by either the TESS-SPOC (Caldwell et al. 2020b) or the MIT Quick-Look Pipeline ) and manually retrieved from the MAST online portal. 14 All photometric measurements listed as NaN are removed. Individual TESS sector light curves are normalized by dividing both the flux and flux uncertainties by the median flux value of that sector. Each complete light curve is then compiled by  (Savitzky & Golay 1964) and only selecting data points in the original light curve that lie within one standard deviation of the flattened light curve. Note that because of flares and outlier photometric points, this threshold is typically larger than what it would be for pure Gaussian noise. Our period measurements (described below) are in good agreement whether or not we include this outlier rejection step. For each processed light curve used in this study, we produce generalized Lomb-Scargle periodograms (Zechmeister & Kürster 2009) over the frequency range 0.0005-100.0 day −1 (0.01-2000 days) to search for any periodic modulation that could be attributed to rotation. A Gaussian is fit to the highest periodogram peak and the resulting mean and standard deviation are adopted for the period and its uncertainty, respectively. For stars with multiple sectors separated by data gaps, aliasing due to the spectral window function causes a "fringe pattern" with a broad envelope for all periodogram peaks. In such cases, the Gaussian is fit to the envelope structure in order to reflect that spread. Each resulting light curve is visually inspected to ensure significant periodicities reflected astrophysical variations instead of spacecraft-related systematic oscillations, for example from optical artifacts or earthshine. These can manifest as gradual low-amplitude variations, sharp discontinuities from spacecraft motion, and scattered light features that rise and fall in sync with the 13.7 day orbit. Period measurements are considered reliable if they show a single strong peak, a phase-folded light curve with consistent structure among phased curves, and a rotational modulation amplitude that is large compared to the limiting sensitivity of TESS for that target brightness. Finally, we limit adopted measurements to periods 13 days, which corresponds to half of the TESS sector baseline; instrumental systematics make it challenging to recover lower-frequency signals in TESS light curves (e.g., Avallone et al. 2022).
Altogether 46 stars have light curves that appear to reflect true rotational modulation. They exhibit a wide range of characteristic behaviors-amplitudes with 20%-level variations (PDS 70); consistent, symmetric, sinusoid-like modulations (2MASS J02155892-0929121, CD-35 2722, and G 196-3); signs of dramatic and rapid starspot evolution (HD 130948, HD 16270, HD 49197, and HD 97334); double-peaked structure (2MASS J01225093-2439505 and PZ Tel); and beating patterns reflecting the superposition of two frequencies, most likely reflecting binarity or differential rotation (TWA 5 Aab). Full light curves, periodograms, and phased light curves are shown in Figures 20-27. Rotation periods and TESS sectors used in this analysis are listed in Table 2. LP 261-75 reflects an interesting example of an M dwarf whose rotation period would have been interpreted incorrectly if only one TESS sector had been available. The periodogram of its full light curve from Sectors 21 and 48 shows two comparably strong peaks at 1.11 and 2.23 days ( Figure 25); one signal is clearly a harmonic of the true rotation period. This confusion is evident in reported photometric rotation period measurements of this star in the literature: Newton et al. (2016) and Irwin et al. (2018) found values of 2.219 days and 2.22 days, respectively, while Canto Martins et al. (2020) listed a value of 1.105 ± 0.027 days. We therefore separately analyzed each individual sector light curve. The Sector 21 light curve shows a strong peak at 1.11 days with smooth modulations and consistent amplitudes. Seven-hundred twelve days later marks the beginning of the next set of observations in Sector 48. In this sector, the regular modulations occur with a periodicity of 2.23 days.
We interpret this longer modulation as the true rotation period. This "period-doubling" behavior is consistent with the following scenario: during the Sector 21 observations, two groups of starspots with comparable surface covering fractions may have been located on opposite sides of the spinning star with a longitudinal phase difference of ∼180°. Typically this setup with two large groups of spots would produce a light curve with a double-peaked structure having two amplitudes, but here the similar spot sizes and 180°phase difference seem to have conspired to mimic a faster period in Sector 21, which only later could be differentiated from the true value when one group of spots evolved or disappeared. This phenomenon is different from typical rotation period confusion from photometric light curves, which are usually caused by sampling and windowing effects from insufficient cadence relative to the rotation period of the star. This is not the case for these fine 2 minute observations from TESS. Instead, LP 261-75 illustrates a (probably rare) instance of bias that can arise if only a single sector is available, which applies to 13 of the 62 host stars with rotation periods in our sample. Note that the deep 1.882 day transit events from LP 261-75 C (Irwin et al. 2018) were removed from the light curves for this analysis. These small regular gaps do not impact the results.

New Projected Rotational and Radial Velocities
Stars that we observed with the Tull Spectrograph range in spectral type from B9 to M5. Projected rotational velocities are determined by either applying a rotational broadening kernel to spectra of slowly rotating RV standards, or by broadening an appropriate model spectrum from the PHOENIX-ACES grid (Husser et al. 2013). Orders with low signal-to-noise ratio (S/ N), strong telluric lines, or few absorption lines (for early-type stars) were avoided. This typically resulted in 20-40 good orders. When using the empirical standards as a reference template, we first cross correlate every standard spectrum observed on the same night with each order of the science spectrum. The resulting individual CCFs for a particular standard are summed, and the reference spectrum that results in the highest CCF value is adopted for the broadening analysis. Each order of the template spectrum is then sequentially broadened through a fine grid of * v i sin values spanning 0-200 km s −1 with a rotational kernel following Gray (2005). A new CCF is computed for each * v i sin value, and the broadening kernel resulting in the highest CCF peak for that order is added in quadrature with the measured * v i sin value of the standard star. This is typically obtained from the compilation of projected rotational velocities by Glebocki & Gnacinski (2005); most of the standards we used have * v i sin values below 3 km s −1 for G, K, and M dwarfs and below 10 km s −1 for F stars. The final adopted projected rotational velocity is then derived from the robust mean and standard deviation of the * v i sin measurements for each individual order following the procedure in Beers et al. (1990). Uncertainties generally scale with the * v i sin value, but most are measured at the s *   Notes. a Total period uncertainty including periodogram measurement uncertainty and a term to account for potential differential rotation. We adopt a solar absolute shear of 0.072 rad day −1 for all stars except 51 Eri. See Section 4.2 for details. b Stellar radius estimates from Stassun et al. (2019) have been inflated by 7% based on a comparison in that study between the original estimates and those from asteroseismology. c GJ 3305 AB is a wide binary to the exoplanet host 51 Eri. We include the light-curve analysis and inclination constraint in this study for completeness, but do not count this system itself as a substellar host. This procedure also results in an instantaneous RV relative to the RV standard. A barycentric correction is applied following the approach in Stumpff (1980), then an absolute RV is computed using the RV of the standard. The total uncertainty is determined from the quadrature sum of our measured relative RV and the published uncertainty for the standard star. Although these measurements are not the primary focus of this study, we report RVs alongside * v i sin values in Table 1.
On nights when an appropriate slowly rotating standard was not observed, or for early-type stars for which fewer standards are available, we used solar-metallicity synthetic spectra from the PHOENIX-ACES model atmosphere grid to determine projected rotational velocities. Model effective temperatures were chosen to the nearest 100 K (for T eff < 7000 K) or 200 K (for T eff 7000 K) using the SpT-T eff scale from Pecaut & Mamajek (2013). The models were first smoothed with a Gaussian kernal to the resolving power of the Tull spectrograph, trimmed to the wavelength range of the individual orders, and resampled onto the same wavelength grid as the science spectra. The same procedure for assessing rotational broadening is carried out for the models as for the standard star observations described above. Fewer orders were generally used for the B, A, and F stars as many orders for these earlytype stars lacked lines altogether.
* v i sin values are determined separately for each order and then a robust mean and standard deviation are computed. Model effective temperatures are listed with the * v i sin results in Table 1. Examples of spectra and broadened templates (both empirical and from model spectra) for a single order are shown in Figure 4. We compiled previously determined values for targets in Table 1 to assess how our measurements compare with those in the literature (see Appendix D and Section 4.3). Results are shown in Figure 5. Above about 6 km s −1 , the agreement between our measurement and previous values is generally good; the Pearson correlation coefficient is 0.99, and the rms of the relative residuals is 20%. Below 6 km s −1 , our measurements using the order-by-order CCF approach generally overestimate reported values. This is likely because we are not carrying out a detailed line profile analysis of high-S/N spectra focused on select individual lines. We are limited by previous determinations of broadening in standard stars, which result in a lower limit on our achieved value, and we do not take into account any additional broadening effects such as micro-and macroturbulence that could contribute to the shape and strength of absorption lines. We therefore caution that low * v i sin values from our analysis may be somewhat overestimated. However, only eight of our 45 measurements are below 6 km s −1 , and among stars that also have rotation periods (which enable an inclination analysis), all but one (1RXS J034231.8+121622) have previous * v i sin measurements in the literature. These previous measurements are taken into account in the final adopted * v i sin value (Appendix D and Table 2).

Interpreting Periodic Modulations as Rotation Periods
Starspots are the main source of large-scale periodic photometric variability in intermediate-and low-mass stars. In the Sun, differential rotation gives rise to surface rotation periods of about 25 days at the equator and 34 days at the poles. The locations of sunspots vary during the solar activity cycle but are generally confined to latitudes of about ±30°. However, for other stars with different convection zone depths, rotation rates, masses, compositions, and evolutionary states, the behaviors of differential rotation, starspot locations, and spot filling factors are expected to fundamentally differ from that of the Sun (Schussler et al. 1996;Barnes et al. 2005;Berdyugina 2005). If spots are positioned at nonequatorial latitudes, differential rotation will result in a rotationally modulated signal that deviates from the equatorial rotation period. This can bias measurements of the stellar equatorial velocity and inclination using the projected rotational velocity, Absolute shear ΔΩ is a common metric to quantify differential rotation. It represents the difference in angular frequency at a star's equator and its pole: Equator Pole min max Reinhold et al. (2013) and Reinhold & Gizon (2015) analyzed thousands of (typically old) stars from the Kepler mission and found that absolute shear as traced by multiple periodogram peaks can span a wide range of values for a given rotation period-a proxy for age for Sun-like stars-and effective temperature, but does not vary strongly as a function of these parameters, at least below about 6200 K (≈F8). The typical range of ΔΩ values spans 0.01-0.1 rad day −1 . At higher effective temperatures, the absolute shear increases to values greater than 1 rad day −1 . For comparison, the Sun's absolute shear is 0.07 rad day −1 . Because differential rotation is likely to play an important role in interpreting the periodicity of light curves from stars in our sample, we inflate the rotation period uncertainties computed from periodograms to more accurately reflect the unknown latitude being tracked by dominant starspot groups. We define an error term associated with the star's shear, σ P,shear , that is designed to conservatively reflect half the difference between the maximum (polar) and minimum (equatorial) rotation periods: Solving for P min and inserting it into Equation (4) allows us to relate σ P,shear , P max (for which we adopt our measured photometric rotation period), and ΔΩ: The Sun, for instance, would have σ P,shear of 2.5 days-about 10% of its actual equatorial velocity. Assuming the same absolute shear, a younger Sun with P rot = 10 days would have σ P,shear = 0.5 day, a 5% relative uncertainty. To estimate σ P,shear for stars in our sample, we assume a constant solar-like absolute shear of ΔΩ = 0.07 rad day −1 .
Our final adopted rotation period uncertainty is the quadrature sum of the uncertainty from the periodogram analysis and from potential differential rotation: The median rotation period precision, σ P,tot /P rot , is 2% with a range of 0.2%-20% across our entire sample of 64 stars.

Compilation of Projected Rotational Velocities
* v i sin values are compiled from our own measurements and those in the literature. There is a wide range in the quality of published projected rotational velocities based on the type of observations (for example, medium versus high spectral resolution) and the approach to the measurement itself (such as detailed treatment of broadening mechanisms). Some studies report uncertainties, but many do not. Furthermore, many measurements of * v i sin for the same star are formally inconsistent.
Our strategy for this study is to incorporate as much information as possible while also balancing the reliability of various measurements made over the past several decades. Robust uncertainties are especially important to ensure the accuracy of our stellar inclination posterior distributions. If we obtained more than one * v i sin measurement from our own Tull spectrograph observations, then these are combined into a single value through a weighted mean and weighted standard deviation. If more than one value was identified in the literature, and uncertainties are reported, then these are treated as separate independent measurements. In cases where measurement uncertainties are not reported, we compute their mean and standard deviation and treat this as an additional measurement. All of these are then combined into a single (presumably more accurate) * v i sin value, which we adopt for this study using weighted mean and standard deviation. This relies on reported uncertainties having been reliably estimated so that they avoid unjustifiably driving the weighted mean value to that measurement. We therefore also set an upper limit of 1% on the relative precision of s * v i sin / * v i sin for any single measurement. Details can be found in Appendix D, and individual measurements and adopted * v i sin values are listed in Table 4.

Stellar Radii
Most stellar radius estimates have been adopted from the Revised TESS Input Catalog (TIC; Stassun et al. 2019) to ensure they are determined in a self-consistent fashion. These values are based on the Stefan-Boltzmann relation and incorporate Gaia-based distances, extinction-corrected Gband magnitudes, and G-band bolometric corrections. Stassun et al. (2019) found that the resulting radii were in good agreement (typically within 7%) with values for the same stars measured through asteroseismology by Huber et al. (2017). We therefore inflate uncertainties from the TIC catalog by adding 7% errors in quadrature with the quoted uncertainties to reflect a more realistic spread in individual radius estimates. TIC radii are adopted for 47 of the 64 stars in our full sample.
For the remaining stars, radii are either taken from other literature sources or individually determined in this study. This is the case when TIC radii are not listed, or if we suspect binarity may be severely impacting the inferred values. When determining radii in this study, we make use of the Stefan-Boltzmann relation, and adopt a solar bolometric magnitude of M bol,e = 4.755 mag and an effective temperature of T eff,e = 5772 K (Mamajek 2012;Pecaut & Mamajek 2013).
To decompose the apparent magnitudes of visual binaries into their individual constituent brightnesses, we make use of the contrast Δm in magnitudes and the integrated-light apparent magnitude m. The apparent magnitude of the primary m A in each system is then and the apparent magnitude of the secondary is m B = Δm + m A . Details for individual systems can be found in Appendix B, and radii are listed in Table 2.

Stellar Inclinations
Our approach to constrain i * relies on the Bayesian probabilistic framework developed by Masuda & Winn (2020), which takes into account the correlation between the projected and equatorial rotational velocities * v i sin and v eq . Assuming normally distributed measurement errors for: * v i sin , P rot , and R * ; sufficiently precise measurements of the rotation period (at the 20% level); and uniform priors for all three parameters, we show in Appendix A that the posterior distribution of i * can be expressed as: The resulting line-of-sight inclination posteriors for all 53 host stars are shown in Figures 6-8, and summary statistics are listed in Table 3. Constraints on the stellar inclination vary dramatically from star to star depending on the precision of the input rotation period, projected rotational velocity, and stellar radius. In some cases when these constraints are poor, the data add very little information, and the posterior on i * resembles that for an isotropic prior, * i sin . Most inclination uncertainties span ≈10°-40°, but in some cases the posterior is very well constrained to within 1°-2°. As expected from random orientations, the majority of inclination angles peak at high values of i * > 45°, with many consistent with equator-on orientations of i * = 90°. BD+21 55 has the most pole-on orientation with =  -+ * i 3.8 1.6 1.2 . For cases where the host star is a binary, it is assumed that the measurements of P rot and * v i sin are for the brighter component, and therefore that the inclination usually reflects the more massive member of the system. In these cases, and especially when the mass ratio is near unity, these constraints should be treated with caution because it is possible that the companion could impact the input values of period, radius, or projected rotational velocity.
The inferred equatorial velocity should always be higher than the projected rotational velocity for any line-of-sight inclinations that depart from 90°. In several instances in Table 2, the * v i sin value is larger than v eq , but for most of these systems, they are consistent to within 1σ-2σ. However, there are a few notable exceptions including 2MASS J01225093-2439505 (v eq = 12.4 ± 1.0 km s −1 ; * v i sin = 18.1 ± 0.5 km s −1 ) and Gl 229 (v eq = 1.0 ± 0.2 km s −1 ; * . This discrepancy may be due to Note that the direction of the angular momentum vector in the sky plane is usually unknown, as is its orientation toward or away from the observer. For comparison, an isotropic prior distribution is shown in Figure 8. overestimated rotational broadening measurements, underestimated radius estimates, or periods that are overestimated. Longer-than-expected rotation periods could originate from stars with solar-like dynamos that have both strong differential rotation and spots located at nonequatorial latitudes such that the observed modulations are not tracking equatorial rotation regions. 2MASS J01225093-2439505 and Gl 229 are individually discussed in more detail in Appendix B. For this work, we treat the posterior distributions of i * at face value when because, when this occurs, the Bayesian framework we make use of yields qualitatively similar results to = * v v i sin ; eq both imply an equator-on orientation with a posterior "pressed up" against i * = 90°.

Orbital Inclinations
Probability distributions for orbital inclinations are assembled from the literature for companions with orbit fits     to various combinations of relative astrometry, absolute astrometry, and RVs. This limits the available sample of obliquity constraints to 21 companions at separations where orbital motion is apparent on timescales of years to a few decades. For our sample, the range of separations spans 10-250 au, with most companions residing between 10 and 100 au. 15 When possible, actual orbital inclination posterior distributions are used in our analysis, but in some instances we make use of our own approximations (for example, assumptions of normality) based on the reported summary statistics such as the mean and standard deviation. Orbital inclination posteriors for eight systems are taken from the uniform analysis in Bowler et al. (2020a), which makes use of the orbitize! orbit-fitting package (Blunt et al. 2017(Blunt et al. , 2020 Table 3.

Minimum Stellar Obliquities
Figures 9 and 10 show the posterior distributions of i * and i o visualized in polar coordinates. Here the observer is looking down along the z-axis, and the x-y sky plane is located at an inclination of 90°. Because the orientations of stellar inclinations in the sky plane are unconstrained, as are the senses of their spin directions (clockwise or counterclockwise from the observer's perspective), the stellar rotational angular momentum vectors may point toward the observer or away into the plane of the sky. When assessing spin-orbit alignmentmutual agreement of orbital and angular angular momentum vectors-the stellar inclination distribution is therefore mirrored about i = 90°. In these polar plots the results often look like boomerang throwing sticks, so we refer to these figures as "boomerang diagrams." The consistency between i * and i o is quantified by sampling from these distributions and evaluating whether their absolute 19.4-35.5 Notes. a Maximum a posteriori (MAP) probability. b Systems are classified as misaligned ("Yes") if the probability that Δi values are greater than 10°is 95% and the MAP value of the Δi posteriors is greater than 10°. If P(Δi > 10°) 80% and the Δi MAP value is >10°, then the system is classified as being "Likely" misaligned. If neither is satisfied, the system is consistent with spin-orbit alignment ("No Evidence"). c Our normal approximation to reported i o . d The stellar inclination of κ And is taken from the solar-metallicity model fits to oblateness measurements in Jones et al. (2016). The slightly asymmetric reported uncertainties are averaged to 4°for these purposes. (This table is available in machine-readable form.) 15 Note that in addition to 20 systems from our main sample with i * constraints determined using rotation periods, radii, and projected rotational velocities, we also include in our obliquity analysis κ And-a rapidly rotating B9 star hosting a substellar companion at ≈55 au ( 15°d epending on whether the stellar inclination is ≈30°and the polar position angle is ≈63°, or whether the stellar inclination is ≈180°− 30°= 150°and the position angle is ≈63°+ 180°= 243°. Regardless, it appears that κ And B is significantly misaligned on a prograde or retrograde orbit around its host star. However, for a fair incorporation with the rest of our sample, which does not possess information about the orientation of the star in the sky plane, we neglect the polar position angle and only include information about the rotational and orbital inclinations in our subsequent analysis of the κ And system. Figure 9. Line-of-sight stellar inclinations (i * , orange) compared with the orbital inclinations of substellar companions (i o , blue). In these "boomerang diagrams," stellar inclinations are mirrored about the x-y sky plane (i = 90°) because angular momentum vectors may point toward or away from the observer looking down along the z-axis (as depicted in the bottom-right panel of Figure 10). If i o agrees with i * (e.g., HD 49197), this implies that the host and companion are consistent with spin-orbit alignment; however, the true spin-orbit angle may nevertheless be nonzero in these cases. If i o and i * disagree (e.g., Gl 504), the system is misaligned by at least the difference in the inclination angles. See Sections 2 and 4.7 for details. difference, Δi, deviates from 0°. A Δi value of 0°is consistent with alignment, 0°< Δi < 90°is consistent with a prograde orbit, Δi = 90°is consistent with a polar orbit, and 90°< Δi 180°is consistent with a retrograde orbit. However, in all cases true obliquity angles (ψ) can be larger. The resulting distributions of Δi for all 21 systems are sorted by their median values in Figure 11. These differenced distributions take on a range of shapes: some are pressed up against i = 0°; some are bimodal, a reflection of the degeneracy of stellar inclination about the sky plane; and others significantly depart from alignment. Many distributions are quite broad and could be consistent with alignment, polar orbits, or retrograde motion. Table 3 summarizes the Δi maximum a posteriori (MAP) values, median values, and credible intervals.
We define two criteria to classify a system as being misaligned based on these Δi distributions. Most of the distribution power must be located at values of Δi away from 0°. But this alone would not be a sufficient metric because a flat distribution from 0°-180°, for instance, could satisfy this criterion but Δi would be unconstrained. So we also require that the MAP value is nonzero and beyond a given threshold. The characteristic uncertainty in stellar inclination measurements is about 10°, so we adopt this as a reasonable value to distinguish alignment from misalignment. Systems for which P (Δi > 10°) 95% and where Δi MAP values are beyond 10°a re deemed "misaligned." Those with P(Δi > 10°) 80% and Δi MAP > 10°are labeled as "likely misaligned"; most of their posterior power lay beyond that threshold. If neither are satisfied, the system is considered to be consistent with alignment. Finally, in addition to the minimum true obliquity distributions, we also overplot π − Δi, the maximum values of ψ, in Figure 11. In most cases these maximum values are near 180°, showing the broad range of potential angular momentum orientations in each system. Two companions in the sample are on significantly misaligned orbits. The orbital plane of TWA 5 B is offset with respect to the rotation axis of at least one component of the host binary TWA 5 Aab by at least 61°± 30°. The 95.4% credible interval of Δi spans 11°-110°. As with most systems, this implies that prograde, polar, and retrograde orbits are allowed. Note, however, that because TWA 5 Aab is a nearequal-mass resolved binary (Macintosh et al. 2001), it is possible that the stellar inclination could be impacted by the blended light curve and unresolved * v i sin measurement. Gl 229 B orbits its host in a nearly face-on orientation (i o = 5°.50 ± 0°.16; Feng et al. 2022), 17 while we find the star is viewed nearly equator-on. The resulting minimum misalignment value is Δi =  -+ 85 16 15 . Although the direction of the stellar rotation is not known, this does not impact the interpretation of the spin-orbit misalignment because in either scenario (whether it rotates clockwise or counterclockwise from the observer's perspective), the inclination is confined close to the sky plane. As a result, the maximum value of ψ (π − Δi) is close to the minimum value, so ψ ≈ 90°. Gl 229 B therefore appears to be on a polar orbit around Gl 229 A. This is the only system with this type of potentially well-constrained geometry in the sample.
We caution that this interpretation for the Gl 229 system is especially sensitive to the * v i sin value of Gl 229 A because its equatorial velocity is so small (1.0 ± 0.2 km s −1 ). It is difficult to reliably determine projected rotational velocities below about 2 km s −1 , so if the true * v i sin value is, for instance, 0.5 km s −1 , that would be challenging to measure. Our adopted * v i sin value is elevated above v eq , indicating that the projected rotational velocity may be overestimated-likely because it is near the systematics-dominated floor for broadening measurements. See Appendix B for a detailed discussion of this system. Additional in-depth rotational broadening analysis would help to clarify the inclination and spin-axis orientation of this system. For the purposes of this study we rely on these available measurements at face value, so a pole-on orbital orientation appears to be the most likely angular momentum architecture for this system, although this could change if the actual * v i sin value of the host star is less than ≈1 km s − i sin measurement, which impacts the width or the stellar inclination posterior, and the constraint on the orbital inclination. Improvements to both of these parameters will help narrow the resulting Δi distributions to more definitively establish the angular momentum geometry.
Altogether, this implies that a total of 11 out of 21 systems (52 -+ 11 10 % of the sample) show signatures of offsets between the rotational and orbital angular momentum vectors. Note that this measurement assumes that host star inclinations are equally likely to point toward or away from the observer (creating the boomerang structures in Figures 9-10). If mutual alignment is a priori expected (from particular disk-based formation channels, for instance), and this nonuniform prior probability were to be imposed on the stellar inclination distributions so as to make them asymmetric about 90°, this would impact the resulting Δi distributions and the implied misalignment fraction. As expected with Bayesian probabilistic inference, the answer and interpretation will depend to some degree on the priors. We simply emphasize here that the priors for the orientation of the rotational angular momentum vectors are identical and are not conditioned on whether the orbital inclination distribution points toward or away from the observer. We discuss the implications of this prevalence of misalignments in Section 5.
The remaining 10 systems are consistent with alignment: 1RXS J034231.8+121622, 2MASS J01225093-2439505, 51 Eri, HD 130948, HD 19467, HD 49197, HR 7672, κ And, 18 PDS 70, and PZ Tel. However, as is evident in Figure 11, the complementary distribution to Δi for the maximum value of ψ implies that there are a wide range of potential true obliquities each system can have. It is challenging to interpret individual systems with Δi distributions consistent with 0°because these represent minimum misalignments; their true obliquities can be much larger. In many cases, obliquities spanning any value from 0°-180°are consistent with the observations. Instead we focus the following discussion on misaligned (and likely misaligned) systems, and we develop simple population-level models to compare with our observed family of Δi distributions.

Comparison with Previous Obliquity Constraints
Several previous studies have explored spin-orbit alignment for individual systems with imaged substellar companions similar to the method we employ here: 1. Bowler et al. (2017) found that the orbit of ROXs 12 B is likely to be misaligned with the spin axis of its host star, with P(Δi > 10°) = 0.94. This is similar to the misalignment probability of 0.95 derived in this study. 2. Bonnefoy et al. (2018) carried out the same analysis for Gl 504 and noted possible signs of misalignment, with P (Δi > 10°) = 0.78. We infer a higher probability of 0.91, most likely because of the more precise * v i sin we adopt as well as the updated orbit constraint. 3. Maire et al. (2020) found that the brown dwarf companion HD 19467 B is consistent with alignment. We reach a similar conclusion, with P(Δi > 10°) = 0.75, which is below the threshold of 80% that we adopt for classifying systems as being "likely misaligned. (2021) analyzed the orbital configuration of GQ Lup B with respect to the host star's spin axis, the companion's spin axis, and the orientation of the circumstellar disk. Although the orbit is not yet well constrained, both studies conclude that nonzero stellar obliquity is possible. We find that spin-orbit misalignment is likely in the GQ Lup system, with P(Δi > 10°) = 0.896.
Altogether our results from this study are in good agreement with previous work focused on individual systems.

Implications of Misalignments
Most of our sample comprises brown dwarf companions ( Figure 12), but there are two stars hosting giant planets (PDS 70 and 51 Eri) as well as a handful of objects whose masses and origin are more ambiguous (such as VHS J1256-1257 b, 2M0122-2439 b, AB Pic b, Gl 504 B, and GSC 6214-210 B). Because the overall sample size is modest, and individual constraints are fairly broad, we consider our sample as belonging to one substellar population of (predominantly) brown dwarf companions for this analysis and discussion about formation channels.
Brown dwarf companions are expected to form like stellar binaries from the turbulent fragmentation of molecular cloud cores or through disk fragmentation (e.g., Low & Lynden-Bell 1976;Bate et al. 2002;Bate 2009;Stamatellos & Whitworth 2009;Jumper & Fisher 2013). In principle, these formation sites-in cores, filamentary structures, or in diskscan result in similar orbital and rotational axes for the host and companion through a localized conservation of angular momentum during collapse. However, a variety of mechanisms have been proposed that can disrupt this alignment during or after the star formation process. This includes nonaxisymmetric collapse from variable accretion, interactions with nearby protostars, or inhomogeneities in the parent cloud core (e.g., There have been few large-scale empirical tests of these various models owing to the difficulty in determining obliquities for stars hosting long-period brown dwarf companions. Observationally, close stellar binaries have been shown to be well aligned in many individual instances (e.g., Sybilski et al. 2018), although there are several notable examples of spinorbit misalignments (e.g., Albrecht et al. 2009). Ensemble obliquity measurements of binary systems have yielded tentative or inconclusive results (Hale 1994;Justesen & Albrecht 2020). Recently, however, Marcussen & Albrecht (2022) found that most close binaries in their sample of 43 systems are consistent with alignment.
The most significant statistical result from this analysis is that spin-orbit misalignments appear to be common among stars hosting wide substellar companions. The (minimum) fraction of systems likely to be misaligned, -+ 52 11 10 %, is higher than for cool stars hosting hot Jupiters (∼5%; Winn et al. 2010;Albrecht et al. 2022) and close-in small planets (Campante et al. 2016;Winn et al. 2017). It also appears to be higher than the misalignment fraction for stars hosting debris disks (10%; Le Bouquin et al. 2009; Watson et al. 2011;Greaves et al. 2014), a comparison that may be more relevant because the spatial scales of tens to hundreds of astronomical unit are closer to the population of companions considered in this study. Since debris disks are expected to trace the sites of planetesimals and planet formation, this disagreement between the incidence of misaligned debris disks and misaligned brown dwarf companions may indicate that companions in our sample did not predominantly form in disks; if this were the case, the alignment frequencies would be expected to be similar. Less is known about the overall coplanarity of systems that contain both spatially resolved debris disks and substellar companions. HD 2562 B (Konopacky et al. 2016;Maire et al. 2018) and HD 206893 B (Delorme et al. 2017;Milli et al. 2017;Marino et al. 2020;Ward-Duong et al. 2021) are notable examples of brown dwarf companions that appear to be aligned with their exterior debris disks-suggesting a disk-based origin for these companionswhereas planets orbiting HR 8799, β Pic, and HD 106906 show a range of orientations with respect to exterior and interior disks (e.g., Dawson et al. 2011;Bailey et al. 2014;Wilner et al. 2018;Nguyen et al. 2021). A larger sample would help clarify the prevalence of aligned and misaligned orientations.
An alternative explanation is that subsequent post-formation scattering with a second object in the system (so as to alter the initial obliquity distribution) may be more common for stars hosting wide substellar companions, perhaps because multiple massive objects formed in the same disk. With a few exceptions (Haffert et al. 2019;Lagrange et al. 2019;Bohn et al. 2020b;Hinkley et al. 2022), however, most follow-up searches have not identified additional giant planets and brown dwarfs that could represent potential scatterers (e.g., Bryan et al. 2016).
It is also possible that some widely separated substellar companions could have been captured at an early age while still embedded in a dense cluster. Close encounters with passing stars can liberate planets and brown dwarfs (e.g., Parker & Quanz 2012;Daffern-Powell et al. 2022); these free-floating planets, and other isolated brown dwarfs formed during the star formation process, could then become captured onto very wide orbits by other stars, especially during the cluster dispersal phase (e.g., Perets & Kouwenhoven 2012;. Without any significant realignment mechanism at these large orbital distances, the stellar obliquity distribution in such cases would be expected to be isotropic. The misaligned systems in our sample are therefore also consistent with this capture process.
There is no simple interpretation of this high misalignment fraction. If circumstellar disks are predominantly aligned with the spin axes of their host stars, this result points to a diversity of formation channels or dynamical processing in some systems. Alternatively, the observed distribution of spin-orbit orientations could reflect the primordial distribution of aligned and misaligned disks. Systems with spin-orbit misalignment are also consistent with the early capture of free-floating giant planets and brown dwarfs in dense clusters. A better understanding of the initial conditions for disks around cool stars would provide additional context to interpret results from this study.

Stellar Obliquities in Systems with Directly Imaged Planets
Alignment rates are less clear among systems with directly imaged planets, largely because long-period giant planets tend to reside around early-type stars whose rotation rates cannot be accurately determined through photometric monitoring in the same way as can be inferred for cool stars . However, in some cases stellar inclinations and true three-dimensional orientations can be constrained using asteroseismology or spatially resolved interferometry. Below we summarize orbital and rotational angular momentum alignment in five planetary systems with previous obliquity constraints: β Pic, HR 8799, 51 Eri, HD 206893, and PDS 70. Note that among these, only 51 Eri and PDS 70 are in our sample of 21 stars with obliquity constraints; the rest have had previous spin-orbit measurements obtained using a variety of alternative techniques. Kraus et al. (2020) analyzed the spatial displacement of the A6 host star β Pic in the Brγ absorption line using VLTI/ GRAVITY spectro-interferometry. When combined with the asteroseismic stellar inclination from Zwintz et al. (2019), they found that β Pic is well aligned (to within 3°± 5°) with the orbit of β Pic b and the outer debris disk. More recently, a second giant planet in the system, β Pic c (Lagrange et al. 2019), was directly imaged and shown to be consistent with this coplanar geometry Nowak et al. 2020;Brandt et al. 2021a;Lacour et al. 2021).
The orbital inclinations (and longitudes of ascending node) of the four planets orbiting HR 8799 (Marois et al. 2008(Marois et al. , 2010 have progressively improved over time with continued orbit monitoring. The precision of the constraint depends on assumptions about resonant orbits, stability, and coplanarity, but in the unconstrained case the four orbital planes lie between about 20°and 40°(e.g., Wang et al. 2018;Sepulveda & Bowler 2022); in the most restricted case, Zurlo et al. (2022) found an inclination of  -+ 26.9 0.17 0.16 for the planetary system. Wright et al. (2011) measured an inclination of 40°for the host star using oscillation frequencies from high-cadence RVs, which would imply a spin-orbit misalignment by at least 10°. However, Sodor et al. (2014) concluded that asteroseismic analysis of HR 8799 as an A5/F0 γ Dor variable is challenging because of resonances and amplitude variations among the pulsation frequencies. The inclination and obliquity of HR 8799 are therefore not yet reliably established, although an estimate of the core rotation period by Sepulveda et al. (2023) implies a stellar inclination in good agreement with the orbits of its four planets. Furthermore, the exterior cold disk in this system appears to be coplanar with the planets' orbits (e.g., Faramaz et al. 2021).
Maire et al. (2019) derived a stellar inclination of the young F0 planet host 51 Eri (Macintosh et al. 2015) using the equatorial and projected rotational velocity method. The equatorial velocity was based on a rotation period of 0.65 ± 0.03 days computed using Hipparcos photometry. They conclude that the stellar inclination (39°-51°or 129°-141°) and orbital inclination of 51 Eri b (126°-147°) are consistent with orbital and rotational alignment. However,  recently analyzed the TESS light curve of 51 Eri and identified it as a γ Dor pulsating star. The variability previously attributed to rotational modulation is more likely caused by pulsation modes. Although the surface rotation rate is not constrained, the authors estimate a core rotation period of -+ 0.9 0.1 0.3 days and use this to deduce a stellar inclination of ∼62°. We make use of this core rotation period estimate here; based on our stellar inclination of  -+ 54 19 13 we conclude there is no evidence of misalignment.
HD 206893 hosts a brown dwarf companion (Milli et al. 2017) and giant planet (Hinkley et al. 2022) on orbits consistent with being coplanar. Delorme et al. (2017) presented ground-based photometric monitoring of the F5 host star and found clear modulations with a period very close to 1 day. They interpreted this signal of 0.996 ± 0.03 days as most likely originating from a rotation period. They then use the equatorial and rotational velocities to infer a stellar inclination of 30°± 5°, which (given the symmetric nature of the distribution about 90°), is consistent with the orbital inclination of 146°± 3°for HD 206893 B and 150°± 3°for c (Hinkley et al. 2022). This is also consistent with measurements of 40°± 3° (Marino et al. 2020) and 45°± 4° (Nederlander et al. 2021) for the orientation of the debris disk. However, Delorme et al. (2017) noted that the photometric periodicity they measured could be attributed to stellar pulsation modes if HD 206893 is a γ Dor variable star. This is certainly possible as F5 is near the boundary of the onset of this instability strip (Kaye et al. 1999;Aerts 2021). While the ≈1 day rotation period is plausible, additional monitoring of this system would be valuable to more robustly establish the nature of the periodic variations. For the time being, we therefore cautiously interpret the angular momentum geometry of this system as being potentially-but not securely-aligned. 19 PDS 70 hosts two accreting giant planets nested inside a transition disk Haffert et al. 2019;Zhou et al. 2021). The orbits of PDS 70 b and c have been constrained to  -+ 133 3 4 and 132°± 3°, respectively, assuming noncrossing orbits and near-coplanarity ). Thanathibodee et al. (2020) measured a rotation period of 3.03 ± 0.06 days for PDS 70 from the TESS light curve and inferred a stellar rotation axis of 50°± 8°(equal to 130°± 8°) -in good agreement with the orbits. We find a similar result consistent with spin-orbit alignment: i * =  -+ 60 15 10 , Δi = -+ 43 43 24 , and P(Δi > 10°) = 0.804. These orbital and rotational orientations are also in agreement with those of the transition disk (Keppler et al. 2019). The PDS 70 system therefore appears to be in a state of organized angular momentum alignment.
In summary, a variety of constraints are now available on the true obliquities or minimum misalignments for a handful of stars hosting directly imaged planets. The β Pic system is unambiguously aligned in ψ. PDS 70 is consistent with spinorbit alignment in Δi. 51 Eri and HD 206893 are potentially aligned with their companions. The stellar inclination and obliquity of HR 8799 is not reliably constrained.
In Figure 12 we compare the minimum misalignment constraint for each system as a function of companion mass and separation. β Pic b and c are included based on their true obliquity constraint from Kraus et al. (2020) and planet properties from Nowak et al. (2020). The sample size is modest, but we note that all four of the companions in three systems with masses below 10 M Jup are consistent with spinorbit alignment, and, similarly, seven of the eight companions with separations less than 20 au are consistent with alignment. At higher masses and wider separations, there is a mixture of aligned and misaligned systems with no clear correlation with increasing brown dwarf mass or orbital separation.
Altogether, this points to an emerging trend that systems with long-period giant planets have low stellar obliquities, although the small sample size of (primarily early-type) host stars with inclination measurements remains very small. This would, however, have important implications for the formation of these companions. A uniform analysis of eccentricities by Bowler et al. (2020a; see also Nagpal et al. 2023) revealed that imaged planets between 5 and 100 au have low eccentricities, whereas brown dwarf companions have a broader range of eccentricities. If both of these trends are validated with larger samples, this would reinforce a distinction in the orbits, obliquities, and formation of these two populations: long-period planets are most consistent with a disk-based origin with minimal dynamical evolution, whereas brown dwarf companions as a population are more aligned with formation from turbulent cloud fragmentation.

The Underlying Distribution of Δi with Hierarchical Bayesian Modeling
Although it is difficult to make a meaningful constraint on the underlying true obliquity distribution, we can nevertheless compare our observed sample of Δi distributions with extreme cases of isotropic stellar orientations and completely aligned spinorbit orientations, as well as a mixture of these two scenarios. We generate three synthetic data sets for these purposes.
For the isotropic case, we randomly draw a stellar inclination from a * i sin distribution from 0 to π/2. Assuming the measurement of i * follows a Gaussian with a standard deviation of 10°, which is a characteristic level of uncertainty from this study, we then mirror the distribution across 90°to reflect the unknown orientation of the spin angular momentum vector. For each system we use the measured i o and our new i * distributions to compute a realization of a sample of Δi distributions assuming isotropic stellar line-of-sight inclinations.
For the aligned case, the procedure is the same as the isotropic case, except instead of selecting a randomly drawn mean value of i * , we use the MAP value of i o . If this lies between 0 to π/2, then i * is mirrored onto the π/2 to π range, and vice versa. The absolute difference then provides a distribution of minimum obliquities assuming perfect spin-orbit alignment and an unknown orientation of ¾  * L to mimic the information available to an observer. For the mixture case, we randomly draw 10 unique distributions from the isotropic sample described above and nine from the aligned sample to create an approximately even mixture of 19 systems. This represents an intermediate scenario between the two extreme examples of completely aligned and independent angular momentum architectures.
The observed and synthetic limiting cases are shown in the upper panels of Figure 13. Δi distributions for the isotropic case are broadly located between about 0°and 130°. As expected, the aligned case is clustered at Δi = 0°with many individual distributions being bimodal because i * is often bimodal. Altogether, the observed sample better resembles the isotropic or mixed cases with many distributions peaking at nonzero Δi values.
To quantify this agreement, we employ hierarchical Bayesian modeling to reconstruct the underlying populationlevel distributions of the real and synthetic samples. Our general strategy follows the framework described by Hogg et al. (2010), in which an importance sampling approximation of the likelihood function enables a phased approach to the problem: individual-level posteriors can be generated first and subsequently sampled to constrain hyperparameters of an assumed population-level model.
Here we adopt a flexible model introduced by Iriarte et al. (2021) called the Lambert generalized bimodal (LGB) distribution, which was developed to introduce skewness to the generalized bimodal distribution. The LGB distribution has four parameters: the location parameter μ (defined for all real numbers), the scale parameter σ (defined for positive numbers), and the shape parameters γ (defined over [0, 2)) and α (defined over (0, e), where e is Euler's constant). Depending on these parameters, the LGB distribution can be unimodal, bimodal, symmetric, or asymmetric; examples of the diversity of shapes can be seen in Figure 1 of Iriarte et al. (2021). γ controls the bimodality of the distribution, and α determines the relative heights of the bimodal peaks or unimodal asymmetry.
This choice of the underlying model is motivated by the bimodal nature of many Δi distributions, which tend to produce power at secondary peaks, as well as the versatility and modest number of hyperparameters of the LGB distribution. It is not, however, otherwise physically motivated, and certainly other distributions or mixture models could be used in its place to capture the behavior of the (unknown) distribution from which this sample was drawn. The exact shape and inferred parameter values are less important for this exercise than a comparison of the results from the observed sample with the isotropic, aligned, and mixture cases.
The LGB probability distribution function for Δi takes on the following form: , and Φ(z) is the cumulative distribution function of the standard normal: ( ( )) + z 1 erf 2 2. Hierarchical Bayesian modeling of all three samples is carried out in a similar fashion as in Bowler et al. (2020a) for the case of the two-parameter Beta distribution, except here the posteriors of four hyperparameters are sampled. We implement a Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm (Metropolis et al. 1953;Hastings 1970) with tunable Gaussian jump parameters and Gibbs sampling. Uniform hyperpriors are adopted for all four hyperparameters: ( ) m  ∼ ( ) -   1000 , 1000 , so as to allow the location parameter to drift beyond the range of the data, . 20 A single Markov Chain ran for 10 5 steps with parameter acceptance rates generally falling between 20% and 50%. Trace plots are visually analyzed to assess convergence, and a burn-in fraction of 30% is adopted.
Results are shown in the bottom panels of Figure 13. The family of hyperparameter posterior distributions is clustered between Δi values of ≈0°-50°for the observed sample Following Nagpal et al. (2023), we define the metric  to quantify the level of agreement between the family of underlying distributions in the observed sample and the three synthetic cases:  is the average area of the absolute difference between random posterior draws (n) from the underlying model of the observed sample, f (Δi) obs,n , compared with the isotropic, aligned, or mixture models, f (Δi) model,n . Lower values of  correspond to better agreement between the two underlying families of distributions. We also quantify the spread in this metric with the standard deviation of these absolute residual areas: We find  and s  values of 0.83 ± 0.30 when comparing the observed and isotropic case, 2.00 ± 0.46 for the observed Figure 13. Top: Δi distributions for our observed sample of substellar companions (left) compared with one synthetic realization for an isotropic stellar inclination distribution (middle left), perfect spin-orbit alignment (middle right), and a random mixture of synthetic isotropic and aligned distributions in equal proportions (right). The observed sample is qualitatively most similar to the isotropic or mixture cases, but differs substantially from a purely aligned scenario. Bottom: reconstructed underlying distributions for each sample using hierarchical Bayesian modeling. Here the flexible Lambert generalized bimodal (LGB) distribution is used for our population-level model, and constraints on hyperparameter posteriors are listed in each panel. The population-level distribution for the observed sample closely resembles a mixture of isotropic and aligned systems. The thick curve shows the LGB distribution using median hyperparameter values, while thin curves show 30 draws from the MCMC chains. Note that the scale has been adjusted in the aligned panel (middle right). 20 Nagpal et al. (2023) analyzed the eccentricity distribution of long-period substellar companions and found that hyperpriors on the two Beta distribution shape parameters can impact the hyperparameter posteriors in counterintuitive ways. Specifically, a uniform hyperprior can result in a biased eccentricity distribution, and the wider the range of the shape parameter hyperprior, the narrower the resulting posteriors. However, in this case for the LGB distribution, the parameters μ and σ directly map the location and spread of the distribution function, and the two shape parameters γ and α by definition have limited ranges. We therefore do not expect uniform hyperpriors to bias our results in the same way as they would for the Beta distribution, although more study on the effect of hyperpriors for this distribution is needed. 21 Although all four cases show some asymmetric structure in the recovered distributions, we do not necessarily ascribe this as a real feature, as this is likely to simply be an artifact of a weakly constrained γ parameter. and aligned case, and 0.67 ± 0.24 for the observed and mixture case. The inferred distribution for the observed sample is in much better agreement with the mixture scenario than either the pure aligned or isotropic cases. Indeed, all four parameters in the LGB distribution model are consistent at the 1σ level with the mixture sample. The aligned scenario has the poorest agreement among the three simulated samples that we tested. Although we do not attempt to constrain the relative proportion of isotropic versus aligned orbits that best fit the observed sample, the good agreement with the equal mixture suggests that there are more aligned cases in our observations than would be expected for a purely isotropic distribution. 22 Overall, there does not appear to be a strong correlation between the orbital and rotational planes for systems with wide substellar companions. These results reinforce the view that misalignments are common among brown dwarf companions. A larger sample and more precise individual constraints will help to clarify this picture in the future.

Sample Biases
There are several observational biases that are important to highlight as they are likely to impact the sample selection from this study with varying degrees of significance. These are individually detailed below and summarized in Section 5.5.8.

Light-curve Selection Biases
Because our strategy to constrain i * relies on rotation period measurements, any biases that impact our ability to detect and interpret periodic photometric modulations will influence the types of host stars we can analyze. For instance, light curves with a limited time baseline are inherently more sensitive to brighter stars with larger-amplitude variations and shorter rotation periods. Some of the TESS light curves we rely on comprise only a single 27 day sector (see Table 2), which means we would only be able to reliably establish periods on timescales of about half that baseline (∼13 days). The Sun has a small surface spot covering fraction, long equatorial rotation period (24.5 days), and brightness variations typically under 0.1% (Reinhold et al. 2020). Similarly old G dwarfs would not be recovered in TESS, so inclination and obliquity measurements for these types of host stars would not be possible using methods that rely on photometric rotation periods. As a result, there is very likely a selection bias in this study in favor of systems with brighter, younger, and later-type host stars in our sample, as has been recently established by Masuda (2022) for Sun-like stars in the Kepler field. This age bias is already present for systems with directly imaged planets, but it is likely reinforced for stars with brown dwarf companions as well by disfavoring older stars.

Pole-on Orientations
We also expect there to be a geometric bias associated with selecting stars based on the presence of light-curve modulations. For stars with solar-like dynamos, spot distributions located in equatorial regions would be more visible and would produce higher-amplitude variations if seen in an equator-on orientation. Pole-on viewing angles are inherently geometrically unlikely, but for Sun-like stars, this dearth of pole-on viewing angles is likely to be reinforced. However, highlatitude and polar spots are commonly observed on the surfaces of rapidly rotating active stars (Strassmeier 2009;Yadav et al. 2015;Roettenbacher et al. 2016) as a result of youth, spin-orbit coupling in close binaries, or low stellar masses with large convective envelopes. If the dominant spot structure is located in a large polar cap-like spot, and the evolution of that spot is slow, the impact on brightness variations in a light curve would be small for a star viewed pole-on compared to a similar spot rotating in and out of view in equatorial regions. There may therefore be a two-fold selection bias against pole-on orientations near i * = 0°as a result of the lack of polar spot structures for Sun-like stars and long-lived polar spots for active stars in the sample.
This possible bias against pole-on systems can be directly tested with our sample. Three out of 53 stars with inclination constraints in Table 3  = 0.015. The expectation value in a sample of 53 stars is 0.8, so we would expect about one star to be viewed this way. This is similar to the number of near-pole-on stars in our sample. We therefore do not find strong evidence for a significant bias disfavoring pole-on orientations.

Differential Rotation
Differential rotation is expected to be common and has been observed in many stars. Starspots located at mid-or highlatitudes can result in a systematically longer inferred rotation period than the true equatorial rotation period. The severity of this bias is increased for spots located at higher latitudes and stars with stronger levels of shear. An overestimated rotation period will underestimate the equatorial velocity, which in turn will overestimate the stellar inclination when using projected rotational velocity to infer i * .
For example, the active regions in the Sun are confined to latitudes of about ±30°where the rotation period is ≈25.7 days (based on the latitude-dependent relations for differential rotation from Snodgrass & Ulrich 1990). This is only 5% higher than the equatorial rotational velocity of 24.47 days, so the effect is not particularly large for old Sun-like stars with solar-like dynamos. However, using asteroseismology of stars in the Kepler field, Benomar et al. (2018) found that some solar-type stars rotate twice as fast at their equators compared to their mid-latitudes. A broad distribution of absolute shear also appears to be present among low-mass stars, even though the overall level of differential rotation decreases with effective temperature (Barnes et al. 2005;Reinhold et al. 2013;Reinhold & Gizon 2015). Although most stars are expected to follow a solar-like differential rotation law with shorter rotation periods at the equator and longer periods at the poles, antisolar rotation (with poles rotating faster than the equator) and cylindrical rotation may also be possible (Brun et al. 2017). Altogether, there is considerable uncertainty in the detailed picture of differential rotation and the overall impact of spot distributions on light curves. While we have made conservative assumptions about the uncertainties associated with our photometric rotation periods (see Section 4.2), we did not make explicit adjustments 22 It is interesting to note that the stars with the two lowest-mass companions in the sample-51 Eri and PDS 70-are consistent with alignment. If these were excluded, it is likely that the underlying distribution for the observed sample would broaden out to better reflect the isotropic distribution. That is, the mixture model might be preferred in part because our sample itself is a mixture of brown dwarf companions and giant planets.
to the period values to correct them for potential systematic overestimates.
One approach to assess whether there is a bias toward high inclinations is to check for an overrepresentation of equator-on systems in our sample in the same fashion as we did in Section 5.5.2. For this analysis, we adopt MAP values of i * above 80°as a threshold for equator-on viewing angles. There are 19 stars in the sample with i * > 80°. The probability of a random orientation falling in that range is ò   = = * * * * i di sin i i 80 90 = 0.174, implying an expectation value of ≈9 for a sample of 53 stars. It does therefore appear that there are more equator-on stars than expected from chance alone. To quantify the significance of this discrepancy, we can treat this outcome as a Bernouilli trial and compute the probability of an event occurring that is at least as extreme as what was measured (k = 19 out of n = 53 systems) using binomial statistics. Here the probability that a success occurs is p = 0.174. This statement is equivalent to computing 1 minus the probability of observing fewer than 19 stars in the 80°< i * < 90°range: P (k 19|p = 0.174, n = 53) = 1 − P(k < 19|p = 0.174, probability of there being so many equator-on systems by chance is therefore very small. This can be interpreted either as confirmation that the overestimated rotation period bias is impacting our sample, or alternatively as a result of the Bayesian formalism we adopt, in which a broad or unconstrained * v i sin value relative to v eq will revert toward the * i sin isotropic prior in which i * is near 90°. Of course, both explanations could be at play here.
To examine whether poorly constrained * v i sin values are driving the high incidence of equator-on stars, we compare the distribution of s * v i sin / * v i sin values from the full sample in Table 2 to the subset of 19 stars with i * > 80°. If the bias in line-of-sight stellar inclinations is caused by poor constraints on projected rotational velocities, we would expect the subsample to be shifted to higher values of s * However, this is not what we find: the median precision of * v i sin for the full sample is 3% (with a standard deviation of 13%), and for the subsample is 5% (with a standard deviation of 13%). This suggests that overestimated rotation periods, likely caused by differential rotation, are artificially increasing some of the stellar inclinations in our sample.

Inclination Discovery Bias
Companions found through high-contrast imaging are only visible outside of an inner working angle (IWA), which represents the closest angular separation at which a detection can be made for a given contrast. The IWA is often defined by the edge of a coronagraph or the radial extent of residual features after primary subtraction. For a star at a given distance and a companion on a circular orbit with a given semimajor axis located just outside of the IWA, a face-on orbital orientation will be more readily discoverable in a blind survey than a companion on an inclined orbit. The same companion on an edge-on orbit will spend the least time outside the IWA. This "inclination discovery bias" will impact systems with angular semimajor axes (α) close to the IWA, or α/IWA ∼ 1 (e.g., Janson 2010). Bowler et al. (2020a) examined the impact of (the reciprocal of) this ratio to assess whether an analogous "eccentricity discovery bias" was present in their sample of companions between 5 and 100 au with orbit constraints. They found that most systems at the time of discovery had α/IWA values between 1.2 and 10, and therefore that this discovery bias did not play a strong role in shaping the eccentricity distribution of their sample. Thirteen out of 21 systems with obliquity constraints in this study overlap with the sample in Bowler et al. (2020a) that was used to analyze the impact of discovery bias. The additional systems in this study are generally at separations slightly beyond the 100 au cutoff used in that study and generally have large α/IWA values much greater than 1.0.
Here we examine whether an inclination discovery bias might have impacted the orbital properties of substellar companions that have been discovered in high-contrast imaging surveys (e.g., Bowler 2016; Nielsen et al. 2019;Vigan et al. 2021), which in turn could influence the minimum obliquity distribution Δi in this study. We simulate the astrometric orbits of companions and randomly sample their projected separations on the plane of the sky as a function of orbital inclination, semimajor axis (in units of α/IWA), and eccentricity. For each orbit, we randomly draw a longitude of ascending node from 0 to 2π, argument of periastron from 0 to 2π, and time of periastron passage. The following α/IWA values are examined: 1.01, 1.1, 1.2, 1.5, 3, and 6. 10 5 synthetic observations of companions on unique orbits are generated for each set of sampled parameters; the proportion of observations for which the companion's projected separation is beyond 1 α/IWA (constituting a "detection") represents the recovered companion fraction. Results are displayed in Figure 14.
Considering the case of circular orbits, there is a steep decline sensitivity for more inclined systems and an increase in detection sensitivity for companions on wider orbits. The recovery rate is nearly 100% for α/IWA values beyond ≈3, with only the most edge-on orbits experiencing some loss in sensitivity. As soon as nonzero eccentricity is introduced, there is a general increase in sensitivity for a given α/IWA value, except for near-face-on orbits where a slight drop is evident. This behavior is readily explainable: companions on eccentric orbits will reach farther apastron distances and will spend most of their time at those locations. But for face-on circular orbits, where detection sensitivity is maximized, even a slight eccentricity will bring part of the orbit inside the IWA. This trend of overall sensitivity with increasing eccentricity persists even to unrealistically high eccentricities above 0.9. There is also a clear inclination bias present: companions on more edgeon orbits are recovered at ≈2-3 times lower rates compared to companions on more face-on orbits, with this drop occurring beyond inclinations of about 40°for low eccentricities and beyond 60°for high eccentricities.
However, these results do not take into account the relative geometric probability of observing face-on versus edge-on systems. To account for this effect, each curve should be adjusted by a factor of * i sin , as shown in Figure 15. The impact is to flatten the recovered companion fraction curves so that for any given semimajor axis, there is a roughly uniform probability of observing any orbital inclination between the * i sin envelope of isotropic orbits to 90°. This behavior persists to about e = 0.5, then begins to more closely resemble an isotropic distribution at higher eccentricities, especially for large values of α/IWA where sensitivity to companions is nearly complete. In this limit, the isotropic distribution is recovered. These trends should better represent more realistic biases on the orbital inclinations of discoveries made through high-contrast imaging. The actual distribution of observed orbital inclinations will depend on the underlying range of semimajor axes, host star distances, eccentricities, and IWAs.
While not the focus of this study, we note that these "truncated isotropic inclination distributions" should be more appropriate to use as Bayesian priors for the orbital plane inclinations ( ( )  i o ) of newly discovered companions located near the observational IWA than the more traditional ( ) =  i i sin o o inclination prior that is widely adopted for orbit fitting. Below is an approximation of these more realistic priors that account for the general behavior of inclination discovery bias as a function of semimajor axis and orbital eccentricity: γ is the inclination in radians beyond which the distribution is constant for increasing values of i o . For this approximation, we implement a break from the sine curve that depends on both the angular semimajor axis α, here parameterized as β = α/IWA, and eccentricity e: 0.5 cosh 1.5 log 1 . 15 1 2 Here the inverse hyperbolic cosine term is a better fit to the γ break values as a function of β parameter than a polynomial. The second term is motivated by a logarithmic-like scaling of γ with eccentricity and β. Coefficients were selected so as to reasonably agree with values in Figure 15 as a compromise between consistency with the numerical simulations and analytical simplicity. In the limit of wide orbits (α  11.6 IWA), γ > π/2 and the standard isotropic distribution is recovered. The choice of e and β may depend on the context of the system. For example, e = 0.0 may be an appropriate choice for a new directly imaged planet, with β corresponding to the projected separation over the IWA, whereas a moderate eccentricity of e ≈ 0.5 may be a better choice for a brown Figure 14. Assessing biases in orbital inclinations for direct imaging surveys. The sensitivity to faint substellar companions on various orbital inclinations depends on orbital eccentricity and the ratio of the angular semimajor axis, α, and the inner working angle (IWA), which establishes whether a companion is detected or obscured (by a circular coronagraph, for example). Imaging is more sensitive to companions with low inclinations, wider orbits (higher values of α/IWA), and higher eccentricities. These trends are corrected for isotropic orbital viewing geometry in Figure 15. dwarf companion Nagpal et al. 2023). Examples of these unnormalized priors are shown in Figure 16.
How do these expectations compare with actual orbital inclinations? A histogram of the median orbital inclinations from Table 3 is shown in Figure 17 compared to an isotropic distribution. Because the sample size is modest, we have reflected the inclination values (which span 0°-180°) about 90°a s the predicted yields are symmetric with respect to the sky plane. The observed inclinations seems to match expectations reasonably well; there is good agreement with the isotropic case at low inclinations, and the overall distribution is approximately constant from ≈40°to 90°. Altogether, this can be explained by a mixture of modest eccentricities of brown dwarfs and the range of α/IWA values from Bowler et al. (2020a).

Orbital Inclination Priors
Another relevant form of inclination bias is related to the challenge of fitting orbits with small phase coverage. Most substellar companions with orbit constraints have had less than 15% of their orbits mapped with relative astrometry ). This makes the inferred orbit solutions and the detailed shape of the orbital element posteriors sensitive to both random and systematic uncertainties in relative astrometry, as well as the choice of parameter priors used in the analysis. Recently, Ferrer-Chávez et al. (2021) analyzed systematic biases associated with fitting orbits to companions with very small coverage (≈0.5%). They found that face-on inclinations can result in significant biases in the posterior median and mode, and that for progressively higher eccentricities, more inclined orbits can be impacted by this bias. A similar conclusion was found by O'Neil et al. (2019); without adopting their observable-based priors, inclinations are broader and less accurate. For our sample in this study, these effects can be exacerbated by the use of different approaches to orbit fitting in the literature and parameter priors that were adopted. It is therefore challenging to determine the exact scale of inclination bias caused by priors, although it is expected to be more severe for orbits with smaller fractional coverage.  Figure 14 but now taking into account the relative probability of observing companions with random orientations. When accounting for this * i sin isotropic inclination distribution, the sensitivity to a companion is roughly constant for a given semimajor axis from the * i sin envelope to i ≈ 90°. This results in a "truncated isotropic inclination distribution" with fewer companions having higher (more edge-on) inclinations and represents a more realistic outcome for the expected orbital inclination distribution of discoveries from blind direct imaging surveys.

Slow Rotation
As discussed in Section 4.7, it is challenging to reliably determine inclinations for slow rotators with equatorial velocities less than a few kilometers per second because this is near the limit where traditional spectroscopic approaches to measuring rotational broadening, macroturbulence, and microturbulence can be difficult to distinguish (e.g., Smith & Gray 1976;Gray 2005;Petigura et al. 2017;Masuda et al. 2021). The result of this is a tendency to interpret slow rotators as being equator-on, even if their actual inclinations substantially depart from 90°. This may contribute to the overabundance of stellar inclinations between 80°and 90°discussed in Section 5.5.3. However, only seven of 53 stars with i * constraints have v eq < 3 km s −1 , so we do not expect these systematic errors to severely impact statistical results from this study.

Relative Biases
In addition to overall biases that may be impacting our sample selection and inferred inclination distributions, relative biases may also be present within our sample that could in principle influence the interpretation of the misalignment results in this study. In particular, the apparent distinction between obliquities of stars hosting long-period giant planets and brown dwarf companions could instead be caused by obliquity evolution, for example through dynamical excitation over time. Imaged planets are typically young (200 Myr), whereas the systems with brown dwarfs span a much broader range of ages because brown dwarfs remain bright enough to readily detect even at field ages of several gigayears. If obliquities are initially low but evolve toward a broader distribution over time, perhaps through dynamical interactions, that could manifest as a distinction between planets and brown dwarf companions as a result of their different characteristic ages. However, our sample only comprises a few field-age brown dwarfs-Gl 229, HR 7672, and HD 19467-so if obliquities are evolving, then this mechanism would likely need to operate on relatively short timescales (1 Gyr) corresponding to the typical age range of systems with brown dwarfs in this study.
Other possible sources of relative biases within our sample could include a dependence of obliquity on stellar host mass, metallicity, or formation environment. All of these could be explored in more detail with larger samples in the future.

Summary of Biases
Below we summarize the most important biases likely to impact stellar inclinations, orbital inclinations, and minimum obliquity values (Δi) in our sample.
1. We expect a strong bias in favor of active stars with younger ages (faster rotation periods) and cooler temperatures (greater starspot covering fractions) as a result of the finite photometric precision and limited baseline coverage of light curves. 2. Both differential rotation and slow rotation can result in a systematic overestimation of stellar inclinations. This may explain the overrepresentation of equator-on orientations that is present with high confidence in our larger sample of 53 host stars with inclination constraints. 3. A discovery bias affecting the orbital inclination distribution is expected for companions found with high-contrast imaging. This distribution is not isotropic but instead more closely resembles a truncated i sin o distribution. This deficit of edge-on orientations is reflected in the orbital inclination distribution of companions in our sample. 4. Bayesian priors may be influencing the detailed shape of orbital inclination posteriors for systems with more faceon orientations and less orbital phase coverage compared to systems with well-mapped orbits. Figure 16. Example of parameterized Bayesian priors for orbital inclinations that take into account the inclination discovery bias caused by the presence of a hard IWA in high-contrast imaging data sets (Equations (14) and (15)). Note that for separations very close to the IWA, the expected orbital inclination prior deviates substantially from a traditional isotropic ( i sin o ) distribution. For circular orbits, this can resemble a uniform distribution. These modified isotropic distributions are unnormalized approximations to the results in Figure 15 and depend on the angular semimajor axis α and orbital eccentricity. Because long-period giant planets and brown dwarf companions have distinct eccentricity distributions Nagpal et al. 2023), a different choice for the characteristic eccentricity might be warranted depending on the properties of a newly imaged companion-for instance, e = 0.0 for planets and e = 0.5 for brown dwarfs. Figure 17. Distribution of orbital inclination MAP values for the 21 substellar companions with obliquity constraints examined in this study. Compared to an isotropic inclination distribution (gray), the observed sample has a deficit of objects with edge-on orientations between ≈50°and 90°. This closely resembles predictions in Figure 15, which take into account inclination discovery bias in direct imaging surveys caused by a hard IWA. 5. We do not find evidence that our sample is biased against pole-on orientations as might be expected for stars with solar-like dynamos and starspots that are confined to equatorial regions. 6. Relative age biases could impact stellar obliquities for systems with giant planets and brown dwarf companions, but only if obliquity evolution occurs on relatively short timescales (1 Gyr).
Altogether, it is difficult to confidently establish how these biases might be combining to affect individual and cumulative results for Δi distributions. This certainly warrants future consideration, especially in conjunction with a better understanding of other unknown astrophysical obliquity trends that we are ignoring. These could include stellar obliquities as a function of system age, companion mass, separation, eccentricity, or multiplicity. We did not consider these in this study because the sample size and individual constraints are generally modest, but in the future these relationships can be explored in more detail with larger samples and more precise constraints on Δi.

Spin-Spin Opportunities
Most companions to the 53 host stars with stellar inclinations do not have constraints on the inclinations of their orbital planes because the physical separations in these systems are prohibitively large. However, this sample may be beneficial to explore the statistical distributions of spin-spin orientations between the host stars (i * ) and companions (i p ). A handful of obliquities have now been measured for individual brown dwarf companions and giant planets by combining rotation periods, radii, and projected rotational velocities (e.g., Bowler et al. 2020b;Bryan et al. 2020aBryan et al. , 2021Zhou et al. 2020;Palma-Bifani et al. 2023), but larger samples are needed to assess the prevalence of stellar-companion spin-spin alignment in order to interpret this in the context of formation scenarios (e.g., Jennings & Chiang 2021). Companions to host stars with existing inclination constraints from this study may be promising targets for follow-up photometric monitoring and high-resolution spectroscopy in the future.

Summary
We have carried out a homogeneous analysis of the rotation periods, inclinations, and obliquities of stars hosting directly imaged brown dwarf and giant planet companions. Below we summarize the main conclusions and interpretation from this work.
1. TESS light curves for substellar host stars were uniformly examined in search of rotational modulations. Our search was isolated to spectral types F5 to avoid confusion with stellar pulsations. Supplementing these rotation periods with values from the literature yields 62 stars with period measurements, 53 of which also have * v i sin values from our own new observations and previous determinations. We derived inclination posteriors using the Bayesian framework from Masuda & Winn (2020) along with general analytical expressions for inclination posteriors given P rot , R * , and * v i sin values in Appendix A. 2. Twenty-one systems with stellar inclinations also have constraints on the inclination of the companion's orbital plane. Together these provide a measurement of Δi, the minimum value of the true spin-orbit obliquity.
Companions in this subsample are predominantly brown dwarfs located between 10 and 100 au with a few companions extending out to 250 au. 3. Spin-orbit misalignments for systems with substellar companions appear to be common, as revealed by Δi distributions. Eleven of the 21 systems in our sample (52 -+ 11 10 %) have Δi values beyond 10°with 80% confidence, indicating that misalignments are prevalent among brown dwarf companions. Depending on the orientation of the star in the sky plane and the longitude of ascending node of the orbital plane, most companions can be on prograde, polar, or retrograde orbits. This high misalignment fraction contrasts with partial or full constraints on stellar obliquities for systems with directly imaged planets, which show an emerging pattern of angular momentum alignment between the star, planet orbital plane, and, when present, the disk. We also find signs of a preference for spin-orbit alignment for companions within 20 au. Brown dwarfs and companions at wide separations show a greater diversity of rotational and orbital angular momentum architectures. 4. When the same systems are compared with stars with random orientations, stars with rotation axes aligned with orbital angular momentum vectors, and a mixture of both scenarios, our hierarchical Bayesian analysis of the underlying distributions of Δi for the observed sample is most consistent with the even mixture and disfavors the aligned case. However, the modest sample size, generally broad Δi distributions, potential for sample biases, and minimal information about the true obliquity angle ψ limits more nuanced conclusions about how obliquities might depend on age, companion mass, and separation. 5. Altogether, these results are consistent with predictions from hydrodynamic simulations that spin-orbit misalignments should be a common outcome for binary companions that form from a fragmenting molecular cloud core (e.g., Bate et al. 2010;Offner et al. 2016).
Although these results could also be explained by dynamical processing of systems with primordial alignments, this would require another massive object to interact with in these systems. Large high-contrast imaging surveys have failed to uncover a population of multiple brown dwarf companions orbiting stars in nonhierarchical configurations. Some misaligned systems may also result from the capture of isolated or ejected giant planets and brown dwarfs in young dense clusters. 6. We have evaluated potential biases that could be impacting the stellar and orbital inclinations in our sample. There is evidence for a higher number of systems on edge-on orbits than would be expected from chance, which may be a result of stars with differential rotation or slow rotation in the sample. More detailed modeling of the impact of these effects is warranted. 7. As part of this analysis we highlight a bias in the orbital inclination distribution of companions found with highcontrast imaging in the presence of an IWA-which may be a coronagraph or residuals from primary star subtraction. This "inclination discovery bias" is most severe for objects close to the IWA and varies with semimajor axis and eccentricity. Its effect is to disfavor edge-on orientations and deviates substantially from an isotropic distribution in which ( ) µ  i i sin o o . The distribution of orbital inclinations for substellar companions is consistent with having been shaped by the presence of this effect. We formulate a parametric approximation for this bias, which may be useful as a more realistic Bayesian prior for orbit fits to new discoveries. Stellar obliquities provide a rich resource to probe the formation and dynamical histories of planets and brown dwarf companions. Results from this study reinforce an emerging distinction between imaged giant planets and brown dwarf companions in terms of their orbits , demographics (Nielsen et al. 2019;Vigan et al. 2021), and mass functions (Wagner et al. 2019): brown dwarf companions are most consistent with formation in a stellar-like pathway from fragmenting turbulent clouds than within aligned disks. In the future, we expect that new discoveries from dynamically informed surveys (e.g., Franson et al. 2023;Hinkley et al. 2022), continued astrometric orbit monitoring of known systems (e.g., Zurlo et al. 2022), precision astrometry from VLTI/GRAVITY Wang et al. 2021), and RVs of companions from instruments such as the Keck Planet Imager and Characterizer (Mawet et al. 2016) will be especially fruitful to build a larger sample of systems well-characterized spin-orbit constraints. In addition, spectro-interferometry capable of spatially resolving offsets from stellar rotation will enable true obliquity measurements for imaged planets orbiting early-type stars . Ultimately, refined orbital inclinations and projected rotational velocities will facilitate more precise population-level obliquity constraints to carry out direct tests for differences in the obliquity distributions between long-period brown dwarf companions and directly imaged planets. has made use of the VizieR catalog access tool, CDS, Strasbourg, France (DOI:10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23.

Appendix A Analytical Formalism to Infer Stellar Inclinations
A.1. Measurements of v sin i * Masuda & Winn (2020) provided a thorough guide to infer the posterior probability distribution function of stellar inclination (i * ) given measurements or estimates of a rotation period (P rot ), stellar radius (R * ), and projected rotational velocity ( * v i sin ) following a Bayesian formalism. The equatorial velocity v eq is simply 2π R * /P rot . If i * , R * , P rot , and However, one of the important insights recognized by Masuda & Winn (2020) is that v eq and * v i sin are not independent, so accurately determining the probability distribution of i * must take into account the correlation between these parameters.
Under the reasonable assumption of independence between the data sets used to measure v eq and * v i sin (so that their likelihood functions are separable), and assuming v eq and i * are similarly independent (so that their priors are separable), the probability distribution for * i cos given a set of data d is: Here  v eq is the likelihood function for the equatorial velocity, is the Bayesian prior on the cosine of the inclinations-equal to 1 for isotropic orbits-and ( )  v v eq eq is the prior on the equatorial velocity. If both priors are constant, then sin eq 2 eq eq Since * v i sin is typically measured from high-resolution spectroscopy, it is usually reported with Gaussian uncertainties so its likelihood function will follow a normal distribution in eq 1 cos 2 sin sin 2 If the equatorial velocity is taken from measurements (or estimates) of R * and P rot , then the likelihood function for v eq may not be straightforward to express analytically. While the stellar radius and rotation period may be normally distributed such that rot , the ratio of two independent normally distributed random variables of the form 2πR * /P rot itself does not generally follow a normal distribution. An analytical expression has been derived but it is not succinct, even in the case when both random variables are uncorrelated (Hinkley 1969;Oliveira et al. 2015).
However, there are certain cases when this normal ratio distribution can be approximated with a Gaussian distribution. For independent random variables ( ) m s  X , 2 , where β = μ x /μ y , δ x = σ x /μ x , and δ y = σ y /μ y . Smaller values of δ x and δ y produce better agreement between the normal approximation and true distribution function. In general, the normal approximation is good if δ y is small (0.2). For the application to v eq , this corresponds to a rotation period measurement at the 20% level, which is satisfied for all targets we consider in this study.
Under these assumptions, the likelihood function for v eq becomes Substituting in the values of a, b, and c gives the following expression for the posterior probability of the cosine of the stellar inclination: The uncertainty in the equatorial velocity v eq is given in Equation (A5). Finally, it can be useful to represent the stellar inclination constraint in terms of inclination i * itself. This can be achieved using a variable transformation of the form f (i * ) di * = g(y) dy.
Here y is a function of i * , ( ) = = * * y y i i cos , and g(y) is the  Under the assumptions outlined above, Equations (A11) and (A12) provide the full (unnormalized) expressions for the distribution function of a stellar inclination and its cosine given measurements of * v i sin , P rot , and R * . In Figure 18 we test these relations with the same examples used in Masuda & Winn (2020). Here we assume R * = 1.0 ± 0.1 R e and P rot = 5.059 ± 0.88 days, which is tailored to give v eq = 10 ± 2 km s −1 . Compared to their Figure 1, our results are in good agreement for two cases of projected rotational velocities: * v i sin = 9.8 ± 2.0 km s −1 and * v i sin = 4.0 ± 2.0 km s −1 .

A.2. Upper Limits on v sin i *
In some cases, only an upper limit on * v i sin is available. This can occur if absorption lines are not fully resolved by an instrument or if a star is rotating slowly. Here we derive analytical expressions for this scenario in a similar fashion as was carried out for normally distributed measurements of For an upper limit l on the projected rotational velocity, its likelihood function is simply a uniform distribution from 0 km s −1 to l km s −1 : 0, sin km s −1 . This is equivalent to placing limits on the definite integral over * v i sin from 0 to l. The likelihood function for the equatorial velocity, , is approximately normally distributed following Equation (A4) if the rotation period is measured to a precision of 20%. Adopting an isotropic prior for ( ) *  i cos and a uniform prior for v eq , Equation ( Therefore, for an upper limit on * v i sin of l km s −1 , the posterior distribution for * i cos is simply proportional to the sum of two error functions. Several examples using Equations (A16) and (A17) are shown in Figure 19 following Masuda & Winn (2020). As in Appendix A.1, we use R * = 1.0 ± 0.1 R e and P rot = 5.059 ± 0.88 days for this exercise. For a star with v eq = 10 ± 2 km s −1 , we vary the upper limit of a hypothetical * v i sin measurement from 2.5 to 30 km s −1 . As the upper limit moves beyond the equatorial rotational velocity, the posterior distribution of * i cos tends to unity, and for i * it tends to an isotropic distribution.

Appendix B Notes on Individual Stars
1RXS J034231.8+121622: This young mid-M dwarf is a likely member of the Hyades cluster (Kuzuhara et al. 2022). Its brown dwarf companion was found by Bowler et al. (2015a), and the latest orbit for the pair is determined in Bowler et al. (2020a). Our 2MASS J01033563-5515561: This low-mass hierarchical triple system comprises a close visual binary (AB) with an unresolved spectral type of M5-M6 and a wide 12-14 M Jup accreting L-dwarf companion, 2MASS J01033563-5515561 b (Delorme et al. 2013;Eriksson et al. 2020;Betti et al. 2022). The light curve from two TESS sectors shows a strong periodic signal at 0.1664 ± 0.0008 days (3.994 ± 0.019 hr), which is consistent with fast rotation expected as a lowmass member of the ≈45 Myr Tuc-Hor young moving group (Kraus et al. 2014b). It is unclear which member of the binary is producing the observed modulations, but their near-equal flux ratios imply that they share similar physical properties.
To estimate their radii, we adopt an effective temperature of 2980 ± 75 K and J-band bolometric correction of 1.99 mag for young stars from Herczeg & Hillenbrand (2015) for 2MASS J01033563-5515561 A, which corresponds to young M5 stars. Assuming the B component is slightly cooler with a spectral type of M5.5, this corresponds to T eff = 2920 ± 75 K and BC J = 2.01 mag. The absolute J-band magnitude of 2MASS J01033563-5515561 A is M J = 7.36 ± 0.05 mag (Delorme et al. 2013), which corresponds to a bolometric magnitude of M bol = 9.35 ± 0.05 mag, a luminosity of  L L log bol = -1.84 ± 0.02 dex, and a radius of R * = 0.45 ± 0.03 R e . The properties are similar for B: M J = 7.56 ± 0.05 mag (Delorme et al. 2013), M bol = 9.57 ± 0.05 mag,  L L log bol = -1.93 ± 0.02 dex, and a radius of R * = 0.43 ± 0.02 R e .
For the inclination analysis, we assume the properties of A; although given the consistent radii, the results would be similar for B. We find a stellar inclination of i * = 8.9  -+ 1.3 0.9 . In the future, this measurement may be compared with the orbital angular momentum of the AB binary and, in principle, both the spin and orbit of the companion.
2MASS J01225093-2439505: This M3.5 star is a member of the AB Dor young moving group and hosts an imaged L4 substellar companion with a mass near the deuterium-burning limit (Bowler et al. 2013;Hinkley et al. 2015).
2MASS J01225093-2439505 has an inferred equatorial rotational velocity of 12.4 ± 1.0 km s −1 , which is inconsistent with the projected rotational velocity of 18.1 ± 0.5 km s −1 . The equatorial rotational velocity is based on a radius of 0.37 ± 0.03 R e from Stassun et al. (2019) and a rotation period of 1.493 ± 0.035 days from our analysis of the TESS light curve spanning two sectors. Note that this rotation period is nearly identical to the value of 1.49 ± 0.02 days determined by Bryan et al. (2020a) from a single TESS sector. The inconsistency between v eq and * v i sin implies that either the radius is too small, the rotation period is overestimated, the * v i sin value is overestimated, or some combination of these. The TESS light curve covers about 26 period cycles and shows regular single-mode modulations in Sector 3 followed by double-peaked modulations in Sector 30 with different peak amplitudes. Because of this double-peaked structure, which seems to trace two dominant groups of starspots at different longitudes, it is unlikely that the signal itself is an alias of a true underlying signal (as was discussed for LP 261-75 in Section 3.3). However, it is possible that mid-latitude starspots coupled with strong differential rotation could cause light-curve modulations that are longer than the equatorial period, which might account for the discrepancy between v eq and * v i sin . It is also possible that the radius is underestimated. If the radius is the sole origin of the discrepancy between v eq and * v i sin , a size of ≈0.53 R e would be needed to reconcile the two values. We recompute the radius using the apparent 2MASS J-band magnitude of 10.084 ± 0.028 mag, a J-band bolometric correction of 1.885 mag (midway between M3 and M4 from Herczeg & Hillenbrand 2015), an effective temperature of 3300 K (Herczeg & Hillenbrand 2015), and a Gaia DR3-based distance of 33.74 ± 0.03 pc (Gaia Collaboration et al. 2022). This yields a bolometric magnitude of M bol = 9.33 ± 0.05 mag, a luminosity of  L L log bol = -1.83 ± 0.02 dex, and a radius of R * = 0.37 ± 0.02 R e -in good agreement with the TIC value from Stassun et al. (2019).
Finally, it is also possible that the * v i sin value of 2MASS J01225093-2439505 may be overestimated. Our adopted value of 18.1 ± 0.5 km s −1 is a weighted average of 14.2 ± 3.2 km s −1 from Malo et al. (2014) and 18.2 ± 0.5 km s −1 from Bryan et al. (2020a). Additional measurements of rotational broadening for this star would help clarify the nature of the discrepancy between v eq and * v i sin . FU Tau: FU Tau A is a young accreting late-M dwarf in the Taurus star-forming complex (Jones & Herbig 1979), which may be driving a large molecular outflow (Monin et al. 2013; although see Wu et al. 2020 for an alternate interpretation of the line emission). Its wide companion, FU Tau B, was discovered by Luhman et al. (2009) and is expected to have a mass near the deuterium-burning boundary.
The inferred mass and radius of FU Tau A is highly sensitive to the measured spectral type, extinction, effective temperature, and system age; values of 0.05-0.2 M e and 0.55-1.8 R e have been previously estimated (Luhman et al. 2009;Stelzer et al. 2010Stelzer et al. , 2013Bulger et al. 2014 2022) is 7.835 ± 0.075 mas, which corresponds to a distance estimate of 126.4 ± 1.0 pc (Bailer- Jones et al. 2021). This is about 10% closer than the 140 pc value that was often used in the past for Taurus members that lacked direct distance measurements, including FU Tau. We therefore re-derive the fundamental parameters of FU Tau A following a similar procedure as in Stelzer et al. (2013) using an (extinction-corrected) J-band bolometric correction and the revised distance.
Here we use a J-band bolometric correction of BC J = 2.045 mag from Herczeg & Hillenbrand (2015), which is intermediate between young M6 and M7 stars and assumes a spectral type of M6.5 (Stelzer et al. 2013;Herczeg & Hillenbrand 2014) for FU Tau A. The optical extinction of A V = 0.5 ± 0.5 mag from Stelzer et al. (2013) corresponds to a Jband extinction of A J = 0.17 ± 0.11 mag based on infrared interstellar reddening law in Tokunaga (2000). This gives an absolute J-band magnitude of M J = 5.10 ± 0.12 mag for FU Tau A, a bolometric magnitude of M bol = 7.15 ± 0.12 mag, a bolometric luminosity of  L L log bol = -0.96 ± 0.05 dex, and a radius of R * = 1.4 ± 0.1 R e . These values rely on an effective temperature of T eff = 2815 ± 75 K for a young M6.5 star (Herczeg & Hillenbrand 2015).
Based on the periodicity at 3.93 days we find from the TESS light curve, we infer an inclination of i * = 75  -+ 5 14 . The millimeter disk around FU Tau was detected by Wu et al. (2020) with the Atacama Large Millimeter/submillimeter Array (ALMA) but was not spatially resolved in that data. In the future, the inclination of FU Tau A's disk and FU Tau B's spin orientation may be compared with the inclination measured here.
FW Tau: FW Tau is a binary pair of young low-mass stars in the Taurus star-forming complex. The pair has equal flux ratio, implying similar physical properties, and has shown significant orbital motion since 1989 (Schaefer et al. 2014). Bowler et al. (2014) found a spectral type of M6 ± 1 and modest extinction of A V = 0.4 -+ 0.4 1.3 mag. A faint companion was found with HST imaging by White & Ghez (2001) and confirmed to be comoving by Kraus et al. (2014a); however, follow-up studies suggest FW Tau C may be a brown dwarf or a low-mass star with an edge-on disk (Bowler et al. 2014;Caceres et al. 2015;Kraus et al. 2015;Mora et al. 2020).
We compute the radii of FW Tau A and B by assuming each component is identical in luminosity and effective temperature. No parallax is available in Gaia DR3 for FW Tau, presumably because of its binary nature, so we adopt a characteristic distance of 140 ± 10 pc for Taurus (Kenyon et al. 1994). Using a bolometric correction of 2.03 mag and T eff = 2860 ± 75 K from Herczeg & Hillenbrand (2015), an optical extinction of A V = 0.4 mag (implying a J-band extinction of A J = 0.11 mag following Tokunaga 2000), and an apparent magnitude decomposed into equal-brightness individual magnitudes of J = 11.1 ± 0.05 mag, we find M J = 5.25 ± 0.16 mag, M bol = 7.28 ± 0.16 mag,  L L log bol = -1.01 ± 0.06 dex, and a radius of R * = 1.28 ± 0.11 R e .
Our stellar rotation period of 0.91 day agrees well with the value of 0.906 day found by Chen et al. (2020) using g-band photometry from the Zwicky Transient Factory. They also report an r-band period of 0.475 day, which is probably an alias of the true rotation period. Combined with the * v i sin of 49.2 ± 0.7 km s −1 from Kounkel et al. (2019), the resulting posterior stellar inclination distribution peaks at i * = 44°, and the 68% credible interval spans 39°-49°. Although the orbit of AB is not yet well constrained (Schaefer et al. 2014), its inclination and the properties of FW Tau C such as its disk inclination and spin inclination can eventually be compared with this stellar inclination to shed light on the formation of this hierarchical multiple system.
Gl 229: This nearby (5.7614 ± 0.0006 pc; Gaia Collaboration et al. 2022) M1 star hosts the first T dwarf companion to have been directly imaged (Nakajima et al. 1994;Oppenheimer et al. 1995;Oppenheimer 2014). With over 25 years of relative and absolute astrometry, the orbit and dynamical mass of Gl 229 B have now been exquisitely measured (Brandt et al. 2020(Brandt et al. , 2021bFeng et al. 2022).
The equatorial velocity of 1.0 ± 0.2 km s −1 that we infer for Gl 229 is slightly (but significantly) smaller than its * v i sin value of 2.8 ± 0.4 km s −1 , which is based on a compilation of projected rotational velocities from the literature (see Table 4). Most of the literature values do not include estimates of the measurement uncertainties, so these are grouped into a single value following the procedure outlined in Section 4.3. Note that our adopted value is dominated by the measurement of 2.6 ± 0.4 km s −1 from Hojjatpanah et al. (2019), which is more heavily weighted because of its small reported uncertainties.
Most of these published * v i sin measurements indicate that Gl 229 is a slow rotator, which is consistent with its old age (2-6 Gyr; Brandt et al. 2020), lack of activity signatures, and long rotation period of 27 days (Suarez Mascareno et al. 2016). Several * v i sin measurements fall below 2 km s −1 (Glebocki & Gnacinski 2005;Reiners 2007;Reiners et al. 2022), which at this level becomes very challenging to determine because it can be difficult to separate contributions to line broadening from rotation, microturbulence, and macroturbulence. The resulting obliquity constraint for Gl 229 may be suffering from a systematic bias because of an overestimated * v i sin value. The disagreement could also stem from the radius estimate or the photometric rotation period measurement. The radius of 0.55 ± 0.06 R e we adopt is similar to several other estimates such as 0.58 ± 0.05 R e from Houdebine et al. (2016) and 0.52 ± 0.02 R e from Schweitzer et al. (2019). An independent rotation period measurement would also be valuable for this star to compare with the 27 days value found by Suarez Mascareno et al. (2016).
SR 12 AB: SR 12 AB (2MASS J16271951-2441403, ROXs 21) is a close resolved binary (Simon et al. 1987) and a wellestablished member of the ρ Ophiuchus cloud complex (e.g., Bouvier & Appenzeller 1992). Kuzuhara et al. (2011) discovered a widely separated companion, SR 12 c, with a mass near the deuterium-burning limit. SR 12 c was later found to have an accretion disk (Santamaría- Miranda et al. 2017;Martinez & Kraus 2022) with a dust mass of about one lunar mass based on a submillimeter detection with ALMA (Wu et al. 2022).
Here we derive the radius of SR 12 A using the Stefan-Boltzmann relation, but results are expected to be similar if the photometric rotation originates from SR 12 B. Schaefer et al.
(2018) measured a flux ratio in the narrow-band Keck/NIRC2 J cont filter of 0.951 ± 0.006, or 0.055 ± 0.006 mag. The apparent integrated-light J-band magnitude of SR 12 is 9.424 ± 0.023 mag (Cutri et al. 2003), which decomposes into apparent magnitudes of J A = 10.149 ± 0.023 mag and J B = 10.204 ± 0.023 mag. A parallax of 8.9034 ± 0.4288 mas (d = 112 ± 5 pc; Bailer- Jones et al. 2018) is listed in Gaia DR2, but no parallax is present in DR3, likely because orbital motion of the binary is evident in the updated astrometry (the excess astrometric noise parameter is 7.121 mas but no RUWE value is listed). We therefore instead use the Gaia-based distance of 138.4 ± 2.6 pc inferred by Ortiz-León et al. (2018) for other members of the L1688 cloud. Adopting a spectral type of M0 from Wilking et al. (2005) implies an effective temperature of T eff = 3900 ± 125 K (Herczeg & Hillenbrand 2015). Here we assume the J-band extinction of A J = 0.5 ± 0.2 mag from Kuzuhara et al. (2011). This implies an extinction-corrected absolute magnitude of M J = 3.95 ± 0.21 mag, a bolometric magnitude of M bol = 5.61 ± 0.21 mag, a bolometric luminosity of  L L log bol = -0.34 ± 0.08 dex, and a radius of R * = 1.48 ± 0.17 R e .
For this study we adopt the rotation period of 3.918 days found by Rebull et al. (2018) from K2, which is similar to the period of 3.927 days from Kiraga (2012) and 3.516 days from Grankin et al. (2008). Note that Jayasinghe et al. (2018) reported a period of 7.849 days using light curves from the All-Sky Automated Survey for Supernovae (ASAS-SN); this is about twice the value from Rebull et al. (2018) so one is probably an alias of the true period. The cadence of ASAS-SN is every few days, whereas K2 is much finer (≈30 minutes) and therefore should be more reliable.
Combining our radius estimate with the rotation period and adopted * v i sin value of 27 ± 11 km s −1 ( Table 4) yields a broad inclination constraint that only slightly deviates from the * i sin prior. Much of this is caused by the large uncertainty in projected rotational velocity. Orbital motion of AB has been detected but is not yet well determined (Schaefer et al. 2018); in the future, the binary orbital plane and the spin inclination of the wide companion (Bryan et al. 2020b) can be compared with this host inclination to test for alignment.
TWA 5 Aab: TWA 5Aab is an M2 member of the ≈10 Myr TW Hydrae association. The brown dwarf companion TWA 5 B was found by Webb et al. (1999) and Lowrance et al. (1999), and the host was later resolved into a close visual binary by Macintosh et al. (2001). Köhler et al. (2013) measured a dynamical mass of 0.9 ± 0.1 M e for the total mass of TWA 5Aab and a mass ratio of 1.3, implying individual masses of about 0.4 M e and 0.5 M e of the Ab and Aa components, respectively-albeit with large uncertainties. The radius of TWA 5 is listed as 0.40 ± 0.07 R e in Stassun et al. (2018), but this value differs substantially from theoretical expectations, perhaps because TWA 5 is a close binary. For example, solar-metallicity MIST models (Choi et al. 2016;Dotter 2016) predict a radius of 0.88 R e for a 0.5 M e star at 10 Myr. We instead adopt the weighted mean of radius estimates from Stassun et al. (2019) for two other coeval stars in the TWA association with the same spectral type of M2, TWA 7 (0.917 ± 0.119 R e ) and TWA 10 (0.860 ± 0.114 R e ). Following the discussion in Section 4.6, we also add a 7% uncertainty in quadrature with the original error estimates. This yields 0.89 ± 0.10 R e .

Appendix C TESS Light Curves
Figures 20-27 show TESS light curves for host stars in our sample with periodic variations consistent with rotational modulation. Full light curves are shown in the left panels with breaks between noncontiguous sectors. The middle panels depict the Generalized Lomb-Scargle periodogram with the most prominent peak highlighted. The corresponding phased light curves and binned phased curves are displayed in the right panels. Starspot evolution is apparent for many targets.
Individual TESS sectors and rotation period measurements can be found in Table 2. A full description of the light-curve extraction and reduction is summarized in Section 3.3. The observations used in this analysis can be accessed from MAST (Caldwell et al. 2020a; Huang 2020).    Table 4. When available, N independent measurements for the same star are combined into a single adopted value by computing a weighted mean and standard deviation: In some cases, projected rotational velocities are reported in the literature at precisions below 1%. These would dominate over other measurements with more typical uncertainties (5%-20%) when computing a weighted mean. Because of the difficulty in separating the effects of microturbulence, macroturbulence, and rotational broadening, we conservatively In many instances, * v i sin values are reported without associated uncertainties. In these cases, when more than one are identified, we merge these into a single measurement by computing a robust mean and standard deviation following Beers et al. (1990).