Disc–star alignment I: pre-main-sequence stellar parameters and the statistical alignment between discs and stellar rotation

ABSTRACT

Astronomers generally assume planet-forming discs are aligned with the rotation of their host star. However, recent observations have shown evidence of warping in protoplanetary discs. One can measure the statistical alignment between the inclination angles of the disc and stellar spin using the projected rotational velocity, radius, and rotation period of the star and interferometric measurements of the protoplanetary disc. Such work is challenging due to the difficulty in measuring the properties of young stars and biases in methods to combine them for population studies. Here, we provide an overview of the required observables, realistic uncertainties, and complications when using them to constrain the orientation of the system. We show in several tests that we are able to constrain the uncertainties on the necessary stellar parameters to better than 5 per cent in most cases. We show that by using a hierarchical Bayesian model, we can account for many of the systematic effects (e.g. biases in measured stellar and disc orientations) by fitting for the alignments of each system simultaneously. We demonstrate our hierarchical model on a realistic synthetic sample and verify that we can recover our input alignment distribution to |$\lesssim 5^\circ$| with a modest (⁠|$\simeq$|30 star) sample. As the sample of systems with disc inclinations grows, future studies can improve upon our approach with a 3D treatment of misalignment and better handling of non-Gaussian errors.

1 INTRODUCTION

Discs of gas and dust surrounding newly formed stars (called ‘protoplanetary’ or ‘planet-forming’ discs) are natural consequences of angular momentum conservation as molecular clouds collapse to form protostars. Young stars are initially embedded in gas-rich envelopes, but after |$\simeq$|1 Myr they become optically visible as the dust and gas settle into a disc and evolve through viscous accretion. The protoplanetary discs around such pre-main-sequence (PMS) stars are thought to form planets within 5–10 Myr, after which their gaseous disc material dissipates, first into dusty debris discs, then into mature planetary systems (Williams & Cieza 2011; Morbidelli et al. 2012).

The existence of protoplanetary discs was implied by the shape of the Solar System; the planets appeared to occupy a flat plane, orbiting in-line with the rotation of the Sun. In fact, the rotational axis of the Sun is only inclined |$6^{\circ }$| relative to the average angular momentum of the eight major planets (Beck & Giles 2005), and the spread in inclinations is no more than |$\sim 7^{\circ }$|⁠. The observed geometry of the Solar System inspired the Kant–Laplace nebular hypothesis, which led to the widely accepted idea that protoplanetary discs should be aligned with the rotation of their newly formed host stars.

Observations of the Rossiter–McLaughlin effect during planetary transits have revealed a significant population of planets with orbits misaligned from their host star’s spin axis (e.g. Winn et al. 2010; Albrecht, Dawson & Winn 2022). This seemingly contradicts the basic planet-formation model, so early explanations often focused on later-stage many-body interactions, such as Kozai–Lidov (Wu & Murray 2003), planet–planet scattering (Ford & Rasio 2008; Naoz et al. 2011), or angular momentum transport within the host star (Rogers, Lin & Lau 2012). Other work has shown that planetary misalignment can arise from early interactions with the protoplanetary disc (e.g. Petrovich et al. 2020), particularly if the disc is misaligned (e.g. Batygin 2012).

A protoplanetary disc can be quickly misaligned from its host due to the influence of a wide binary companion, a dense stellar environment (stellar flybys), or interactions with the surrounding star-forming cloud (e.g. Bate, Lodato & Pringle 2010; Fielding et al. 2015; Takaishi, Tsukamoto & Suto 2020; Kuffmeier et al. 2021). The inner disc may stay bound to the orientation of the host star, but hydrodynamical simulations suggest the outer disc (⁠|$\gtrsim$|10 au) can end up with a near-random orientation (e.g. Bate 2018).

We define the disc-star alignment angle (⁠|$\alpha$|⁠) as the difference between the inclination angles of the disc orbital axis and the stellar spin axis. This is not to be confused with the sky-projected spin-orbit angle in transiting exoplanet studies (i.e. the sky-projected obliquity |$\lambda$|⁠; see e.g. Fabrycky & Winn 2009; Dong & Foreman-Mackey 2023) which is perpendicular to |$\alpha$|⁠. We show a diagram of a misaligned protoplanetary disc in Fig. 1.

Recent studies using observations from the Atacama Large Millimeter/submillimeter Array (ALMA) have found misaligned protoplanetary discs that are pole-on (⁠|$\sim 0^{\circ }$|⁠; e.g. Kennedy et al. 2019; Ansdell et al. 2020) or have potentially retrograde (e.g. Kraus 2020) orbits. High-resolution imaging has also revealed shadows in a significant fraction of the so-called transition discs that can be best explained by a misalignment between inner and outer discs (e.g. Avenhaus et al. 2014; Casassus et al. 2018). How often such misalignments occur is unclear.

There have only been a few statistical surveys of disc–star alignment, primarily focusing on small (15–20) samples of debris (Watson et al. 2011; Greaves et al. 2014) or protoplanetary (Davies 2019) discs. These have generally found that most systems are aligned or weakly misaligned (⁠|$< 30^{\circ }$|⁠). One more recent study found 6 of 31 debris discs were significantly misaligned (Hurt & MacGregor 2023).

These studies often used previously published stellar parameters which can be unreliable for cool (e.g. Mann et al. 2015; Newton et al. 2016) and PMS stars (e.g. Kraus et al. 2015; Rizzuto et al. 2016). They were also limited in sample size and stuck to a basic comparison between the inferred stellar and disc inclinations without accounting for the statistical biases and co-dependence of these inputs. While such studies are valuable, these limitations make it difficult to draw statistical conclusions.

This work generally relies on computing the stellar inclination (⁠|$i_{\star }$|⁠) from the combined rotation period (⁠|$P_{\textnormal {rot}}$|⁠), stellar radius (⁠|$R_{\star }$|⁠), and rotational spectral broadening (⁠|$v\sin i_{\star }$|⁠). This is effectively comparing the line-of-sight broadening to the equatorial velocity (⁠|$v_{\textnormal {eq}}$|⁠) as in Campbell & Garrison (1985):

This method has been similarly used for studying the alignment between the orbits of (transiting) planets and the rotation of their stellar hosts (e.g. Mann et al. 2022; Wood et al. 2023a). However, the approach is far more complicated than implied by equation (1) due to the dependent and non-Gaussian probability distributions and (large) measurement uncertainties (see Morton & Winn 2014; Masuda & Winn 2020, for further discussion). Further, young disc-bearing stars tend to have faster rotation (and hence more easily measured |$v\sin i_{\star }$| and |$P_{\textnormal {rot}}$|⁠), but the presence of the protoplanetary disc and the difficulties obtaining accurate radii of PMS stars make this more challenging than studying mature planet hosts.

In this paper, we explore the challenges associated with estimating precise stellar parameters and a means of combining them with the disc inclination (⁠|$i_{\textnormal {d}}$|⁠) to study the disc–star alignment distribution. Our goal is to provide an approach to measure the statistical projected alignment between stars and their protoplanetary discs with realistic uncertainties. In Section 2, we outline our methodology for estimating the disc–star alignment angle. We test part of our basic methodology to measure stellar parameters for a sample of PMS stars lacking protoplanetary discs, which we discuss in Section 3. In Section 4, we test the full method on a synthetic population of disc-bearing stars. Finally, in Section 5, we summarize the results of the various tests and offer concluding remarks on how our methods can be applied to a real disc-bearing star population.

2 METHODS

Most simply, we need |$i_{\star }$| and |$i_{\textnormal {d}}$| for a sample of targets. For |$i_{\star }$|⁠, we need |$v\sin i_{\star }$|⁠, |$R_{\star }$|⁠, and |$P_{\textnormal {rot}}$|⁠, and realistic uncertainties, as well as a framework to handle complications when applying equation (1). We describe each of these parameters in the subsections below, followed by a description on how one can combine the sample within a hierarchical Bayesian framework, as well as a discussion of complications from missing spatial information and multistar systems.

2.1 Projected stellar rotation velocity

|$v\sin i_{\star }$| is most commonly derived from the broadening of spectral lines in high-resolution spectra. |$v\sin i_{\star }$| is observationally derived as a single parameter despite containing information on both the rotation and inclination of the star. A particular problem for this program, |$v\sin i_{\star }$| can be biased in the youngest stars due to degeneracies with pressure and magnetic broadening.

We adopted the method used by Kesseli et al. (2018), which is similar to methods used in some previous studies (e.g. West & Basri 2009; Muirhead et al. 2013; Reiners et al. 2018). In brief, we compared a given target spectrum to a slowly rotating template spectrum of a star with the same spectral type as the target. The templates have negligible |$v\sin i_{\star }$| (⁠|$\ll$|2 km s|$^{-1}$|⁠) compared to the target spectrum and instrumental broadening. We artificially broadened the template with a grid of |$v\sin i_{\star }$| values (from 2 to 62 km s|$^{-1}$|⁠), and cross-correlated each broadened template with the original, unbroadened template. Then, we cross-correlated the target spectrum with the unbroadened template spectrum. Finally, we measured the full-width-at-half-maximum (FWHM) of the target’s cross-correlation function (CCF), and we linearly interpolated the FWHM onto the grid of FWHM values from the broadened template CCFs to determine the |$v\sin i_{\star }$| value. In Fig. 2, we show an example of this method.

Artificial broadening requires a limb-darkening coefficient, which we calculated using the Python Limb-Darkening Toolkit (LDTk; Parviainen & Aigrain 2015), which uses models from Husser et al. (2013), and stellar parameters which we describe in Section 2.2. Varying limb-darkening for a range of stars suggests uncertainties in limb-darkening parameters have a negligible impact on the final |$v\sin i_{\star }$| when compared to other sources of uncertainty.

For our tests, we opted to use K-band spectra from the Immersion Grating Infrared Spectrograph (IGRINS; Park et al. 2014), as were used in Kesseli et al. (2018). IGRINS spectra are high-resolution (⁠|$R\simeq 45\, 000$|⁠), more than sufficient for the expected |$v\sin i_{\star }$| seen in disc-bearing stars (10’s of km s|$^{-1}$|⁠). IGRINS also covers the near-infrared (NIR) CO bands, which have relatively equally spaced and well-separated lines which are resistant to pressure and magnetic broadening that can be degenerate with |$v\sin i_{\star }$| (see also Lavail, Kochukhov & Hussain 2019; López-Valdivia et al. 2021). NIR spectrographs are also favourable because of high extinction in many star-forming regions.

A downside of using the CO bands is that the lines are weak in warmer (⁠|$\gtrsim$|5000 K) stars. At ages where discs are present (0–5 Myr Mamajek 2009), only stars |$\gtrsim$|1.5 |$M_{\star }$| are this warm. For a typical initial mass function, this is only |$\simeq$|6 per cent of stars; these targets are also more challenging in terms of assigning stellar parameters.

IGRINS has archival spectra from which we can draw slow-rotating templates.¹ We also tested using PHOENIX (Allard et al. 2013) model spectra, and found they often yield lower |$v\sin i_{\star }$| values (by 0.1–1 km s|$^{-1}$|⁠) and larger order-to-order uncertainties. This may be due to inaccuracies in the instrumental broadening and/or effects of missing astrophysics in the models. While the offset is small, because it is systematic, we opted for empirical templates.

For each target, we estimated |$v\sin i_{\star }$| from echelle orders [4, 5, 6, 11, 12, 13, 14] (2.09–2.36 |$\mu$|m), sometimes excluding an order due to a strong emission line and/or no significant absorption lines. We combined the resulting |$v\sin i_{\star }$| values using a weighted mean where the weights were chosen as the maximum values of the normalized CCF peaks for each order and the errors were given as the standard error of the weighted mean. Typical uncertainties for this method are |$\lesssim$|1 km s|$^{-1}$|⁠, regardless of |$v\sin i_{\star }$| or spectral type. When one or multiple orders give |$v\sin i_{\star } \le 2$| km s|$^{-1}$|⁠, we treat the final combined |$v\sin i_{\star }$| value as an upper limit.

2.2 Stellar radius

Stellar radii of young (single) stars are commonly estimated using the Stefan–Boltzmann relation (e.g. Mann et al. 2016; Davies 2019), the scale factor between a template or model spectrum and an absolutely calibrated spectrum of the star (the infrared-flux method; Blackwell & Shallis 1977; Casagrande et al. 2010; Newton et al. 2019), and/or interpolation from a grid of stellar models (e.g. Muirhead et al. 2012; Mayo et al. 2018; Loaiza-Tacuri et al. 2023). The first two methods tend to be applied when handling a single star or a set of similar systems (e.g. just field M dwarfs). The reason is because they depend on the choice of template, model, and available photometry, which varies significantly with age, effective temperature (⁠|$T_{\mathrm{eff}}$|⁠), surface gravity, and extinction. When working with a more diverse set of stars, it is preferable to use evolutionary models, where one can achieve results across the mass function at scale.

The challenge for evolutionary models is they tend to underpredict the radii of cold (e.g. Mann et al. 2015) and/or young stars (e.g. Rizzuto et al. 2016). However, models from the PAdova and TRieste Stellar Evolution Code (PARSEC; Bressan et al. 2012) contain an empirical correction for this. Models from the Dartmouth Stellar Evolution Program (DSEP; Chaboyer et al. 2001; Dotter et al. 2008) also reproduce observed properties of young stars (e.g. Kraus et al. 2015; David et al. 2019b) after adding in effects of magnetic fields (e.g. Feiden 2016).

For this work, we focus on estimating |$R_{\star }$| from a stellar evolution model given multiband photometry. To this end, we have developed stelpar, a Python-based pipeline and analysis tool for estimating |$R_{\star }$| and other stellar parameters from a general set of observational photometry (and a parallax) and an input model grid.² The code is similar to others (e.g. Koposov 2023), but ensures a homogeneous treatment across stellar types.

2.2.1 Stellar evolutionary models

We consider three evolutionary model grids: non-magnetic DSEP, a DSEP-based grid with magnetic enhancement (Feiden & Chaboyer 2012), and a PARSEC grid. Both DSEP-derived models cover ages across the entire pre-main- and main-sequence (1 Myr to 10 Gyr) and stellar masses (⁠|$M_{\star }$|⁠) from 0.09 to 2.45 |$\rm M_{\odot }$|⁠, depending on age. The PARSEC grid has ages ranging from 1–500 Myr.³

Prior work on young stars often used the MESA Isochrones & Stellar Tracks (MIST; Choi et al. 2016), SPOTS models (Somers, Cao & Pinsonneault 2020), or BHAC15 (Baraffe et al. 2015). The first does not include the effects of magnetic fields, spots, or other activity. As such, they tend to give ages significantly lower than other methods, like lithium depletion (Malo et al. 2014; Wood et al. 2023b). The SPOTS models only go to 1.3 |$\rm M_{\odot }$|⁠, while our sample would go to at least 1.5 |$\rm M_{\odot }$| (age-dependent cutoff). The BHAC15 models grid is more coarse and has a more limited set of photometry, which severely limits the precision of the output parameters.

Studies of star-forming regions near the Sun have found most associations are Solar or slightly sub-Solar metallicity (Spina et al. 2017). We only consider Solar metallicity for these models. Tests on slightly sub-Solar metallicities did not significant impact the final radii.

To lessen the computational cost of interpolating the model grid at every iteration, we used DFInterpolator from the isochrones Python package (Morton 2015) to pre-interpolate each DSEP model bilinearly in age and mass. This results in a grid spacing of |$\le$|10 per cent in age (e.g. 0.1 Myr from 1–10 Myr, 1 Myr from 10–100 Myr, etc.), and 0.005 |$\rm M_{\odot }$| for |$M_{\star } < 0.1\, \mathrm{M}_{\odot }$| and 0.01 |$\rm M_{\odot }$| for |$M_{\star } \ge 0.1\, \mathrm{M}_{\odot }$| in mass. We did not repeat this for the PARSEC model since we could choose the age spacing when we downloaded the grid. As a result, the mass spacings are not the same between all the models, but spacings were still smaller than the measurement uncertainties.

2.2.2 Input photometry

For all fits, we download (when available) photometry from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006); the Gaia mission (Prusti et al. 2016) Data Release 3 (DR3; Vallenari et al. 2023); the Sloan Digital Sky Survey (SDSS) Data Release 16 (DR16; Ahumada et al. 2020); B and V from the AAVSO Photometric All-Sky Survey (APASS; Henden et al. 2009) Data Release 9 (DR9; Henden et al. 2015); H|$_{\textnormal {p}}$| from the Hipparcos mission (ESA 1997); and B|$_{\textnormal {T}}$| and V|$_{\textnormal {T}}$| from the Hipparcos TYCHO-2 catalogue (Høg et al. 2000). We did not include the Pan-STARRS1 survey (PS1; Chambers & Pan-STARRS Team 2018) photometry for now, as our test samples are too bright (beyond PAN-STARRS saturation limit) and due to evidence of offsets for cooler stars (Kado-Fong et al. 2016). We also found that most fits are not limited by the photometry; Pan-STARRS and similar surveys (e.g. SkyMapper; Wan et al. 2018) can be included following the same method if needed. We also excluded photometry from NASA’s Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), which would likely be contaminated by emission from the disc.

2.2.3 Photometry model construction

We compared the model to observed photometry within a Markov chain Monte Carlo (MCMC) framework with the emcee Python package (Foreman-Mackey et al. 2013). At each iteration, we interpolate the evolutionary model grid by age and |$M_{\star }$|⁠, neglecting metallicity (with Solar metallicity the model point is uniquely determined by only age and |$M_{\star }$|⁠).

We used models with close grid spacings (either from pre-interpolation or custom download; see Section 2.2.1), so there is no advantage to bilinear interpolation (e.g. with DFInterpolator from isochrones; Morton 2015) which comes with a significant decrease in computation speed. Therefore, we choose to interpolate by searching for the nearest-neighbour value in age, and linearly interpolating in |$M_{\star }$|⁠. This interpolation method can create a bias towards grid points, since we are choosing discrete age values. However, this does not change any of our results because the age and |$M_{\star }$| uncertainties are much larger than the grid spacing.

In total, stelpar contains four free parameters: age, |$M_{\star }$|⁠, extinction (⁠|$A_V$|⁠), and f. The first two are set by the models. For |$A_V$|⁠, we use the synphot Python package (STScI Development Team 2018) following the extinction model of Cardelli, Clayton & Mathis (1989). The final parameter, f is a factor describing the underestimation of the errors on the measured photometry and/or uncertainties in the model, measured in magnitudes. In practice, f would be more accurately modelled as a vector in wavelength and age, which would account for the fact that models tend to struggle more at optical wavelengths due to missing opacities in cool stars (Mann, Gaidos & Ansdell 2013), stronger variability in the blue (Gully-Santiago et al. 2017; Mori et al. 2024), and accretion (Herczeg & Hillenbrand 2008). Until such external constraints are available, we kept this as a single number applied to all photometry for a given star.

When comparing the model to photometry, stelpar assumes Gaussian uncertainties. By default, the code includes uniform priors within a set of user-configurable bounds for each fit parameter. However, the user may optionally choose to define Gaussian priors for any of the fit parameters and/or |$T_{\mathrm{eff}}$|⁠, the latter of which acts as a prior on both age and |$M_{\star }$|⁠. Since the parameters above uniquely determine the model selection, the full posterior distributions for age and |$M_{\star }$| can be turned into posteriors on |$R_{\star }$|⁠, |$T_{\mathrm{eff}}$|⁠, log-surface gravity (⁠|$\log g$|⁠), log-luminosity (⁠|$\log L$|⁠), and stellar density (⁠|$\rho _{\star }$|⁠).

2.2.4 Comparison to empirical values

To test our method and check for additional systematic uncertainties, we compared the inferred |$\rho _{\star }$| to those for seven young transiting systems. In cases of low eccentricity, the transit duration can yield a strong constraint on the stellar density, |$\rho _{\star }$| (Seager & Mallén-Ornelas 2003). The transiting sample encompasses a reasonable range of stellar masses (M dwarfs to F stars) and includes stars with transitional and debris discs (3–20 Myr).

A comparison in |$\rho _{\star }$| can be misleading because inaccurate |$R_{\star }$| can be masked by matching inaccurate |$M_{\star }$| values. However, since |$\rho _{\star }$| is scaled steeply with |$R_{\star }$| (cubic), even (modest) 10–15 per cent precision on |$\rho _{\star }$| would only lead to 3–5 per cent errors in |$R_{\star }$|⁠.

The systems were selected based on their age and planet multiplicity. The youngest systems are expected to have low eccentricity (due to dampening from the protoplanetary disc; Papaloizou & Larwood 2000). Multiplanet systems tend to have low eccentricities and yield stronger constraints on stellar density (Van Eylen & Albrecht 2015).

We show the comparison in Fig. 3 with supplementary information in Table 1. Most of our results agreed with observations regardless of the model. The DSEP-magnetic model result for K2-136 shows a significant offset, but this is as expected. K2-136 is the oldest target in the sample (⁠|$\simeq$|800 Myr) and is relatively quiet; magnetic enhancement is probably not required. Indeed, the non-magnetic DSEP model gave better agreement. We did not use the PARSEC model for K2-136 because our custom PARSEC grid only goes to 500 Myr in age.

Results for HIP 67522 also show disagreement. This is unlikely to be due to our assumption that the planetary orbits are circular; a second planet has since been identified in that system (Barber et al. 2024b) and the lack of dynamical interactions between the planets (transit timing variations) suggest the masses and eccentricities of both planets are low (Thao et al. 2024). However, a large part of this disagreement is because the models yield exceptionally small uncertainties, sometimes better than 2 per cent on radius. If we adopt a flat 5 per cent uncertainty on radius (Tayar et al. 2022), there is agreement across the sample.

One surprising piece is that the DSEP non-magnetic model is consistent too. On further inspection, this is because we had to use |$\rho _{\star }$| rather than a direct comparison in |$M_{\star }$|⁠. Non-magnetic DSEP is yielding radii that are, on average, smaller than the values from the other two grids and smaller when compared to radii measurements from direct spectral analysis in the discovery papers. However, the masses are also smaller, which masks some of the difference. Another effect is that we placed priors on the age, derived from the parent population. Upon removing this prior, the DSEP non-magnetic model tended to yield discrepant ages, as seen in more extensive studies (Rizzuto et al. 2016).

Overall, we find that our model-based radii are reliable, although the uncertainties are underestimated. This holds down to the 3 Myr system with a protoplanetary disc (IRAS 04125+2902; Barber et al. 2024a). In absolute terms, the agreement suggests we can achieve a |$\simeq$|5 per cent uncertainty on |$R_{\star }$|⁠. There is also no evidence of a systematic offset. The DSEP-magnetic model appears to perform the best, although PARSEC does similarly well. While the sample is small, they are broadly consistent with prior comparisons in mass–radius space using young eclipsing binary (EB) systems (e.g. Kraus et al. 2015; Gillen et al. 2017; David et al. 2019b). We note that such studies often find issues with the |$T_{\mathrm{eff}}$|–luminosity scale, but our method appears less impacted by this offset.

2.3 Stellar rotation period

Young stars (⁠|$\lesssim$|10 Myr) tend to have |$P_{\textnormal {rot}}$| values from 0.1 to 30 d. The long-period end of this distribution is challenging for the Kepler-K2 Mission and the Transiting Exoplanet Survey Satellite (TESS) due to the narrow observing windows (27–80 d). However, the peak in the period distribution is |$\simeq$|2 d, and only a small percentage of stars have |$P_{\textnormal {rot}}$| values longer than 10 d (Rebull et al. 2018). Therefore, we used a search grid from 0.2 to 20 d, which captures the vast majority of rotators. We adopted the |$P_{\textnormal {rot}}$| value that corresponds to the highest peak in the LS power spectrum, with an eye inspection to check for multiple periods and aliases.

For each star we generated TESS light curves by using a causal pixel model (CPM) as implemented in the unpopular package (Hattori et al. 2022) and detailed in Barber et al. (2022). unpopular creates light curves using the pixels outside the target aperture to model the pixel response of the pixels within the aperture. We subtracted the systematic model from the raw aperture light curve, which results in the CPM light curve.

Boyle et al. (in preparation) estimated uncertainties on |$P_{\textnormal {rot}}$| values measured in this way. Their method was to compare periods estimated from different TESS sectors as well as between K2 and TESS data of the same set of stars. In principle, this method includes effects like spot evolution and differential rotation, provided the gap between sectors and between K2 and TESS is long compared to the evolutionary time-scale. They found that |$P_{\textnormal {rot}}$| can be determined to better than 5 per cent for stars with rotation periods below 10 d, and |$\simeq$|2 per cent for typical rotators in the disc-bearing sample (⁠|$P_{\mathrm{rot}} \simeq 2$| d). Thus, we expect |$P_{\textnormal {rot}}$| will have little impact on the total error budget of the final result.

2.4 Stellar inclination

Deriving |$i_{\star }$| from the |$v\sin i_{\star }$|⁠, |$R_{\star }$|⁠, and |$P_{\textnormal {rot}}$| values we measured/inferred requires more attention beyond a direct application of equation (1). It is possible, for instance, to measure |$v\sin i_{\star } > v_{\mathrm{eq}}$| (which is non-physical) simply due to random measurement uncertainties (see e.g. Morton & Winn 2014).

Masuda & Winn (2020) presented a framework for measuring |$\cos i_{\star }$| (rather than |$i_{\star }$|⁠) that accounts for the statistical co-dependence of |$v_{\textnormal {eq}}$| and |$v\sin i_{\star }$|⁠. Briefly, the authors state that it is incorrect to sample |$v_{\textnormal {eq}}$| and |$v\sin i_{\star }$| independently. Instead, it is more appropriate to sample |$v_{\textnormal {eq}}$| and |$\cos i_{\star }$| independently and derive |$v\sin i_{\star }$| from the combination of the two. Following Morton & Winn (2014) and Masuda & Winn (2020), we estimate |$\cos i_{\star }$| within an MCMC framework using emcee.⁴

Uncertainties on |$\cos i_{\star }$| can be large, corresponding to inclination uncertainties of 10°–30|$^\circ$|⁠, with significant variation with both |$v_{\textnormal {eq}}$| and |$i_{\star }$|⁠. In this sense it is challenging to determine |$i_{\star }$| of any given target beyond broad categories. However, we can still combine many such measurements to study the overall alignment frequency.

2.5 Disc inclination

|$i_{\textnormal {d}}$| derived from ALMA observations are widely available in the literature. To determine typical uncertainties (⁠|$\sigma _{i_{\mathrm{d}}}$|⁠) from such measurements, we drew from multiple prior surveys which used ALMA. Specifically, we drew |$i_{\textnormal {d}}$| measurements from Ansdell et al. (2016, 2020), Barenfeld et al. (2017), Tazzari et al. (2017), Tripathi et al. (2017), Huang et al. (2018), and Long et al. (2019). Collectively, this provided 138 measurements from 107 discs spanning the protoplanetary disc lifetime (⁠|$\sim$|1–10 Myr). The systems were surveyed from four regions: the |$\rho$| Ophiuchus star-forming region (⁠|$\sim$|1 Myr; Andrews & Williams 2007), the Lupus star-forming complex (⁠|$\sim$|3 Myr; Alcalá et al. 2017), the Taurus Molecular Cloud (⁠|$\lesssim$|6 Myr; Krolikowski, Kraus & Rizzuto 2021), and the Upper Scorpius OB association (⁠|$\sim$|10 Myr; David et al. 2019b).

We show the results of our analysis in Fig. 4. For 68 per cent and 95 per cent of measurements, uncertainties were |$\lesssim 4^{\circ }$| and |$\lesssim 21^{\circ }$|⁠, respectively. There was no significant trend in the uncertainties with disc orientation or age, although the overall detection rate drops with increasing age. There is an almost uniform spread of measurements between face-on and edge-on discs, as expected for (nearly) complete surveys. The exception is a slight deficit of systems with |$80^{\circ } < i_{\mathrm{d}} < 90^{\circ }$|⁠; this is likely due to the fact that edge-on discs will occult their host and therefore be less likely to be included in the ALMA observing list.

As a check on the literature |$i_{\textnormal {d}}$| measurements, we compared 26 systems that had independent |$i_{\textnormal {d}}$| measurements for the same system in at least two or more surveys. These are not necessarily independent measurements; many are using the same underlying ALMA data. Rather, this is a test of sensitivity to the fitting method. We show the results of our literature comparison in Fig. 5. Of the 26 overlapping systems, 24 had measurements that agreed across multiple surveys within 1|$\sigma$|⁠, and the other systems’ measurements agreed within 2|$\sigma$|⁠.

The distribution seen in Fig. 4 shows a deficit of edge-on discs and a surplus of intermediate-angle (⁠|$i_{\mathrm{d}}\simeq 45^\circ$|⁠) discs compared to the expected uniform |$\cos {i_{\mathrm{d}}}$| distribution. The former is due to observational bias. Discs flare outwards by 10°–25|$^\circ$|⁠, so those with |$i_{\mathrm{d}}\gtrsim 70^\circ$| are more likely to be occulted/extincted by the disc. The lack of a visible optical component means these were more likely to be skipped in the sub-mm or Spitzer survey that initially identified the infrared excesses. This bias against edge-on discs can be included in the model (Section 2.6).

The excess at moderate inclinations appears to be a combination of disc measurements with high uncertainties and random noise. Given Poisson errors, the observed and expected distributions are marginally consistent with each other out to |$i_{\mathrm{d}}=75^\circ$|⁠, particularly if we only consider discs with |$\sigma _{i_{\mathrm{d}}}< 10^\circ$|⁠.

Overall, this indicates |$i_{\textnormal {d}}$| is a small component of the total error budget (particularly compared to |$i_{\star }$|⁠), and that using literature determinations from ALMA data and assuming the reported Gaussian uncertainties will not introduce a bias in the overall disc–star alignment measurement. Although the bias against edge-on systems likely needs to be included in the model.

2.6 Fitting for global alignment

As noted above, uncertainties on |$i_{\star }$| can be large in any individual system. However, we can still use these to draw conclusions about the overall alignment distribution within a population, provided we can control for systematic effects. To this end, we combine |$i_{\star }$| and |$i_{\textnormal {d}}$| for all systems within a hierarchical Bayesian model (HBM), where each disc–star offset (⁠|$\alpha ^{\prime }_n$|⁠) is a free parameter evolving under the observational constraints (e.g. |$v\sin i_{\star }$|⁠, |$R_{\star }$|⁠, |$i_{\textnormal {d}}$|⁠) simultaneously with global parameters that describe the population-level disc–star alignment.

The method is flexible to the assumed model of disc–star alignment. For simplicity, we assume the overall (global) alignment can be modelled as a Gaussian described by two parameters (⁠|$\mu$| and |$\sigma$|⁠). Similarly, we also assume the alignments of the individual systems can be modelled as Gaussians (see Section 4 for more detail). Combined with the individual systems the assumed Gaussian global alignment yields |$N+2$| fit parameters. Physically, this corresponds to an offset in the alignment (⁠|$\mu$|⁠) and a spread around that offset (⁠|$\sigma$|⁠). It is just as simple to describe it with a Fisher distribution following Morton & Winn (2014), a uniform distribution (e.g. random alignment), or a combination of the two with a mixture amplitude (e.g. a mix of randomly aligned and aligned systems).

We define the log-probability (⁠|$\log \mathcal {P}$|⁠) in two parts that are summed together: the sums of the log-likelihoods (⁠|$\log \mathcal {L}$|⁠) of observing the measured ensemble of alignments given |$\alpha ^{\prime }_n$| and of observing |$\alpha ^{\prime }_n$| given |$\mu$| and |$\sigma$|⁠:⁵

where |$i_{\star , n}$| and |$i_{d, n}$| represent the measured inclinations of each system and N is the total number of systems.

In the first part of equation (2), we compare the estimated |$i_{\star }$| (converted from |$\cos i_{\star }$|⁠; Section 2.4) and |$i_{\textnormal {d}}$| (Section 2.5) with |$\alpha ^{\prime }_n$| and scale them by the combined uncertainties on |$i_{\star }$| and |$i_{\textnormal {d}}$|⁠, i.e.

In the second part, we compare |$\mu$| to |$\alpha ^{\prime }_n$|⁠, scaled by |$\sigma$|⁠, i.e.

In this way, |$\mu$| and |$\sigma$| are constrained by |$\alpha ^{\prime }_n$|⁠, which themselves are modulated by the previously derived |$i_{\star }$| and |$i_{\textnormal {d}}$|⁠.

2.6.1 Additional effects

Our treatment above does not consider the full 3D alignment, only the line-of-sight-projected alignment. This issue is discussed in prior studies on both disc–star alignment (Davies 2019) and more extensively on planet–star alignment (e.g. Dong & Foreman-Mackey 2023). Since there is a missing angle, a misaligned system may appear aligned because most/all of the misalignment is hidden in the unseen angle. Similarly, systems mostly aligned in three dimensions might look preferentially misaligned if all of the misalignment is between the observed |$i_{\star }$| and |$i_{\textnormal {d}}$|⁠.

While this limits our ability to measure overall alignment in any given system, we could infer that a 2D treatment would only affect global |$\sigma$| in an ensemble analysis. Assuming there is no preferred direction to the misalignment or the observer angle, then the amount that any misalignment is hidden in the unseen angle should be effectively random. Thus, when |$i_{\star }$| and |$i_{\textnormal {d}}$| are aligned, any misalignment must come from the difference in sky-projected stellar and disc angles (⁠|$\lambda _{\star }$| and |$\lambda _{\rm d}\lambda$|⁠, respectively). Indeed, modelling |${\rm d}\lambda _{\star }$| and |$\lambda _d$| (or their difference) does not meaningfully change our HBM results except that the true misalignment uncertainty is larger. The impact is also smaller than many of the other effects discussed in this paper (e.g. inclination uncertainties).

This breaks down if there is a preferred direction of misalignment, such as if disc warping were driven by the overall angular momentum of the parent star-forming cloud. In that case, the degree of misalignment in the observable direction is correlated between systems. The solution here is to observe over many star-forming regions and compare, since any preferred direction would be different between the populations.

As we show in Fig. 6, measurements of both the disc and star projected onto the sky may lead us to thinking misaligned systems are aligned. For the star, |$v\sin i_{\star }$| only measures the tilt of the star towards or away from the observer, not the direction. For the disc; there is no way to tell which part of an inclined disc is in front of versus behind a star, so identical protoplanetary discs with line-of-sight-flipped orientations will have the same measured |$i_{\textnormal {d}}$| values.

For an ensemble analysis, the orientation mirroring has a minor impact on our ability to measure the alignment distribution, other than requiring a larger sample. A more subtle but related effect is the barrier at inclinations of |$0^\circ$| and |$90^\circ$| (both star and disc). This could make systems look more aligned than they are. This effect is included in our HBM example, and the major impact is to systematically overestimate |$\mu$| and underestimate |$\sigma$| (see Fig. 10). The size of this bias is small compared to other effects discussed in this paper, but larger than the uncertainties for sample sizes more than |$\simeq$|30. The solution for this problem is either to model a simulated population and apply a correction, or include the effect as a parameter of the simulation.

2.7 Multistar systems

Multistar systems can complicate estimating |$R_{\star }$| (unresolved systems will appear brighter), |$v\sin i_{\star }$| (unresolved but separated lines will increase apparent broadening), and |$P_{\textnormal {rot}}$| (there may be two rotation signatures in the light curve). Modelling these effects can be difficult because the impact depends on the separation and contrast of the companion. Fortunately, these generally only impact binaries where the components are unresolved and have low mass ratios where there is still significant flux from the companion. These represent a small fraction of star systems. Double-line spectroscopic binaries (SB2s) and high-order multiples, for example, make up only |$\sim$|3 per cent of main-sequence stars (Kounkel et al. 2021). Equal-mass binaries, while over-represented in the binary population (El-Badry et al. 2019), are also the easiest to identify.

Separately, the disc–star alignment distribution is likely different for close binaries as it is for single stars and wide binaries. Discs in close and intermediate binaries are expected to be shorter lived, and their orientations should be heavily influenced by the companion through disc truncation (e.g. Artymowicz & Lubow 1994; Andrews et al. 2010; Williams & Cieza 2011; Jang-Condell 2015), increased accretion (e.g. Artymowicz & Lubow 1994; Jensen et al. 2007), and increased photoevaporation (e.g. Alexander 2012; Rosotti & Clarke 2018). Discs orbiting both binary components (circumbinary discs) are also believed to be rare and may only live for a few Myr (e.g. Akeson et al. 2019; Czekala et al. 2019; Offner et al. 2023).

Fortunately, if the disc is not present it would automatically not be included in a disc–star alignment sample; the input is necessarily the list of targets with a resolved disc. Thus, most of the binaries that would impact our stellar parameters will be removed because the disc has dissipated below detectable levels prior to any observations.

The best solution is to remove all close-in binaries from the sample, although complete removal is challenging. Gaia Renormalized Unit Weight Error (RUWE) is an indicator of unresolved companions (e.g. Stassun & Torres 2021). RUWE is most sensitive to a specific range of binary orbits (Wood, Mann & Kraus 2021), but this includes most near-equal binaries that are problematic for our work. Young stars have higher RUWE than their older counterparts (Fitton, Tofflemire & Kraus 2022), so we adopt a more generous cutoff of |$\mathrm{RUWE}< 10$|⁠. Separately, we remove any system that shows two sets of lines in the spectrum or two peaks in the CCFs from the |$v\sin i_{\star }$| analysis (Section 2.1). We remove any target where stelpar yields an age inconsistent with the group age; this is typically a sign of an elevated colour–magnitude diagram (CMD) position due to an unresolved binary. For stars with ALMA data, we can also remove any target with a clearly resolved tight companion.

To explore the effectiveness of these cuts, we adopt the RUWE cut above and ALMA sensitivity and beam size typically used for disc morphology (Ansdell et al. 2017). We generated a population of binaries using the MOLUSC code (Wood et al. 2021), which generates a realistic sample of stellar companions and determines which survive a set of input data (including RUWE). The problematic binaries that may hamper measurements but not deplete the disc or be detected in the above data are generally those with orbits |$<$|0.1 au, which are most likely to be detected as spectroscopic binaries. These also represent a small fraction of the overall binary population (⁠|$<$|1 per cent). We conclude these cuts are sufficient to make a clean sample.

3 APPLICATION TO PRE-MAIN-SEQUENCE STARS WITHOUT DISCS

As a test, we calculated stellar parameters (following Sections 2.1–2.3) for a sample of members mostly in the |$\beta$| Pictoris Moving Group (⁠|$\beta$|PMG), as well as the Tucana-Horologium Moving Group (Tuc-Hor) and the Carina-Extended Association (Car-Ext; Luhman 2024), and compared them to literature estimates. Members of |$\beta$|PMG, Tuc-Hor, and Car-Ext are ideal to test our methods for stellar parameters because they are PMS with ages |$\sim$|11–26 Myr (Couture, Gagné & Doyon 2023), |$\sim$|40 Myr (Kraus et al. 2014), and |$\sim$|34–44 Myr (Wood et al. 2023b; Luhman 2024), respectively. These stars are nearby, and many have data from IGRINS. Additionally, stars in this association should not have protoplanetary discs, eliminating the bias against high |$i_{\textnormal {d}}$| values. The result is that |$\cos i_{\star }$| should be uniformly distributed.

We chose our sample based on available IGRINS data. We downloaded archival IGRINS K-band spectra for 25 stars with spectral classes ranging from early G- to late M-type. We provide general information for the sample in Table 2. We removed 5 of 25 stars that were SB2s and/or had very high RUWE values (⁠|$>$|20).

We calculated |$v\sin i_{\star }$| and |$R_{\star }$| values for the 20 remaining stars following Section 2.1 and Section 2.2, respectively. We generated light curves using TESS data downloaded from the Barbara A. Mikulski Archive for Space Telescopes (MAST) for all but one star in our sample (HD 358623), which did not have TESS data, and calculated |$P_{\textnormal {rot}}$| values and errors (following Section 2.3). We provide our estimated stellar parameters alongside literature values in Table 3.

We compare our results to those from the literature in Fig. 7. Specifically, we show the calculated equatorial velocities (⁠|${v_{\mathrm{eq}} = 2\pi R_{\star } / P_\mathrm{rot}}$|⁠), using |$R_{\star }$| and |$P_{\textnormal {rot}}$| from the literature and from our own estimates, and compared them to literature |$v\sin i_{\star }$| and our estimated |$v\sin i_{\star }$|⁠, respectively, via fractional residual |$(v \sin i_{\star } - v_{\mathrm{eq}})/v_{\mathrm{eq}}$|⁠.

Since there are more edge-on inclinations than face-on ones, we expect there to be more stars near zero (equality) than at |$\simeq -1$| and some stars may be over one due to random errors. However, the literature results are extreme even considering this.

As a test, we generated a fake sample of systems assuming random |$\cos i_{\star }$| and properties (including measurement uncertainties) matching the sample in Fig. 7. For our measurements, the K-S test gives a p-value of 90 per cent, suggesting values are consistent with the expected distribution. For the literature, the result depends on how we select the sample. Using just the best measurements yield a p-value of just 0.3 per cent (inconsistent). Combining all measurements gives 9 per cent although this double counts the same stars (not done in the synthetic sample; see Section 4). Averaging multiple measurements for the same star (many of which do not agree with each other) dropped this to 0.2 per cent.

This effect has been noted in prior studies (e.g. Newton et al. 2016). Literature |$v\sin i_{\star }$| values tend to be overestimated, likely due to incomplete consideration of effects discussed in Section 2.1. This would lead to an overestimate in the number of edge-on stars. Our measurements do not seem impacted by this, which offers an observational confirmation of the methods outlined here.

4 APPLICATION TO A SYNTHETIC SAMPLE

As a test of the issues discussed in Sections 2.6 and 2.6.1, we constructed a synthetic sample of systems with realistic stellar and disc parameters and applied the methods of Sections 2.4 and 2.6. First, we created a parent alignment distribution, which we assumed to be a normal distribution centred at a true mean (⁠|$\mu _\mathrm{true}$|⁠) and scaled by a true intrinsic scatter (⁠|$\sigma _\mathrm{true}$|⁠), both of which we manually selected. These are defined such that a sample with perfect disc–star alignment would have |$\mu _\mathrm{true}=\sigma _\mathrm{true}=0$|⁠.

Next, we drew a sample of N systems from the parent population. We defined a |$\cos i_{\star }$| distribution as a uniform distribution of N points in the interval [-1, 1], and we calculated |$i_{\textnormal {d}}$| values from the difference between the generated |$i_{\star }$| values and the Gaussian alignment distribution. Since |$i_{\star }$| and |$i_{\textnormal {d}}$| are measured from 0–90|$^{\circ }$| (0–1 in cosine), as we discussed in Section 2.6.1, we forced |$i_{\star }$| and |$i_{\textnormal {d}}$| to 0–90|$^{\circ }$| by taking the absolute value of the cosine of the angle.

We defined |$R_{\star }$| and |$P_{\textnormal {rot}}$| distributions as log-uniform distributions of N points in the intervals [0.1, 1.5] |$R_{\odot }$| and [0.2, 12] d, respectively. While not a perfect match to an observed star-forming population, the results were not sensitive to changes in the |$R_{\star }$| and |$P_{\textnormal {rot}}$| distributions. We then calculated |$v_{\textnormal {eq}}$| values from the assigned |$R_{\star }$| and |$P_{\textnormal {rot}}$| values. By chance, some stars had unphysically large |$v_{\textnormal {eq}}$| (above breakup speeds), so these were adjusted manually. Last, we calculated |$v\sin i_{\star }$| from a combination of |$\cos i_{\star }$|⁠, |$R_{\star }$|⁠, and |$P_{\textnormal {rot}}$|⁠.

At this point all assigned parameters have no uncertainties (i.e. these are true values instead of measured ones). So, we assigned each measured parameter an uncertainty based on the empirical tests in prior sections. Following Section 2.1, we used 1 km s|$^{-1}$| uncertainties for all |$v\sin i_{\star }$| values. For |$R_{\star }$| we used 5 per cent uncertainties. For |$P_{\textnormal {rot}}$| we calculated uncertainties following Boyle et al. (in preparation). We gave all |$i_{\textnormal {d}}$| values an uncertainty of 4|$^{\circ }$|⁠, following our analysis in Section 2.5.

4.1 Test: stellar inclination

One concern from our approach is that our likelihood assumes the |$\cos i_{\star }$| distributions can be approximated as Gaussians. We can test this with the synthetic sample. If we do not remove edge-on discs, the resulting stellar inclinations should be uniformly distributed in |$\cos i_{\star }$|⁠, while a non-uniform distribution would suggest problems (especially near the extreme inclinations) with the Gaussian assumption.

We drew 25 synthetic systems from the above population, from which we derived |$\cos i_{\star }$| values. We chose 25 systems for this test arbitrarily, and we anticipated similar results for a larger number of systems, which we discuss further in Section 4.3.

We compared our estimated values to a uniform distribution using empirical cumulative distribution functions (ECDFs). We took 1000 random draws from each |$\cos i_{\star }$| posterior and calculated the ECDF at each draw (i.e. we had 1000 ECDFs, each derived from 25 random |$\cos i_{\star }$| values). Similarly, we took 1000 random draws of 25 data points from a uniform distribution and calculated the ECDFs at each draw. We calculated the mean and standard deviation of each ensemble of ECDFs, which we show in Fig. 8.

For comparison, we also calculated the ECDF of the originally defined |$\cos i_{\star }$| distribution (before any perturbation), which was essentially a single random draw of 25 data points from a uniform distribution. Fig. 8 shows that the ECDFs calculated from the |$\cos i_{\star }$| posteriors and the uniform distribution are in agreement with one another, which is what we expect.

4.2 Test: alignment distribution

As a test on our HBM (Section 2.6), we used the same setup as in the previous test and inferred |$\mu$| and |$\sigma$| for the parent alignment distribution, and the alignment values (⁠|$\alpha _n^{\prime }$|⁠) for the individual systems. In this case, we did not account for the fact that edge-on discs are disfavoured observationally.

In Fig. 9, we show our result for the inferred |$\mu$| and |$\sigma$| as a 2D contour, accompanied by 1D probability densities for each |$\alpha _n^{\prime }$|⁠. The estimated |$\mu$| and |$\sigma$| agreed with the initial alignment distribution to within 1|$\sigma$|⁠.

We show the HBM fit to the individual systems in Fig. 9, where they were all in agreement with the alignments derived through the combination of |$i_{\star }$| (following Section 2.4) and |$i_{\textnormal {d}}$| (following Section 2.5). Note that the method assumes Gaussian distributions, which is not accurate for some systems. As we discuss below, this does not prevent recovering the global |$\mu$| and |$\sigma$|⁠, but likely leads to underestimated uncertainties on these parameters.

4.3 Dependence on the number of systems

To explore what kind of sample size is required to retrieve the disc–star alignment distribution, we varied the number of systems from 5 to 150 and re-ran our HBM to estimate |$\mu$| and |$\sigma$|⁠. For this, we used the same setup as in previous tests and kept the input |$\mu$| and |$\sigma$| constant (except we increased |$\sigma$| from previous tests), although the result shows only weak dependence on the exact |$\mu$| and |$\sigma$| used. The resulting underestimation is more easily seen at larger N where the uncertainties on |$\mu$| and |$\sigma$| are smaller.

We show these results in Fig. 10. At lower N, the variation around the input values is expected given the uncertainties. However, past |$N\simeq 50$|⁠, the recovered |$\mu$| and |$\sigma$| vary around the input values by more than the expected uncertainties. This is driven by the assumption that the individual |$\alpha ^{\prime }_n$| estimates are Gaussian, while many are asymmetric (Fig. 9).

The bias caused by enforcing inclinations to 0–90|$^{\circ }$| (see Section 2.6.1) is most apparent in |$\sigma$|⁠, where points are statistically |$\simeq 3^\circ$| below the input values. On a real data set, this can be corrected either by including it in the model, or generating a synthetic sample like this one and applying a correction. In either case, the bias depends on the underlying model. A tight distribution of aligned systems, for example, shows an almost negligible bias. This suggests exploring a few different possible distributions in the final fit as a test of sensitivity to such assumptions.

5 SUMMARY AND CONCLUSIONS

In this paper, we considered the important factors that contribute to measuring the distribution of disc–star alignment angles for a large number of stars hosting protoplanetary discs. We explored what uncertainties are realistic for the input parameters, including |$R_{\star }$|⁠, |$v\sin i_{\star }$| |$P_{\textnormal {rot}}$|⁠, and |$i_{\textnormal {d}}$| using existing observations of infant stars. We applied our methodology to both real and synthetic data sets to evaluate the impact of assumptions in our methods.

We summarize the main results as follows:

• |$v\sin i_{\star }$|⁠: Following Kesseli et al. (2018) we were able to estimate |$v\sin i_{\star }$| with uncertainties as good as |$\simeq$|1 km s|$^{-1}$| for G- through M-type PMS stars, including disc-bearing stars.

• |$R_{\star }$| and stelpar : With stelpar, we are able to estimate |$R_{\star }$| for PMS stars to |$\simeq$|5 per cent. Given the complexities of estimating the fundamental properties of young stars, this is surprisingly precise. However, this method reproduces empirical densities from transits and eclipsing binaries, including for at least one disc-bearing star (IRAS 04125+2902).

• Archival |$i_{\textnormal {d}}$| Measurements: ALMA |$i_{\textnormal {d}}$| measurements already exist in the literature (e.g. Huang et al. 2018) spanning the protoplanetary disc lifetime (e.g. |$\rho$| Ophiuchus, Lupus, Taurus, Upper Scorpius). Of these, 68 per cent had uncertainties |$\lesssim 4^{\circ }$|⁠. Overlapping measurements between surveys agree within reported uncertainties. While there is an observational bias against edge-on discs, the |$i_{\textnormal {d}}$| values follow the expected distribution for |$i_{\mathrm{d}}< 75^\circ$|⁠. We conclude that most previous |$i_{\textnormal {d}}$| measurements are precise and consistent between studies, and uncertainties/systematics are small compared to measurements related to |$i_{\star }$|⁠.

• PMS star literature comparison: We calculated stellar parameters for 20 young (⁠|$\sim$|11–44 Myr) PMS stars (without discs) within |$\beta$|PMG, Tuc-Hor, and Car-Ext. We found that literature |$v\sin i_{\star }$| estimates are systematically overestimated, yielding a non-random stellar inclination distribution. Our estimates reproduce the expected distribution.

• Global alignment distribution: Based on a synthetic sample, we find that an HBM analysis is sufficient to reproduce the input values within uncertainties, despite the necessary assumptions made along the way. Modest samples (⁠|$N\simeq 20$|⁠) are sufficient to identify any significant population of misaligned systems. At |$N \gtrsim 70$|⁠, systematics dominate over random samples, requiring more sophisticated modelling (e.g. consideration of non-Gaussian uncertainties).

We aim to use this method to estimate the disc–star alignment distributions for nearby young populations with protoplanetary discs. With increasingly available data from ALMA, TESS, K2, and Gaia, we need only IGRINS or similar NIR spectra of the sample. Many such samples already exist in the literature (e.g. López-Valdivia et al. 2021), suggesting this study may be possible with largely existing data.

Right now results on the alignment between discs and stars are limited primarily by the quality of the stellar parameters, the small sample of targets with high-quality data (high-resolution spectra with high signal-to-noise ratios and precise disc inclinations, e.g. from ALMA), and the methods used to turn these into population statistics. Our work has focused on the larger sources of bias and uncertainty, hence we can make significant improvements over prior studies (e.g. Davies 2019; Hurt & MacGregor 2023) while still making a number of simplifying assumptions in our HBM (e.g. ignoring the missing angle and assuming Gaussian uncertainties on |$i_{\star }$|⁠). As more data becomes available and astronomers develop better methods to measure fundamental stellar parameters, it will be more important to properly account for the non-Gaussian posteriors on |$\alpha$| and use a full 3D misalignment (i.e. marginalize over the missing angle).

DATA AVAILABILITY

The IGRINS spectra are available in the Raw & Reduced IGRINS Spectral Archive at https://igrinscontact.github.io/.

Photometry from 2MASS, Gaia, SDSS, APASS, TYCHO, and Hipparcos are available in online data bases, e.g. in the VizieR archive at https://vizier.cds.unistra.fr/viz-bin/VizieR-2.

TESS photometry is available in MAST at https://archive.stsci.edu/.

|$i_{\textnormal {d}}$| measurements were taken from literature references (cited above), but the underlying data is available in the ALMA archive at https://almascience.nrao.edu/aq/ or https://almascience.eso.org/aq/.

All new data derived in this work are available in the tables herein.

ACKNOWLEDGEMENTS

The authors thank the two referees for their comments; their feedback helped improve the paper. The authors thank Bernie and Bandit for their useful contributions. This work used The Immersion Grating Infrared Spectrometer (IGRINS). MJF acknowledges support from the North Carolina Space Grant Graduate Research Fellowship. This research was completed with funding from the National Aeronautics and Space Administration (NASA) Exoplanet Research Program (XRP) award #80NSSC21K0393.

Footnotes

REFERENCES

APPENDIX A: stelpar COMPUTATIONAL COST CONSIDERATIONS

Following Section 2.2.3, the pure-synphot extinction calculation can be computationally expensive. To work around this, we implemented a method which uses synphot to create the blackbody spectrum and the extinction model, but all of the relevant calculations are performed ‘by hand’ with numpy (Harris et al. 2020). The pared-down numpy approach already makes a noticeable speed improvement, but the extinction calculation could be performed tens of times for every iteration of the simulation, so it is important that our methods are as computationally inexpensive as possible.

To this end, stelpar employs numba just-in-time (JIT) compilation (Lam, Pitrou & Seibert 2015) to as many aspects of the extinction procedure as possible. JIT is a method which converts Python code to optimized machine code when first compiled. The first function call can oftentimes be more computationally expensive than the original (non-JIT) function, but every subsequent call will show vast speed improvements. Thus, it is a much more advantageous method to use when a function is called multiple times, as is the case with the extinction calculation. From the pure-synphot method to the combined numpy and numba method, we estimate a 10–15x increase in compilation speed depending on the number of bandpass filters used in the calculation.⁶

APPENDIX B: stelpar EXTINCTION CALCULATION

Following Section 2.2.3, the ‘by-hand’ numpy-based extinction calculation procedure is designed to reproduce the values derived by the pure-synphot method. The extinction value in a particular bandpass filter (⁠|$A_{\lambda }$|⁠) is calculated by

where |$\mathcal {S}_{\mathrm{eff}}^{\prime }$| and |$\mathcal {S}_{\mathrm{eff}}$| are the extincted and pure blackbody spectra, respectively. Effective stimulus is the flux density an observer would measure given a certain amount of flux generated from the synthetic (blackbody or extincted) source spectrum through a particular bandpass filter (in |$F_{\lambda }$| units; erg s|$^{-1}$| cm|$^{-2}$| Å|$^{-1}$|⁠). Effective stimulus is calculated by

with flux |$F(\lambda )$| in |$F_{\lambda }$| units, wavelength |$\lambda$| in Å, and the bandpass filter response function |$e(\lambda )$| (also called ‘bandpass transmission function’) in dimensionless fractions between 0 (no transmission) and 1 (full transmission).

Equations (B1) and (B2) give their results in |$F_{\lambda }$| units. If |$A_{\lambda }$| needs to be presented in different units (e.g. magnitudes, as is the case by default within stelpar), both |$\mathcal {S}_{\rm eff}^{\prime }$| and |$\mathcal {S}_{\rm eff}$| need to be converted individually before being used in equation (B1). When converting to AB magnitudes, for example, |$F_{\lambda }$| must be converted to |$F_{\nu }$| via

where c is the speed of light and |$\lambda _{\textnormal {piv}}$| is the pivot wavelength for a particular bandpass filter. CCDs follow the ‘equal-energy convention’ for |$\lambda _{\textnormal {piv}}$| given as

whereas the synphot method follows the ‘quantum-efficiency convention’ given as

Both pivot wavelengths lead to values of extinction that differ to |$\ll$|1 per cent in magnitudes.