The Comprehensive Archive of Substellar and Planetary Accretion Rates

Accretion rates ( Ṁ ) of young stars show a strong correlation with object mass (M); however, extension of the Ṁ–M relation into the substellar regime is less certain. Here, we present the Comprehensive Archive of Substellar and Planetary Accretion Rates (CASPAR), the largest compilation to date of substellar accretion diagnostics. CASPAR includes: 658 stars, 130 brown dwarfs, and 10 bound planetary mass companions. In this work, we investigate the contribution of methodological systematics to scatter in the Ṁ–M relation and compare brown dwarfs to stars. In our analysis, we rederive all quantities using self-consistent models, distances, and empirical line flux to accretion luminosity scaling relations to reduce methodological systematics. This treatment decreases the original 1σ scatter in the logṀ–logM relation by ∼17%, suggesting that it makes only a small contribution to the dispersion. The CASPAR rederived values are best fit by Ṁ∝M2.02±0.06 from 10 M J to 2 M ⊙, confirming previous results. However, we argue that the brown-dwarf and stellar populations are better described separately and by accounting for both mass and age. Therefore, we derive separate age-dependent Ṁ–M relations for these regions and find a steepening in the brown-dwarf Ṁ–M slope with age. Within this mass regime, the scatter decreases from 1.36 dex to 0.94 dex, a change of ∼44%. This result highlights the significant role that evolution plays in the overall spread of accretion rates, and suggests that brown dwarfs evolve faster than stars, potentially as a result of different accretion mechanisms.


INTRODUCTION
In the classical picture of star formation, molecular cloud cores collapse under gravity to form new stars.As the cores collapse, rotational velocity and conservation of angular momentum causes infalling material to settle into a circumstellar disk (Hartmann 1998).These primordial disks have been found to have a lifetime of ∼1-10 Myr (e.g.Strom et al. 1989;Armitage et al. 2003;Sicilia-Aguilar et al. 2006;Li & Xiao 2016), during which time they provide the material essential for both planet formation (Mamajek 2009) and stellar accretion (Hartmann et al. 2016).The evolution and dispersal of the disk unfolds through several processes.Planet formation occurs through core accretion of planetesimals in the inner few astronomical units (Safronov & Zvjagina 1969;Hayashi et al. 1985;Pollack et al. 1996) leading to terrestrial planet formation, while in the outer disk, core accretion, fragmentation, and instabilities are theorized to form giant gaseous protoplanets (Kuiper 1951;Cameron 1978;Boss 1997;Bate et al. 2002Bate et al. , 2003;;Johnson & Gammie 2003;Rafikov 2005;Cai et al. 2006).
Through the T Tauri phase, the disks are also actively accreting material onto the star, enabled by a combination of viscous accretion through the magnetorotational instability (Hawley & Balbus 1991) and/or MHD disk winds, depending on physical conditions in the disk (Pascucci et al. 2023;Lesur et al. 2023).Additionally strong near-UV, far-UV, and/or X-ray radiation from the star and its accretion shock can heat the gaseous disk surface leading to thermally driven photoevaporative winds beyond the gravitational radius, which likely account for the final clearing of the disk gas (Alexander et al. 2014, and references therein).
Within a few stellar radii of the star (traditionally assumed to be ∼ 5 R ⊙ ; Gullbring et al. 1998), the disk is interrupted by strong stellar magnetic fields and disk material flows to the stellar surface along accretion columns following magnetic field lines resulting in strong shocks on the stellar surface.Resultant emission from the accretion onto the star includes broad emission lines in the free-falling magnetospheric flows (Muzerolle et al. 2001;Hartmann et al. 2016), though Dupree et al. (2014) suggested that in the case of hydrogen, the broad lines could be formed in a turbulent postshock region, and forbidden lines from accretion shocks and winds (Hartigan et al. 1995).When the gas shocks on the stellar photosphere, the already fully ionized gas heats to 10 6 K.The optically thin preshock region is seen primarily as Balmer continuum excess (Valenti et al. 1993;Gullbring et al. 1998Gullbring et al. , 2000;;Calvet & Gullbring 1998), while the optically thick post-shock region emits Paschen continuum excess1 (Calvet & Gullbring 1998).These sources of excess continuum emission result in veiling of photospheric absorption lines.Accretion rates measured from both optically thin and thick shock regions can be inferred from the total excess luminosity produced by accretion; however, there is currently no direct method to measure the mass accretion rate ( Ṁ ) from the emission produced in the accretion flows.Instead, a scaling relation must be applied to relate a single emission line luminosity to a mass accretion rate.
Comprehensive multiwavelength studies have found that Ṁ decreases with decreasing stellar mass (M ) (Muzerolle et al. 2003;Calvet et al. 2004;Herczeg & Hillenbrand 2008;Alcalá et al. 2014;Manara et al. 2017b), following a power law of Ṁ ∝ M 2 in the stellar regime.This mass accretion rate−mass ( Ṁ − M ) relation has been assumed to extend from the stellar to the substel-lar (M ≤ 0.075M ⊙ ) regime with no variation in slope (e.g., Muzerolle et al. 2003;Mohanty et al. 2005;Muzerolle et al. 2005), though some studies suggest a break to a steeper relation around 0.2 M ⊙ for older star forming regions (Alcalá et al. 2017;Manara et al. 2017a).Additionally, at all masses, there is significant 1-2 dex scatter in accretion rates.Within the stellar regime, Manara et al. (2023) assert that the majority of this scatter results from physical variation and not observational uncertainty.
Various studies have looked at possible physical mechanisms responsible for this dispersion in the Ṁ − M relation.These include: intrinsic variability and the decrease of Ṁ with age (e.g.Natta et al. 2006;Costigan et al. 2014;Venuti et al. 2014;Hartmann et al. 2016), differences in the properties of star-forming cores (e.g.Dullemond et al. 2006;Clarke & Pringle 2006;Ercolano et al. 2014), competition among accretion mechanisms (viscosity or gravitational instability, e.g.Vorobyov & Basu 2008, 2009;DeSouza & Basu 2017), multiplicity (e.g.Zagaria et al. 2022) and, for PMCs, differences in the instability threshold and reservoir for accretion as a result of disk fragmentation (e.g.Stamatellos & Herczeg 2015).However, systematic studies of large numbers of accreting substellar objects are lacking, and it is not clear if this proposed explanation holds in this low mass regime.
Additionally, scaling relationships between line emission and accretion luminosity have been empirically developed and calibrated for stars for a wide variety of emission lines (Alcalá et al. 2014(Alcalá et al. , 2017;;Rigliaco et al. 2011;Natta et al. 2004).It is not clear to what extent empirical scaling relations are valid in the substellar regime, where potential differences in accretion (magnetospheric, planetary shock) and physical parameters (energy loss, magnetic field strength, temperature of accreting gas, gravitational potential, disk mass), could alter the relationship between line luminosity and mass accretion rate (Thanathibodee et al. 2019;Aoyama et al. 2020).In order to understand the origin and accretion of both bound planetary mass companions (PMCs; which include protoplanet candidates (e.g., PDS 70b and c, Delorme 1 (AB)b, LkCa 15b; Haffert et al. 2019;Wagner et al. 2018;Sallum et al. 2015;Eriksson et al. 2020;Ringqvist et al. 2023;Betti et al. 2022) and brown dwarf (BD) companions, objects in bound orbits around a higher mass host, (which we define as bound objects below M < 30 M J ; c.f. Martinez & Kraus 2019)), and young BDs (hereafter, all considered substellar objects), we first must characterize the physical (e.g., variability, age), and systematic (e.g., accretion rate tracer, evo-lutionary model) properties that affect accretion rate estimates.
The aim of this paper is threefold: a) to provide the largest compilation to-date of brown dwarf and protoplanet accretion rates derived under a uniform methodology, b) investigate methodological differences in the scatter of the Ṁ − M relation between stars and brown dwarfs, and c) determine whether the statistics of Ṁ measurements suggest accretion differences between these mass populations.In Section 2, we give an overview of the Comprehensive Archive of Substellar and Planetary Accretion Rates (CASPAR).In Section 3, we rederive object properties (e.g.mass, distance, temperature) in a consistent manner.In Section 4, we detail the technique we applied to derive linear fits for the Ṁ − M relation.We also present an updated Ṁ − M relation and discuss the role of methodology in producing scatter in Section 5.In Section 6, we quantify the contribution of various drivers producing the physical scatter observed in the overall Ṁ − M relation, and in Section 7, we focus on the BD population.In Section 8, we discuss how these phenomena affect our interpretation of accretion in the substellar regime.Results are then summarized in Section 9.

OVERVIEW OF THE DATABASE
We have assembled accretion rates for young planetary mass companions, brown dwarfs, and Classical T Tauri stars (CTTS) from large surveys of accreting objects, as well as individual object papers.This database consists of two parts: a compilation of published accretion properties, unmodified from their source publications (hereafter the Literature Database), and a unified re-derivation of accretion properties from these studies, CASPAR 2 .
We have focused on collecting properties for known accreting substellar objects.Within the Literature Database, 86 objects are considered substellar (below the hydrogen burning limit, M < 0.075 M ⊙ , e.g.Mohanty et al. 2005), of which 10 are PMCs (5 are protoplanets, 5 are M < 30 M J brown dwarfs).The database also includes a substantial compilation of CTTS accretion rates for stars later than G spectral type.We exclude Herbig stars, as detailed accretion census papers for this population already exist (e.g., Guzmán-Díaz et al. 2021;Vioque et al. 2022) and we are particularly focused on substellar accretion.To-date, we have compiled data for 798 objects from 46 studies for a total of 1058 independent accretion measurements spanning 24 2 CASPAR is openly available on Zenodo: https://doi.org/10.5281/zenodo.8393054.
years, from 1998 to early 2022.The list of references is given in Table 1.Sky positions for all objects are shown in Fig. 1.As many of the objects are in associations and clusters with small angular scales, they appear as a single point.We show zoomed-in views of six of the regions as insets.
As this sample is compiled from many individual studies, it is an incomplete survey of nearby objects in both mass and volume.In Fig. 2, we show the mass function of our sample with uniformly rederived masses (as discussed in Sec.3), colored by age and compared to the Chabrier (2005) initial mass function (IMF) and Kirkpatrick et al. (2021) field brown dwarf mass function, all normalized so the integral over mass is one.The mass distribution of objects in CASPAR is consistent with the IMF, though there is a difference at 0.1 − 0.3 M ⊙ , likely due to undersampling of objects at all ages in this mass range.A majority of CASPAR objects are young (< 3 Myr), as expected since disk fraction rapidly declines after 2.5 − 3 Myr (Mamajek 2009).
From the Chabrier (2005) IMF, if all of the star forming regions follow the same IMF, we would expect ∼ 20% more substellar objects.These missing objects have either a) not been surveyed or b) were initially missed when compiling CASPAR.For example, when we compare CASPAR objects in ρ Ophiuchus to the census from Esplin & Luhman (2020) (complete up to spectral type earlier than M6), we find (as shown in Fig. 3) that CAS-PAR includes 86% of all known ρ Ophiuchus substellar objects (32/37; M < 0.075 M ⊙ ; ≤ M5.5) with optically thick disks (which we use as a proxy for potential to be accreting).
Broad classes of object properties included in the database are summarized in Table 2, with individual column headers listed in Appendix A. Both the Literature Database and CASPAR have identical columns.Each accretion rate is assigned its own row and a unique number identifier, identical between the two databases.Therefore, an object observed at multiple epochs has multiple unique number identifiers.Each object is also identified with a unique name.Mass accretion rates have been measured from four broad accretion diagnostics families, namely: continuum excess, line luminosity, Hα photometric luminosity, and line profile.In Appendix B, we discuss in more detail the process of compiling the literature database, kinematic, photometric, and age information for each object.
We show all literature database accretion rates as a function of mass colored by their accretion diagnostic in Fig. 4. Overall, we see that the accretion rates vary by 2 orders of magnitude while still following an empirical power law relationship between accretion rate and mass,   (Doi et al. 2015;Takita et al. 2015).Insets: Enlarged views of several nearby star forming regions.Ṁ ∝ M 2.15 (black line) similar to other empirically derived relations (Muzerolle et al. 2003(Muzerolle et al. , 2005;;Mohanty et al. 2005;Alcalá et al. 2017).However, as shown in the botton panel of Fig. 4, we do find systematic offsets and variable slopes in the Ṁ ∝ M x relation when we fit by accretion diagnostic (see Section 4 for fit details), which will be discussed in Section 5.2.

UNIFIED DERIVATION OF QUANTITIES
The studies in the literature compilation come from a variety of instruments, analysis pipelines, and accretion tracers, which likely contributes to the wide dispersion of accretion rates at a single mass (i.e., the scatter).Additionally, the dependence of mass and radius estimates on the application of a variety of different evolutionary models and spectral fitting tools also introduces scatter (see gray dashed lines in the top panel of Fig. 4).In order to remove these effects, we first investigate the dispersion introduced by methodology by re-estimating object and accretion parameters under a unified set of assumptions and by comparing them to Literature Database values.Estimates of PMC spectral types and masses are highly uncertain and have been estimated from a variety of methods including: kinematics, orbital fitting, and spectral fitting.Due to the larger uncertainties in deriving accurate masses for the lowest mass objects, we focus here only on stars and BDs.Updating the PMC entries so that their masses are uniform and comparable to the full CASPAR sample will be the subject of future work.As a result, in the remainder of this study, we use the Literature Database mass and accretion rate estimates for PMCs in our plots and calculations.While we fit this PMC population and report it here, we do not attempt to infer any trends, as these accretion rates and masses are not rederived under a unified system.
In this work we update and unify literature values by performing the following modifications: • adopt Gaia DR2/EDR3 (Gaia Collaboration et al. 2018, 2021) distances when available (N=604, Bailer-Jones et al. 2018, 2021) • adopt single ages for each star forming region (or subregion, where available; see Section B.1), • adopt the spectral type -temperature conversion from Herczeg & Hillenbrand (2014) • extract mass, luminosity, radius, and surface gravity using the MIST MESA models (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015;;Dotter 2016;Choi et al. 2016) K spectral types, and Baraffe et al. (2015) evolutionary models for all others • calculate accretion luminosities using the Alcalá et al. (2017) L acc − L line scaling relationships.For excess continuum based accretion luminosity estimates, we scale by d 2 from Gaia DR2/EDR3.
• From accretion luminosity, calculate the mass accretion rate as where Ṁ is the mass accretion rate, R is the stellar radius, R in is the truncation radius, which we assume to be 5 R ⊙ (Gullbring et al. 1998), M is the stellar mass, G is the gravitational constant, and L acc is the accretion luminosity.
We describe the unified methodology for deriving each parameter in full detail in Appendix C. We show comparisons of rederived CASPAR parameters and the literature parameters in Fig. 5.For all rederived quantities, 50-89% are within a factor of two of the literature value, indicative of the relatively small effects of these updates.The change from individual or previous estimates of star forming region ages to a uniform set of ages result in the large dispersion seen between the original and rederived values.
We find 60 objects that were originally considered low mass stars with masses > 0.075 M ⊙ are now classified as brown dwarfs with masses below the hydrogen burning limit (HBL; conversely six of the brown dwarfs are now classified as stars).We highlight these in the lower right quadrant in Fig. 5.Of those, ∼ 20 are within 1σ of the HBL.These mass shifts result from using the Herczeg & Hillenbrand (2014) spectral type to temperature conversion and updated ages that were not used by the original references.Of the 16 studies where objects traverse below the HBL, 11 were published before the Herczeg & Hillenbrand (2014) spectral type-temperature conversion; therefore, their methods of calculating temperature and spectral types differed (with the majority using the conversion of Luhman et al. (2003)) leading to differences when reevaluating them.The four published after Herczeg & Hillenbrand (2014) utilized their own spectral fitting and used the spectral type-temperature conversion of Luhman et al. (2003) for the M dwarfs in their sample.When we look at the temperature -spectral type conversions from Herczeg & Hillenbrand (2014) and Luhman et al. (2003), we find that they diverge ∼ 100 − 150 K for M dwarfs, a change of ∼ 0.06 M ⊙ , with the Herczeg & Hillenbrand (2014) temperatures found to be lower for the same spectral type.This results in the decrease in mass seen in CAS-PAR for these objects.The objects are from 12 different star forming regions and from 16 different original references, indicating no preferential biases in deriving masses.Under this uniform derivation, we find 658 stars, 130 brown dwarfs, and 10 PMCs (see Table 3).The residuals between CASPAR accretion rates and their literature-derived values are shown in Fig. 6 for the full sample and each accretion diagnostic.Overall, we find that 1σ standard deviation in the residuals for all rederived accretion rates is 0.38 dex over the range of [−3.8, 1.15] dex (mean = −0.05dex).We find that 909 (88%) of the accretion rates change by less than 0.5 dex, indicating the vast majority of the objects have not markedly changed from their literature value.We find that the most disparate changes in accretion rate measurements are by Hα photometric luminosity and continuum excess, with CASPAR measurements larger than previously calculated.When we look at those accretion rates calculated from excess continuum, we find a median difference between the rederived and literature radius/mass ratio of 0.3 dex.As these accretion rate measurements are dependent on R/M ( Ṁ ∝ R/M ), this in turn leads to a median difference between the rederived and literature Ṁ of ∼ 0.1 dex, indicating the majority of accretion rates derived from excess continuum are higher than originally estimated.

LINEAR FITTING TECHNIQUE
With CASPAR, we can better investigate the causes of scatter in accretion rates.In the following sections, we discuss the relationship between accretion rate and parameters such as mass and age.Unless otherwise stated, these relationships are derived using the hierachical Bayesian linear regression routine linmix (Kelly 2007) in Python3 to determine the slope, intercept, and intrinsic scatter around relations of the form: y = mx + b. linmix allows for heteroscedastic and correlated measurement errors, and takes upper limits into account.It assumes that x and y variables are drawn from a 2D Gaussian distribution, and the covariance matrix is composed of the uncertainties in x and y.Calculated regression coefficients and uncertainties are derived from the posterior probability distributions of model parameters computed using Markov Chain Monte Carlo (MCMC).Not every object in the Literature Database, and therefore CASPAR, has a reported Ṁ uncertainty.For objects with no literature uncertainty, we assume the average Ṁ uncertainty from all reported and rederived measurements in CASPAR (∼ 0.36 dex) in the fit.
For each fit performed, we recover the posterior distribution of the slopes and intercepts and calculate best fit parameters from the median and 1σ uncertainties.
Additionally, we calculate the Pearson correlation coefficient, R, between the x and y parameters.Best fit coefficients are recorded in Table 4.
In Appendix D, we discuss additional linear regression methods that were investigated, such as weighted least squares, ordinary least squares bisector, and orthogonal distance regression.We find that including upper limits into the fits does not significantly affect the best-fit co-  For masses, we show the hydrogen burning limit at M = 0.075 M⊙ with dashed black lines; under a uniform methodology, objects in the lower right quadrant decreased in mass from stellar to substellar, while in the upper left, the reverse occurred.We find 50−89% of CASPAR parameters are within a factor of two of the literature value.efficients.However, the inclusion of x axis uncertainties can significantly affect the resulting best-fit coefficients.For fits that do not follow y = mx + b (such as y = e x ) and therefore cannot be fit with linmix, we utilize the orthogonal distance regression code scipy.odr,which allows for fitting of non-linear functional forms while taking into account both x and y measurement uncertainties, but not upper limits (though this should not greatly affect the fit as discussed above).

EFFECT OF METHODOLOGY ON ACCRETION RATE SCATTER
CASPAR spans a wide range of masses from ∼10 M J to ∼ 2 M ⊙ and compiles accretion rates calculated from four diagnostics.We first investigate a) the extent to which a uniform derivation reduces the extensive scatter in Ṁ , and b) whether or not the accretion diagnostic influences this scatter.

Scatter in Accretion Rate after Re-derivation
The relationship between mass accretion rate and mass follows a power law of the form Ṁ ∝ M α , where previous measurements of α range from 1.0-2.8(e.g.Calvet et al. 2004;Natta et al. 2004;Mohanty et al. 2005;Muzerolle et al. 2005;Herczeg & Hillenbrand 2008;Zhou et al. 2014;Hartmann et al. 2016).Typical dispersions are 1-2 dex (Manara et al. 2023, and references therein).The top panel of Fig. 4 shows the best fit to the Literature Database as log Ṁ = 2.16(±0.08)log M − 8.03(±0.06),with a 1σ dispersion of 1.00 dex and correlation of 0.73.The wide range of previously measured stellar slopes is consistent with our Literature Database fit, and is likely driven by uncertainties in the masses and differences in sample ages and sizes (Hartmann et al. 2016).
Therefore, we first affirm that the rederived values in CASPAR are consistent with previous slope estimates and discuss how the scatter changes for a larger sample and under a uniform methodology.We show the best fit model for CASPAR Ṁ − M in Fig. 7 and find a linear trend of: with a 1σ dispersion of 0.85 dex and a strong correlation coefficient of 0.76.CASPAR was based on previous studies; therefore, a consistent slope with literature estimates is unsurprising.A uniform derivation reduces the 1σ dispersion by 0.15 dex, indicating that methodology accounts for only 17% of the original scatter.This implies that the remaining scatter is due to underlying physical mechanisms rather than methodological systematics, and that methodology (e.g., object mass and radius estimation technique, varied distance references, etc.), while accounting for some of the dispersion in the scatter, cannot explain all of it, especially as large scatter is seen in uniform star forming region surveys (e.g.Manara et al. 2015;Alcalá et al. 2017).
In Fig. 8, we also compare the best fit Ṁ −M relations from Zhou et al. (2014) and Hartmann et al. (2016) to CASPAR.Though we do not expect significant variation between our fit and previous estimates, the samples from Zhou et al. (2014) and Hartmann et al. (2016) were primarily composed of stars between 0.1 and 1 M ⊙ , while our fit includes objects down to 0.01 M ⊙ .Therefore, we investigate whether the substellar population is different than previous stellar fits.
We compare the best fit relation of (a) the CASPAR total sample and (b) stellar sample to the fits of Hartmann et al. (2016) and Zhou et al. (2014).We find that the Zhou et al. (2014)  and CASPAR single population fit shows a positive skew only for PMCs (due to the stellar population dominating the fit).In the histograms of the residuals in Fig. 8, we indicate the center of the distribution with a narrow gray line.We find that while the four fits overlaps in the stellar regime, they deviate at substellar masses, with our fits accounting for most of the scatter, within fit uncertainties.

Accretion Diagnostic Systematics
To assess the consistency among various observational methods, we assume that all methods should produce consistent Ṁ values of estimating accretion rates.In Fig. 9a, we show CASPAR Ṁ − M statistics separated by accretion diagnostic.We bin accretion rates by mass to facilitate comparison.To do this, we assume every detection is a Gaussian probability density distribution (PDF) with a mean of the accretion rate and standard deviation of the uncertainty.For non-detections, we assume a half normal distribution with the cut-off at the upper limit value.We then ran a Monte Carlo simulation drawing random values from the PDFs within each mass bin and took the median of the values.We find consistent median accretion rates for each mass bin across the four diagnostics, with the medians closely following the best fit line found for continuum excess.Fig. 9b shows the fit residuals.Though within uncertainties, Hα photometric luminosities produce higher (∼1 dex) Ṁ s, while line profiles produce lower (∼1 dex) ones.However, the majority of the Hα photometric measurements are for objects with ages < 1 Myr, biasing the results (see Section 6.1).Overall, line luminosity and continuum excess methods produce the smallest residuals in the stellar regime (residuals < 0. diagnostics.As there are only 10 PMCs, in which the masses and accretion rates are not uniformly derived, we cannot conduct a similar analysis, though we show them in Fig. 9.We find line profile derived accretion rates lie consistently below the best fit line for the substellar mass objects, while continuum excess is consistent to within ∼ 0.21 dex across all masses.For line luminosities, we find relatively small residuals (up to 0.8 dex) in the substellar regime.Though within uncertainties, line luminosity and line profile residuals trend upward from the stellar to substellar regime, while this is not seen for the accretion rates derived from the continuum excess tracers.The standard deviation of the residuals between line luminosity and excess continuum derived rates increase from 0.1 dex at 0.5 M ⊙ to 0.65 dex at 0.03 M ⊙ .As line luminosities are dependent on scaling relations (derived from excess continuum) to compute accretion rates, we expect similar trends among excess continuum and line luminosity derived Ṁ s.While they are consistent within uncertainties, the average accretion rates do show differences at low masses, potentially a result of utilizing stellar scaling relations in the substellar regime.We discuss this hypothesis further in Section 7.
Accretion rates measured from line emission are known to have significant uncertainty, as scaling relations are empirically derived and many accretion tracing lines can also be produced by other physical processes, such as winds and chromospheric activity (Jayawardhana et al. 2003;White & Basri 2003).Several criteria have been established to try to separate accretion and chromospheric activity using Hα equivalent widths (< 200 km/s; Jayawardhana et al. 2003;White & Basri 2003) and L acc /L ratio as a function of temperature (and spectral type) for emission lines (−3.19 ± 0.15 for M6 dwarfs; Manara et al. 2017b).
In Fig. 10, we group line fluxes by wavelength, namely: a) Balmer series, b) infrared hydrogen series (Paschen, Brackett and Pfund), c) Helium i emission lines (He i λ 4026, He i λ 4471, He i λ 4713, He i λ 5016, He i λ 5876, He i λ 6678, He i λ 7065), and d) Calcium ii emission lines (Ca ii K, Ca ii H, Ca ii λ 8498, Ca ii λ 8542, Ca ii λ 8662), in order to analyze trends among them.We first find best fit Ṁ − M relations for each line flux group (given in Table 4) and compare these fits to the CASPAR best fit relation.We use the Akaike's information criterion (AIC) to access the performance of the linear fits in explaining the variation in the data4 Since our two fits are independent (i.e., non-nested) with the same number of parameters, this criterion estimates how well the model reproduces the data from the maximum likelihood.
If the AIC values for two fits are within 10%, we consider the fits to be comparable.If the fit to a line flux group is not significantly more descriptive of the variance than the overall best fit relation, we conclude that offsets by method do not contribute to the scatter.In the stellar regime, we find the AICs for each line flux group are all within 10%.
In the substellar regime, the best fit line to the infrared hydrogen and He i line measurements is offset from the overall best fit with a percent difference between AIC ∼ 180%.Additionally, at the deuterium burning limit, the NIR best fit line is offset from the CASPAR best fit by 0.9 dex (compared to < 0.4 dex for other lines).This could indicate that the NIR line flux scaling relations overestimate accretion rates in the substellar regime (i.e.their scaling relations are different), or this may result from small number statistics.
Finally, we estimate the extent to which accretion rates derived from Hα 10% widths, Hα line profiles, and UV excess differ in order to probe their effect on scatter.When accretion rates are derived from different emission lines or continuum at the same epoch, they should not be subject to intrinsic accretion variability.This makes contemporaneous measurements an excellent probe of systematics.While Ṁ s from Hα luminosity rely on L line − L acc scaling relations, the Ṁ −Hα 10% width scaling relation of Natta et al. (2004) relates the 10% width directly to the accretion rate.However, Alcalá et al. (2014) found that it can underestimate accretion rates by almost 0.6 dex for widths < 400 km/s in Lupus (corresponding to M < 0.3 M ⊙ ) compared to excess continuum.We find a similar result when we compare the CASPAR Hα 10% widths to accretion rates derived from excess continuum.
In Fig. 11, we show the residuals in CASPAR accretion rates derived from excess continuum (top) and Hα 10% width (bottom) compared to simultaneous measurements from Hα luminosity as a function of mass.In the stellar regime, accretion rates derived from these three quantities do not significantly differ within 1 dex.In the substellar regime, accretion rates derived from Hα luminosity are systematically high when compared to  Line luminosities (including Hα photometric) show (significant) deviation from the continuum excess with decreasing mass.Hα photometric luminosities appears significantly offset from the CASPAR continuum excess best fit; however, the majority of the stellar accretion rates are from objects with ages < 1 Myr, which biases the high mass results (see Sec 6.1).
excess continuum, which could be indicative of an overestimation of accretion luminosity for the lowest mass BDs and PMCs.We find a large offset in accretion rates derived from Hα 10% width compared to Hα luminosity.Several other processes, including chromospheric activity, outflows, and hotspots, contribute to Hα emission, potentially inflating its width and therefore increasing its inferred accretion rate.Low mass substellar accreting objects can have line widths below the traditional 200 km/s threshold for accretion.Such observations have led previous work (e.g.Alcalá et al. 2014Alcalá et al. , 2017) ) to discourage the use of Hα 10% width as an accretion tracer.We confirm this offset between Hα 10% width and Hα luminosity within the substellar regime and find it can lead to a difference of almost 2 dex in calculated accretion rate near the deuterium burning limit, producing vastly overestimated accretion rates.

Summary
We find that methodological differences in estimates of mass accretion rates, such as differences in evolutionary models and estimated distances (used to scale accretion luminosities), account for only ∼ 17% of the scatter in the Ṁ − M relation, indicating that the remaining scatter is from either observed accretion diagnostic systematics or physical differences (e.g., variability, disk mass, stellar mass, system age).
When we separate accretion rate estimates by accretion diagnostic, we find that accretion rates do not vary significantly (< 1 dex) in the stellar or substellar regime.However, we do find systematically higher accretion rates for Ṁ s derived from NIR line luminosities (0.9 dex offset in best fit line at the deuterium burning limit; Fig. 10b) and Hα luminosity relative to continuum estimates (top panel of Fig. 11).We will discuss the physical and diagnostic drivers of this scatter in the following sections.

DRIVERS OF PHYSICAL SCATTER IN MASS ACCRETION
As shown in the previous section, methodological systematics cannot fully explain the scatter in the Ṁ − M relation for either stellar or substellar objects.Recent work has suggested that multiplicity may affect accretion rates, with binaries accreting at a higher rate than isolated objects (Zagaria et al. 2022;Gangi et al. 2022).While CASPAR currently does not contain many objects in multiple systems (46/798), they appear consistent with isolated object accretion rates and do not have a significant effect on the Ṁ − M scatter.Below, we focus on the relationship between age and intrinsic vari- We find that the NIR line flux fit shows the most deviation in best fit; it is offset by 0.9 dex compared to the best linear fit at the deuterium burning limit (compared to < 0.4 in the other lines).
ability in the observed Ṁ − M scatter 5 , and what this may tell us about the evolution of accretion activity.

Variation in Accretion with Age
Circumstellar disk fraction in young star forming regions has been found to decrease exponentially with 5 Though disk mass has a known correlation with accretion rate (Manara et al. 2023, and references therein), and modeling work has suggested variations in accretion rates are due to differences in disk mass (Vorobyov & Basu 2009) (Hartmann 1998;Manara et al. 2016b;Mulders et al. 2017) are generally explained through a combination of viscous evolution (Lodato et al. 2017;Rosotti et al. 2017;Mulders et al. 2017), disk photoevaporation (Sellek et al. 2020) multiplicity (Zagaria et al. 2022).Accretion rates decrease with time as t α , where α = −1.6 to −1.2 (Hartmann et al. 2016, and references therein).This decrease has been attributed to viscous evolution, though observations of the POISSON sample found higher Ṁ than expected from pure viscous models (Antoniucci et al. 2014).By establishing a "uniform" estimate of ages, we can study the correlation of Ṁ with age, and its impact on Ṁ − M scatter.We exclude PMCs in this analysis as a) they have not been uniformly re-derived, and b) theoretically modeled accretion (Aoyama et al. 2018;Thanathibodee et al. 2019;Aoyama et al. 2020) and formation (Stamatellos & Herczeg 2015) mechanisms posit that stellar age trends could not hold for PMCs.We divide CASPAR into the following four age bins for this analysis: • ≤ 1 Myr: includes the Lagoon Nebula • 1 < t/Myr ≤ 3: includes Chamaeleon I, Taurus, and Lupus • 3 < t/Myr ≤ 8: includes IC 348 and Perseus • > 8 Myr: includes Upper Centaurus Lupus and η Chamaeleontis.
See Table 7 for full list of ages for each cluster and association within these broad groups.
In Fig. 12, we show accretion rate as a function of age, colored by stellar mass.Following Hartmann et al. (2016), we scale to the accretion rate to M * = 0.7 M ⊙ in order to remove the dependence on mass.The best fit to these data using linmix is:  2016) consisted of 148 objects whose accretion rates were computed from continuum excess and line emission with ages of 5 < log t/year < 6.This smaller sample size and age range could account for the difference in fit.
As a simple test of the effect of age, we begin by assuming that objects, no matter their age, share the same Ṁ ∼ M 2.02 slope.For older objects, the accretion rates should be lower resulting from a decrease in available disk material.We can model this simplified assumption as a decrease in the intercept of the log Ṁ − log M relation with increasing age.This allows us to quantify the extent to which age affects the dispersion in the Ṁ − M relation.
As shown in Fig. 13, we find that the intercept decreases by 1 dex from < 1 Myr to > 8 Myr.Subtracting age best-fits from each age population and over plotting the residuals gives an indication of the effect of age on overall scatter.If scatter in Ṁ for objects of the same mass is a result of age, then we expect the overall dispersion to decrease when age is accounted for in this way.However, we find from these residuals that the 1σ scatter only decreases by 0.06 dex (bottom panel of Fig. 13).While there is almost 1 dex decrease in best fit intercepts by age bin, the median scatter remains high.Therefore, as the dispersion within each age bin is roughly the same as the overall scatter, this normalization by age has an insignificant effect.
In the substellar regime, age has an even less effect on the overall scatter; we find the standard deviation of the residuals decreases by 7%.From Fig. 8, the residuals in the substellar regime are positively skewed, exacerbating this effect in the age Ṁ − M residuals, especially for the < 1 and (3-8] Myr regimes.This could indicate that substellar objects are following a different trend with mass or age.We investigate this hypothesis in Sec 7.

Multi-epoch variability
Previous surveys (Biazzo et al. 2014;Costigan et al. 2014;Venuti et al. 2014) have primarily focused on dayto-day intrinsic accretion variability utilizing one accretion diagnostic.This variability produces ∼0.4 dex dispersion, smaller than the 1-2 dex observed in the Ṁ − M relation.The variability is therefore not the dominant source of scatter in the Ṁ − M relation for a given mass.However, CASPAR contains objects that have been observed from many different tracers across months and years.In this section, we quantify how this longer timescale and methodological accretion variability affects the dispersion we see in the CASPAR Ṁ − M relation.
We first look at all sources of variability, including variation in accretion rates from different lines in the same epoch, and variation in accretion rate from the same line over multiple epochs.The median separation in measurements is 3 years, with a maximum separation of 21 years from Gullbring et al. (1998) to Alcalá et al. (2021).Overall, we find a spread of 0.63 dex in the stellar regime that includes several high variability outliers from multi-line measurements from a single epoch (see Fig. 14a), and that increases in the substellar regime to a median of 1.33 dex.We compute an independent two-sample student's t-test, and find a significant (t(134)= 2.08, p = 0.03) difference in variability between the stellar (mean = 1.00 dex, σ = 1.21 dex) and substellar (mean = 1.35 dex, σ = 0.78 dex) populations.This increase in variability could be indicative of either intrinsic variation resulting from rotation of sunspots or accretion flows in the line-of-sight or differences in diagnostics.
We also analyze two main sources of observed variation separately, namely: a) variation among accretion rates measured from multiple lines at one observational epoch ("methodological"; 241 objects, Fig. 14b), and b) variation among accretion rates for one line observed at multiple epochs ("intrinsic"; 67 objects, Fig. 14c).See Table 5 for the full breakdown of objects.We find both the median methodological and intrinsic variation in accretion rate is on average 0.86 dex, double the amount of variability found in previous surveys.The methodological variability in Ṁ determination is relatively consistent with mass, however, we do note some strongly variable outliers.We find these objects show small (< 1.5−2 dex) variability over multi-decade timespans using a single accretion tracer, but have large (> 3 dex) variation between accretion rates measured for different lines, with the largest outliers found in the Ophiuchus star forming region (3.5 − 5.5 dex).The object with the most extreme multi-line variability is DO Tau, which has a relatively stable accretion rate of 10 −8 M ⊙ yr −1 .However, Alcalá et al. ( 2021) found an Hα line flux of 1.92(±0.07)e−17erg s −1 cm −2 , corresponding to an accretion rate of 6.45 × 10 −15 M ⊙ yr −1 in CASPAR, a > 7 dex difference in Ṁ .
For intrinsic variation, we find Paβ and Ca ii K produce the highest amounts of variability (1.26 dex and 1.23 dex, respectively), though small samples sizes affect the measured ranges of variability.For Paβ, we find the median variability increases from stellar (0.33 dex) to substellar (0.57 dex) regimes, though this difference is not statistically significant (t = −1.42,p = 0.16), and as shown by Claes et al. (2022), measurements found from line luminosity can underestimate the variability for CTTS.
Though no significant mass trends were seen for either methodological or intrinsic variation, considered separately, when all sources of variation are considered together, we find a significant increase in the variability of Ṁ in the substellar regime.As we show in Fig. 9b and 10, while excess continuum accretion measures remain relatively constant with mass, accretion rates derived from line fluxes, particularly Paβ, appear to deviate more from the Ṁ − M relation as mass decreases.
For objects derived from both excess continuum and line flux, there could be as much as 0.9 dex difference with Ṁ estimates.This may point to physical differences in accretion processes in the substellar population (e.g. more luminosity in certain emission lines compared to the expectation from stars) that are not accounted for in current L line ∼ L acc scaling relations.This could have a large impact on our interpretation of ongoing accre-tion and variability in the lowest mass objects (Betti et al., in prep).

ROLE OF ACCRETION ON SUBSTELLAR FORMATION
In Sections 5.2 and 6.1, we find offsets and skewed residuals in Ṁ − M relation for the substellar regime according to age, line tracers, as well as from the overall best linear fit, which might point to underlying differences in the Ṁ − M relations needed to describe the stellar and substellar populations.Therefore, we explore whether the accretion rates of brown dwarfs are statistically distinct from stars and whether this can be connected to their accretion or formation mechanisms.If BDs are accreting differently, this could appear as differences in their calculated accretion rates.
When we fit each population (stars, brown dwarfs, and planetary mass companions; hereafter mass populations) separately, we find three distinct relationships, shown in Fig. 15 and described by the following best-fit relations (hereafter three population fit).For the stellar population, we find: log Ṁ = 2.18(±0.11)log M − 7.97(±0.05), with a residual standard deviation of 0.74 dex.The accretion rate, Ṁ , is in M ⊙ yr −1 and M in M ⊙ .For the brown dwarfs: log Ṁ = 3.19(±0.57)log M − 6.54(±0.80), with a residual standard deviation of 1.36 dex.Finally, for the planetary mass companions, we find: log Ṁ = 0.25(±2.91)log M − 10.66(±5.64),( 6) with a residual standard deviation of 0.81 dex.The stellar fit (slope = 2.18) is similar to those found by other authors (e.g.Calvet et al. 2004;Natta et al. 2004 For BDs, the slope of the relation steepens to 3.19 in the substellar regime, while the PMCs are best modeled by a flat dependence with mass.In Appendix E, we compare the three population fits to the single population best fit described in Section 5.1, and find that separate fits are statistically favored over the single population fit with lower AIC statistics. This steeper slope ("knee") in the substellar/low mass star regime has only been seen observationally in older individual star forming systems (Alcalá et al. 2017;Manara et al. 2017a).As CASPAR includes multiple SFRs at different ages, this steeper BD slope could be a result of these older ages having a stronger effect on the best fit Ṁ − M relation for this mass regime.Therefore, in the following sections, we explore the effects of mass with age and mass population in explaining both the scatter and evolution of BD accretion.

Effect of Age on Brown Dwarf Properties
Following the process outlined in Section 6.1, we explore the effect of age in explaining the Ṁ − M relation for substellar objects.We combine the original < 1 Myr and (1-3] Myr populations due to small numbers within the < 1 Myr bin.As shown in Fig. 16, we find consistent decreases in the best-fit y intercept with age for the substellar regime (as well as the stellar regime).Residuals of these age fits show a significant decrease in dispersion as a result of fitting the BD and stellar populations separately; the 1σ standard deviation of the residuals for the BD population is 1.02 dex.We find a decrease of ∼33% from the 1.36 dex found in Section 7.
In the substellar regime, we find higher accretion rates at the HBL compared to what is expected from the stellar best-fit for each age.However with the steeper slope, by the deuterium burning limit (DBL), these relationships predict accretion rates just slightly lower than those extrapolated from the stellar best-fits.Higher mass BDs additionally appear to accrete faster at older ages compared to very low mass stars at similar ages, as seen at the HBL in Fig. 16.
In order to access the rate of decay of accretion for the substellar population compared to the stellar population, we extract log Ṁ at the HBL and DBL for the stellar and substellar fits.We find the following exponential best fits to these data: where Ṁ is in M ⊙ yr −1 .We show these fits at the HBL and DBL as a function of binned age in Fig. 17.The exponential trend in the stellar regime is similar to the exponential disk fraction decay timescale (τ = 2.5 Myr) found by Mamajek (2009), while the decay rate for BD accretion remains high compared to the τ = 3 Myr timescale for disk fraction decay (Mamajek 2009), with objects still accreting quickly at older ages.

Relationship between Accretion Rate and Mass, Age, & Mass Population
To test whether the relations among mass accretion, object mass, and age show evolutionary trends that could explain the (lack of) knee in the accretion-rate timescale, we fit the Ṁ − M relation with age (t) and mass as free parameters.More specifically, we fit the model for each mass population, where Ṁ is in As shown in Fig. 18, in the substellar regime, we find a clear mass and age dependence on accretion rate, with slope increasing with age, while this trend flattens in the stellar regime.In older systems, we see a systematic steepening in the Ṁ − M slope compared to the stellar regime.This is suggestive of a faster evolutionary timescale for accretion onto brown dwarfs (especially low mass BDs), as they accrete material at a relatively higher rate ( Ṁ ∼ 10 −11 ) at young ages depleting their disks quickly.Higher mass BDs and stars instead accrete Figure 16.top: CASPAR accretion rates vs mass colored by age with best linear fits for each age and mass population (brown dwarf, stars).We combined the brown dwarf ≤ 1 and (1 − 3] Myr populations to improve statistical number counts (brown pentagons).We extend the stellar fits into the BD regime indicated by dashed lines.bottom: Residuals for the best linear fits for each age and mass populations.The colors and markers are as in Fig. 13. for longer at a high accretion rate ( Ṁ ∼ 10 −7 − 10 −10 ), indicative of relatively slower disk depletion at younger ages.This confirms previous work in individual SFRs showing a shallower substellar trend for younger systems (Manara et al. 2015;Fiorellino et al. 2021) and a steeper trend in older systems (Alcalá et al. 2017;Manara et al. 2017a).
When we fit to both mass and age, the standard deviation of the best fit residuals for all objects is 0.78 dex, a 9% decrease from the single population best fit standard deviation (0.85 dex).Compared to the three population fits, the residual scatter decreases by 6% (0.71 dex from 0.75 dex) and 44% (0.94 dex from 1.36 dex) in the stellar and substellar regime.When we remove the outlier upper limits, these decrease by 17%, 22%, and 58%, for the total, stellar, and brown dwarf populations respectively, indicating both mass and age have a profound effect on the rate of accretion for BDs.This effect is also apparent when we also bin by age (as discussed above), and have slope as a free parameter (see By Age and Mass in Table 4).For the stellar fits, the slope becomes shallower with age, while the BD fits become steeper.We also see this trend when we examine individual star forming regions in Appendix F.

DISCUSSION
In Section 5.1, we find that uniformly deriving physical and accretion properties reduces the scatter in the CASPAR Ṁ − M relation by 17%, indicating that methodology plays a role in the observed scatter when comparing accretion rates across multiple detection and analysis techniques.However, since a large (∼1 dex) scatter in the Ṁ − M relationship has been seen even in surveys using one detection technique (e.g.Muzerolle et al. 2001;Herczeg & Hillenbrand 2008;Alcalá et al. 2017), methodology cannot fully explain the observed dispersion around the relation.
We first investigate the roles of accretion diagnostic, intrinsic variability, and age on the Ṁ − M scatter.In the stellar regime, the scatter is consistent among diagnostics, especially between line luminosity and continuum excess.This is expected as the L acc − L line relationships are derived for young stellar objects (Alcalá et al. 2014(Alcalá et al. , 2017;;Rigliaco et al. 2011).We do start to see an increase in the line luminosity scatter in the substellar regime, with the median Ṁ increasing ∼ 0.8 dex off the single population Ṁ − M relation (Fig. 9), likely driven by NIR emission (Fig. 10).Though we do see intrinsic accretion variability in Ṁ over multiple lines at the same epoch or over multiple epochs (∼ 0.8 dex), as only 241 and 67 objects in CASPAR are contributing to this variability, respectively, they are likely not driving   the scatter.Indeed, this spread is in line with previous surveys of accretion variability, and smaller than the total observed scatter (see Manara et al. 2023, and references therein).Finally, when we just consider age as a driving factor, we find that the Ṁ scatter in each age bin is roughly the same as the overall scatter (see Figs. 13 and 16) and normalizing by age has little effect.Instead, we posit that the measured accretion rates, dispersion, and variability behavior of BDs is distinct from the stellar regime.We find that three population (stars, BDs, PMCs) fits are statistically favored over fits to a single population for objects of all masses, and that fitting mass populations separately results in the greatest reduction in residuals.This is most clearly seen in the brown dwarf regime, where fitting for both mass and age reduces the residual scatter by 44-58%.When the substellar and stellar populations are fitted separately, age effects are compensated for, and a uniform methodology is applied to derive accretion rates, the total scatter in the Ṁ −M relation across all mass regimes decreases from 1.0 dex to 0.78 dex, a decrease of 28%.
The hypothesis of a different Ṁ − M relation in the substellar regime is not unique to this work.Vorobyov & Basu (2009) predicted a bi-modality in the Ṁ − M relation; in particular, they predicted a steepening at lower masses.This Ṁ − M bi-modality was verified observationally by Alcalá et al. (2017) and Manara et al. (2017a) in the Chamaeleon I (1.5 Myr) and Lupus (2.5 Myr) regions, with a break at M ∼ 0.2 M ⊙ .However, this break was not seen in other regions of similar or younger ages (probing the stellar and high mass substellar regime Manara et al. 2015;Fiorellino et al. 2021).Manara et al. (2023) suggests that it could be an evolutionary effect, wherein low mass stars accrete their disk mass at a faster rate during the late stages of formation.From Fig. 18, we see a similar trend, where the slope steepens with age for substellar masses, showing that this evolutionary effect is universal for BDs even in different star forming regions.
Theoretical studies of brown dwarfs suggest a variety of formation mechanisms, from more planetary processes such as disk fragmentation (e.g., Bate et al. 2002Bate et al. , 2003) ) to protostellar embryo ejection (e.g., Goodwin et al. 2005;Hubber & Whitworth 2005) and turbulent fragmentation (e.g., Kirk et al. 2006).Previous observational surveys of brown dwarfs have found that their disk properties follow trends similar to stars, with similar disk fractions (Luhman et al. 2005), correlations be-tween disk and stellar mass (Testi et al. 2016;Ward-Duong et al. 2018;Sanchis et al. 2021;Rilinger & Espaillat 2021), and an inverse correlation between disk mass and age (Rilinger & Espaillat 2021).The lowest mass BDs should have small disks; however, at the youngest ages, they appear to accrete material at similar rates to higher mass BDs with larger disks (within an order of magnitude; see Fig. 18).In other words, the slope of the Ṁ − M relation for the youngest brown dwarfs is shallow.The steepening of the relation in Ṁ − M with age could be a result of the lowest mass brown dwarfs accreting "too rapidly" at younger ages and depleting their disk material.
Combining both physical and systematic offsets in the BD regime, evidence from CASPAR points to a two-fold complication when deriving accretion rates for the lowest mass objects.First, as discussed above, age and mass play a significant role in accretion rates (see Fig. 18).Additionally, accretion rates derived from line luminosities are calculated using empirically derived L acc − L line scaling relationships for stars.However, if stellar and substellar objects follow different Ṁ − M relationships, there is no guarantee that the same L acc − L line relationships accurately represent single (or bound) brown dwarfs.By using stellar relationships in the substellar regime, we could be artificially overestimating both accretion luminosities and accretion rates derived from line luminosities.Systematic variations in accretion rate scalings may play a role in the apparent variation among diagnostics for substellar objects, with up to 0.8 dex difference in inferred accretion rates for certain tracers (e.g., line fluxes vs. excess continuum).Forthcoming work will explore this issue (Betti et al., in prep), which may have a profound impact on the interpretation of substellar accretion.

Planetary Mass Companion Population
In this section, we describe briefly trends seen in the PMC population.Though the sample of accreting PMCs is small and has not been rederived consistently in mass or accretion rate, we note that these objects appear to follow a much flatter and higher relation compared to BDs.Stamatellos & Herczeg (2015) modeled disk and accretion properties of bound companions formed via disk fragmentation, allowing a gravitationally unstable disk with a mass of 0.7 M ⊙ around a 0.7 M ⊙ star to grow until it started to fragment.They predict that the disks around PMCs are more massive than expected for objects of the same mass forming in isolation from a collapsing core.This is due to the fact that PMCs are forming within the larger circumstellar disk.Before they separate from the disk (become dynamically inde-pendent), they are able to accrete more gas than a BD whose only material reservoir is its natal core (Stamatellos & Herczeg 2015).Stamatellos & Herczeg (2015) predict no strong correlation between object mass and disk mass under this scenario, leading to a flatter slope in Ṁ vs. M .The difference that we observe between high PMC accretion rates and low accretion rates for brown dwarfs of similar mass could be either: a) an observational bias in observed PMCs, such that objects with accretion rates below 10 −12 M ⊙ yr −1 have simply gone undetected, b) PMCs may be fundamentally different from BDs, for example by forming via a disk fragmentation-like mechanism as opposed to core collapse, or c) using incorrect scaling relations pushes our Ṁ estimates of all substellar objects up or increases the dispersion in the relation.
Recent theoretical work by Aoyama et al. (2018), Aoyama et al. (2020), andMarleau et al. (2022) also show that bound PMCs may have a higher fraction of line emission contributing to their total accretion luminosity than accreting stars.These models predict significantly larger (∼2 dex) accretion rates for planetarymass objects than those derived from stellar magnetospheric empirical relations (Alcalá et al. 2014(Alcalá et al. , 2017)), which would drive these Ṁ s to even higher rates.As BDs are traditionally used to place the accretion of PMCs in context, more care will have to be taken in calculating and interpreting accretion signatures from PMCs.

CONCLUSION
In this paper, we introduce CASPAR, the Comprehensive Archive of Substellar and Planetary Accretion Rates, the largest database of substellar and planetary accretion rates to-date.The physical and accretion properties in CASPAR have been rederived under a consistent methodology (Gaia distances, consistent evolutionary models and scaling relations).The goal of this effort was to investigate the contribution of systematic offsets among methods to overall scatter in the Ṁ − M relation.Using the rederived database, we investigate the dispersion in the Ṁ − M relation, and the physical and systematic processes that contribute to it.We also explore variation among stars, brown dwarfs and planetary mass companion populations.We find: • Rederiving all physical and accretion properties using the same methodology decreases the 1.04 dex of scatter about the single population Ṁ − M fit by 17%.The best single population linear fit for CASPAR, log Ṁ = 2.02 log M − 8.02, is consistent with previous estimates from smaller samples, suggesting that methodological differences in deriva-tion play a small role in the slope or scatter of the Ṁ − M relation.
• The choice of accretion diagnostic additionally contributes to the overall scatter at substellar masses, with estimates from line luminosities leading to an average of 0.8 dex variation for a single object.Within the stellar regime, accretion rates are consistent among tracers.Unlike line luminosity, excess continuum derived estimates are consistent to within 0.21 dex of the overall best linear fit across both substellar and stellar mass regimes (0.1 ≲ M/M ⊙ ≲ 2).
• We also find consistent multi-epoch and multitracer variability of ∼ 0.6 dex in the stellar regime, consistent with previous estimations.This variability increases to 1.33 dex in the substellar regime.We posit that this increase is due to either higher variability seen at lower masses or stellar scaling relations being invalid in substellar regimes, leading to offsets in derived accretion rates.
• We investigate the effect of age on dispersion around the relation and find a 1 dex decrease in the best fit intercept between ∼ 1 − 10 Myr.However, the scatter within age bins is ∼ 1 dex, leading to little change in the overall scatter compared to the scatter from the overall best fit residuals.
• We argue that the majority of the scatter can be explained by modeling the star, brown dwarf, and PMC populations by separate Ṁ − M relations, accounting for both mass and age.We find the brown dwarf Ṁ − M scatter decreases by 44% as a result (58% when upper limit outliers are excluded from the residuals).Additionally, we show that the BD Ṁ − M relation steepens with age, while the stellar relation flattens.
We posit that there is a two-fold issue when deriving accretion rates for low mass objects.First, accretion rates are expected to depend much more on age and mass than in the stellar regime.Secondly, accretion rates derived from stellar scaling relations likely overestimate BD accretion rates, contributing to the overall scatter in this population.
• Bound planetary companions seem to follow a flatter Ṁ −M relation compared to brown dwarfs and stars.This may be a result of differences in either their formation or accretion paradigms.Accretion measurements for a larger population of PMCs and individual modeling of these systems will help reveal the underlying physics governing them.
CASPAR is an evolving database and with future/ongoing surveys (e.g., Pittman et al. 2022;Gangi et al. 2022), protoplanet detections (e.g., Ringqvist et al. 2023), and derived scaling relations (e.g., Marleau & Aoyama 2022), it will continue to be updated.All updates and additions will be found on Zenodo at the following: https://doi.org/10.5281/zenodo.8393054.Suggestions for additions to CASPAR can be made to the lead author.

A. CASPAR COLUMN NAMES
In Table 6, we give all columns within the database and their description.These are identical between the Literature Database and CASPAR.For each object, we compiled kinematic information (Right Ascension, Declination, parallax, distance, proper motion, and radial velocity) from Gaia DR2 (Gaia Collaboration et al. 2016a, 2018) and Gaia EDR3 (Gaia Collaboration et al. 2016b, 2021), in case an ob-ject was not observed in one of the data releases.We queried the Gaia archives in order to find the Gaia object associated with each database entry.Of the 793 unique objects, we retrieved kinematic information for 670 that have Gaia observations (of the substellar object, 65/87 have Gaia observations).If the parallax is considered reliable (e.g., its parallax error is less than 25% and parallax > 0.167 mas; Huang et al. 2022), we use the geometric distances found for DR2 and EDR3 by Bailer-Jones et al. (2018) andBailer-Jones et al. (2021), respectively.Near infrared (NIR) photometry for each object is compiled from 2MASS, while the accretion literature reference Hα photometry is included for objects whose accretion rates are measured from this photometry.
For objects with D < 150 pc, we ran their kinematic information through the Banyan Σ tool (Gagné et al. 2018) in order to determine the most likely host association.If the membership probability is > 70%, we assume membership in that association and assign the object the corresponding age from Table 7.For those objects D > 150 pc or whose membership probabilities were lower than 70%, we searched the literature for population census studies that may have determined the membership in an association or a cluster.Finally, for those objects not in Banyan Σ or census papers, we determine if a) it is a known field star, or b) if its kinematic information is close to an association.If the latter, we assign it to the closest association with the caveat that the association is only assumed, indicated by a "*" next to the assigned association in the database.
In order to estimate object ages, we compile a list of the star forming regions (SFR), clusters, associations, and clouds associated with each object.Star formation is not instantaneous within molecular clouds, but occurs on individual, sub-group, and regional levels, leading to gradients and even separate age populations within these associations.Additionally, effects such as extinction, accretion history, and binarity affect the observational uncertainty when determining ages (Krolikowski et al. 2021).These can lead to vertical scatter on Hertzsprung-Russell diagrams (HRD; Baraffe et al. 2017), a result of the variation in the radius-to-mass ratio driven primarily by variations such as age (Pecaut et al. 2012;Soderblom et al. 2014;Rizzuto et al. 2016).These large uncertainties make individual age estimates hard to derive from HRDs and rife with large uncertainties (Pecaut et al. 2012;Malo et al. 2014;Feiden 2016;Rizzuto et al. 2020).Additionally, other methods of calculating ages (i.e.lithium burning, kinematics) all have intrinsic assumptions and systematic effects on the error (see Soderblom et al. 2014, for a detailed review of these effects).
Therefore, recent studies of SFRs have relied on robustly estimating ages of groups within SFRs, which reduces the vertical spread resulting from assuming a single age as well as the large uncertainty from assuming individual ages (e.g.Galli et al. 2020Galli et al. , 2021;;Krolikowski et al. 2021;Esplin & Luhman 2020).This method relies on placing objects on a HRD and comparing to known isochrone and evolutionary track models to estimate age and mass, respectively.Though there are systematic and observational uncertainties due to potential differences in models used, and assumptions of temperature and luminosity, this method has been employed previously for all of the associations found in CASPAR, providing a (semi-)uniform method to derive ages.We investigate if assuming a single age for the whole SFR affects the resulting masses and accretion rates (and the effect on the Ṁ − M scatter).Using Taurus, as there is a large number of subregions and ages, we compare the radius-to-mass ratio derived from assuming a single age (2 ± 1 Myr) and ages from the separate regions (see Table 7).We find that assuming a single age for a whole SFR only affects the radius/mass ratio and accretion rate by less than < 20% compared to using region ages.This is generally less than the uncertainty on both individual and SF ages (∼ 50%) and likely does not greatly affect the derived accretion rates.We therefore use the average age from the most recent census papers for each region (see Table 7) corresponding to the (sub)cluster/association that the object is a member.Individual object ages estimated from isochrone fitting from the literature reference are also included in the database.

B.2. Literature Database
For each accretion rate, we compile the literature reference, and physical and accretion properties used to calculate Ṁ .This includes the association, age, and distance specific to the literature reference, as well as the spectral type, mass, luminosity and radius.If any quantity is not specified in its reference, it is left blank.We also record, if given, the reference used to convert from spectral type to temperature, the evolutionary model used to derive estimates of physical parameters, and the emission line luminosity to accretion luminosity (L line −L acc ) scaling relation.For each accretion rate, we list the specific accretion signature or continuum band used to calculate Ṁ .We list the main accretion diagnostic for each reference in Table 8.We also compile all individual line emission quantities (EWs, line fluxes, accretion luminosities, and accretion rates) reported in the paper.In cases where the accretion rate is not presented in the reference, but can be calculated, we have done so and included them in the database.Finally, the accretion luminosity and rate found from the main reported literature reference are recorded.Objects with multiple measurements (from multiple studies, tracers, or epochs) are connected with a gray dashed line, and highlight the variations in published estimates of mass  2021), who estimate the 1σ span of the highest density interval on the posterior probability density used to determine the distance.This 1σ span sets the lower and upper bounds of the distance and is not assumed to be symmetric around the median distance.
Age-We assume the age and uncertainty of the star forming region or association as described in Section B.1.We refer to the individual papers for full details on age and uncertainty determination.
Spectral Type-For objects that have only one measured accretion rate in the literature database (N = 646), we assume the literature reference spectral type.For objects with multiple measured accretion rates (N = 152), we assume the most recently measured spectral type.For objects with no spectral types listed in their reference paper, we searched VizieR (Ochsenbein et al. 2000) and SIMBAD (Wenger et al. 2000) for a spectral type.If none was found, but a temperature was given, we either a) calculated the spectral type using the spectral type to temperature conversion if stated in the literature reference, or b) calculated the spectral type using the spectral type to temperature conversions of Herczeg & Hillenbrand (2014) if there was none stated.For spectral types without listed uncertainties, we adopt the spectral class uncertainty of Herczeg & Hillenbrand (2014): M dwarfs: 0.2 subclass, K8-M0.5:0.5 subclass, G0-K8: 1 subclass.
Effective Temperature-Temperature is calculated using the spectral type to temperature conversion of Herczeg & Hillenbrand (2014).The uncertainty is found by calculating the temperature at the upper and lower limits of the spectral type estimate.We then take the average value of the difference between the bounds and the given temperature as the uncertainty.
Mass, Luminosity, Radius, Surface Gravity-To consistently estimate mass, luminosity and radius, we use the evolutionary models of Baraffe et al. (2015) and the MIST MESA models (Dotter 2016;Choi et al. 2016;Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015)).Using object age and temperature, we interpolate over the isochronal models to determine the mass, luminosity, radius, and surface gravity.We heavily modified the isochrone6 python package (Morton 2015) to work with the Baraffe et al. (2015) models and to interpolate between temperatures (the package currently only interpolates between masses).After interpolating over age and effective temperature, the best fit mass, luminosity, radius, and surface gravity can be extracted.To determine the uncertainty on these quantities, we find the lower and upper limits on the age and temperature; this produces four bounds (1-low age/low temperature, 2-low age/high temperature, 3-high age/low temperature, 4high age/high temperature).We then take the average of the difference between the bound values and the given value as the uncertainty.We assume that all reported spectral types and line luminosities have been corrected for extinction.

C.2. Accretion Parameters
Below we detail the unified methodology used to derive accretion rates for various accretion diagnostics.
Continuum Excess-Accretion Luminosities are determined from the total excess luminosity (derived from spectral template fitting), which has traditionally been assumed to result primarily from the continuum excess.From Herczeg & Hillenbrand (2008), L acc ∝ d 2 .Therefore, using the original literature accretion luminosities, we scale them by the Gaia distances, and then derive an accretion rate using the updated accretion luminosity, mass, and radius.
Line Luminosities-For papers that report line fluxes (excluding the Mohanty et al. (2005) Ca ii λ 8662 line, see below), we calculate mass accretion rates from the reported line flux under a single set of scaling relations.For references that only report the accretion luminosity or accretion rate of the line(s), we use the L line − L acc scaling relations and distance given in the paper to infer the measured line flux.For studies that only give accretion rates, we first calculate L acc using the reference mass and radius and then proceed as above.Once we have derived a line flux, we calculate the accretion rate as follows: 1.We update the line luminosity using the Gaia DR2 or EDR3 distance and assuming isotropic emission: L line = 4πd 2 F line .2. We convert line luminosity to accretion luminosity following log(L acc /L ⊙ ) = a × log(L line /L ⊙ ) + b) using a single set of independently derived scaling relations (those of Alcalá et al. 2017;Salyk et al. 2013).3. We convert to accretion rate using our rederived mass and radius with Equation (1).For objects with given line flux uncertainties, we propagate this error forward using the uncertainties on Gaia distance, rederived mass and radius uncertainties, and the uncertainties on the scaling relationship.
The Ca ii 8662 Å line fluxes calculated by Mohanty et al. (2005) are a unique case, as the line fluxes are calculated assuming a best fit modeled continuum, causing systematic offsets between their values and "true" values (for full details see Mohanty et al. 2005).Therefore, we use their scaling relations (their eqns 1. and 3.) in step 2 above.
Hα Photometric Luminosities-Narrowband Hα photometry has frequently been used to calculate mass accretion rates for substellar objects (e.g., Kalari et al. 2015;Sallum et al. 2015;Wu et al. 2015Wu et al. , 2017;;Close et al. 2014;Wagner et al. 2018).Once the Hα luminosity is determined, an L Hα − L acc scaling relation can be applied.Therefore, in order to recalculate the accretion rate, we follow the same procedure as our rederivation of spectroscopic line fluxes (see above), substituting recalculated physical parameters (distance, mass, radius) and adopting a single scaling relation (Alcalá et al. 2017).
However, for 16 objects in the Kalari et al. (2015) Lagoon Nebula sample (D ∼ 1326 pc), Gaia EDR3 distances put them at ∼ 200 − 700pc.Henderson & Stassun (2012) proposed a that the cluster might be 15% closer than previous estimates; however, even with a closer distance, these objects are still well below any assumed distance to the cluster.Additionally, ten of the objects have proper motions with Right Ascension and Declination dispersion greater than best-fit values found for the Lagoon Nebula (σ µα ∼ 4.06 km/s, σ µ δ ∼ 2.8 km/s, Wright et al. 2019) assuming a distance of 1326 pc.Therefore, we assume these objects are not true members of the Lagoon Nebula, and their true ages are unknown.Therefore, their masses cannot be estimated properly and we exclude them from further analysis.
Line Profiles-Accretion rates found by modeling the Hα emission profile with radiative transfer models of the magnetosphere rely on the line-of-sight inclination angle of the disk and velocity field.Gas velocities are sensitive to mass (Muzerolle et al. 2001).Though the input mass range will vary according to the evolutionary tracks used, the best fit model diverges if the mass is varied by a factor of two or more (due to significant variation in the model gas velocity; v gas ∝ √ M ).For object radii, uncertainty results from the spectral type-temperature conversion, with a 200 K error (equivalent to 1 spectral subclass) equal to 30% error in radius.Nonetheless, as the modeled gas density is not strongly dependent on radius (is depends on the system geometry), the model is less sensitive to variation.Therefore, as long as the rederived masses are approximately within factor of two, temperatures within 200 K (1 spectral subclass), and/or radii within 30% of the literature values, no rederivation is needed.From Fig. 19, we find that the average difference in temperature is 94 K (corresponding to half a spectral class) and the masses are within a factor of 0.83 and therefore the accretion rates do not have to be recalculated.
Accretion rates have been found to scale directly with Hα 10% width (Natta et al. 2004).Therefore, accretion rates originally found using the Natta et al. (2004) scaling relation do not have to be recalculated; this encompasses all objects with Hα 10% width within CASPAR.

D. OTHER LINEAR FITTING TECHNIQUES
It is well known that different linear regression procedures should be used based on the data being considered.Significant work (Isobe et al. 1990;Feigelson & Babu 1992) has explored how these different algorithms affect astronomical data and results.As we explore the best linear fits for a variety of quantities in CASPAR, reliable fitting is necessary.linmix, a Bayesian linear regression routine (Kelly 2007), takes into account both x and y measurement errors as well as upper limits.In order to determine the extent to which these quantities affect the linear fits, we compare this fit to other fitting techniques that take into account a variety of these parameters.See Table 9.We use python to derive the best fits for each technique:   10 show the results.Algorithms which include x errors but do not properly take into account upper limits (ODR, BCES) still closely match linmix, while those that do not take into account errorbars do not reproduce linmix.This is most prominent in the substellar regime, where the fits diverge significantly.Those without x errors consistently overestimate Figure 19.Absolute difference between literature and CAS-PAR derived temperatures/spectral types as a function of CASPAR temperature/spectral type (top) and the ratio between literature and CASPAR derived masses as a function of CASPAR mass (bottom) for objects with accretion rates derived from Hα emission profile modeling.Estimated accretion rates do not vary with rederived masses and radii if the spectral types are within 1 subclass (< 200 K difference) and masses are within a factor of two of the original values.We find all values in CASPAR to be within those parameters.
the intercept of the fit compared to those with take both x and y into account.

E. SEPARATE MASS POPULATION FIT STATISTICS
In Table 11, we compare the mean residual from each fit, the AIC, and Akaike weights (w) for each fit.As our models have the same number of parameters, the AIC will inform the goodness of fit between the two models, with the minimum AIC corresponding to the preferred   Brown Dwarfs: In the substellar regime, we find a small, though not statistically significant, difference in the residuals between the CASPAR total fit and the brown dwarf-only fit, with the mean decreasing from 0.25 to 0.23.However, from the AIC statistics and weights shown in Table 11, we can show that the BDonly fit is more significantly preferred than the CASPAR total best fit.
PMCs: The best fit clearly showed non-normal residuals for the PMCs.The PMC-only best fit is significantly preferred with the mean of the residuals decreasing from 0.82 to 0.05.

F. SEPARATE MASS POPULATION FITS FOR INDIVIDUAL STAR FORMING REGIONS
We examine substellar and stellar population best fits for individual star forming regions in CASPAR, removing the need to fit for age.We select the ρ Ophiuchus (∼ 2 Myr), Taurus (∼ 2.5 Myr), and Upper Scorpius (∼ 10 Myr) regions as they span a wide range of ages.As shown in Fig. 21, from just these three regions, we see a steepening slope (α = 0.52, 1.61, 4.66) with age in the substellar regime, while the stellar regime appears to slightly flatten (α = 1.32, 2.53, 0.94) for ρ Ophiuchus, Taurus, and Upper Scorpius, respectively, similar to the trends seen by Alcalá et al. (2017) and Manara et al. (2017a).

Figure 1 .
Figure1.All-sky map indicating all objects in CASPAR colored by star-forming region or association overlaid on a 65 µm all-sky map(Doi et al. 2015;Takita et al. 2015).Insets: Enlarged views of several nearby star forming regions.

Figure 5 .
Figure 5.Comparison between the Literature Database physical parameters and those rederived in CASPAR colored by accretion diagnostic.The black line indicates 1-1.For masses, we show the hydrogen burning limit at M = 0.075 M⊙ with dashed black lines; under a uniform methodology, objects in the lower right quadrant decreased in mass from stellar to substellar, while in the upper left, the reverse occurred.We find 50−89% of CASPAR parameters are within a factor of two of the literature value.

Figure 6 .
Figure 6.Residuals between reference and CASPAR accretion rates as a function of CASPAR accretion rate for different accretion diagnostics.88% of all rederived values changed by less than 0.5 dex, with the line luminosity accretion diagnostic showing the most change in calculated accretion rate.right: Histogram of difference between literature and rederived accretion rates for each tracer.The thin gray lines indicate the FWHM.Colors as in Fig. 5.

Figure 7 .
Figure 7.The CASPAR Ṁ −M relation for stars (black circles), brown dwarfs (cyan diamonds), and PMCs (magenta squares).The gray dashed lines show accretion rates derived for the same object, and the downward arrows show accretion rate upper limits.The black line and shaded region shows our best linear fit, log Ṁ = 2.02 log M − 8.02, and 1σ dispersion, σ = 0.85 dex ( Ṁ in M⊙ yr −1 and M in M⊙), to all accretion rates in CASPAR.

Figure 8 .
Figure 8. left: The CASPAR Ṁ − M best-fit relation for overall CASPAR (black solid line), CASPAR stars-only (black dashed dot line), Zhou et al. (2014) (gray dotted line), and Hartmann et al. (2016) (brown solid line).The 1σ scatter is shown as shaded regions for each fit.The residuals between each linear fit and CASPAR are shown below.right: Histogram of residuals for each population.The thin gray lines shows the mean of each distribution.Markers and colors are as in Fig. 7.

Figure 9 .
Figure 9. Panel a. Accretion rate vs mass for individual accretion diagnostics colored as in Fig. 4. The large markers indicated the binned accretion rates by mass.Panel b.Residuals between binned accretion rates and CASPAR continuum excess fit.Line luminosities (including Hα photometric) show (significant) deviation from the continuum excess with decreasing mass.Hα photometric luminosities appears significantly offset from the CASPAR continuum excess best fit; however, the majority of the stellar accretion rates are from objects with ages < 1 Myr, which biases the high mass results (see Sec 6.1).

Figure 10 .
Figure10.Accretion rate vs mass for a range of accretion-tracing lines, namely: a) optical Hydrogen Balmer series lines b) NIR Hydrogen Paschen, Brackett and Pfund series lines, c) Helium i emission lines and d) Calcium ii emission lines.The gray line shows the best linear Ṁ − M fit to the overall CASPAR database, while the colored lines show the best linear fit for each emission line.The histograms show the posterior distributions for the masses (top), accretion rates (top right), and accretion rate residuals (bottom right) for each panel.We find that the NIR line flux fit shows the most deviation in best fit; it is offset by 0.9 dex compared to the best linear fit at the deuterium burning limit (compared to < 0.4 in the other lines).
with a scatter of 0.77 dex and Pearson correlation coefficient of −0.31.The slope found by Hartmann et al. (2016) (slope = −1.07)shows a faster decline in accretion rate with age compared to the CASPAR slope.The sample of Hartmann et al. (

Figure 12 .
Figure 12.CASPAR accretion rate vs age colored by mass.The accretion rates have been scaled by the Ṁ − M relation to M ∼ 0.7 M⊙ following Hartmann et al. (2016).The black line and shaded region show the best linear fit and 1σ scatter.

Figure 13
Figure 13.top: CASPAR accretion rate vs mass colored by ages: ≤1 Myr (gold circles), (1-3] Myr (orange stars), (3-8] Myr (magenta squares), >8 Myr (purple triangles).The solid lines are the best linear fits to each age population with a fixed constant slope of 2.02.We find a 1 dex decrease in the best-fit intercept with increasing age.bottom: Residuals for the best linear fit for each age population.

Figure 14 .
Figure14.Range of measured accretion rates for objects in CASPAR with a) multi-epoch/multi accretion tracer measurements, b) all accretion rates measured for all emission lines in one epoch, and c) all accretion rates in all epochs for individual emission lines.In panels a) and b), the ranges are colored by stellar/substellar populations.In panel c), the ranges are colored by the individual line tracers.The right panels shows the probability density functions for each population.Overall, we find increased variability in the substellar regime (1.33 dex), but consistent 0.8 dex variability with population for emission lines and epochs.

Figure 15 .
Figure15.CASPAR Ṁ − M best fit relations for stars (black circles), brown dwarfs (cyan diamonds), and planetary mass companions (magenta squares).The solid lines and bands show the best linear fits and 1σ standard deviation of the fit for each mass population, while the dashed lines show extrapolations beyond the bounds of the fit regions.The bottom panel shows the residuals for each best fit.Downward arrows indicate accretion rate upper limits.

Figure 17 .
Figure17.Best fit log Ṁ as a function of age at the hydrogen burning limit (M = 0.075 M⊙; solid line and filled markers) and at the deuterium burning limit (M = 0.012 M⊙; dashed line and open markers) for stars (black circles) and brown dwarfs (cyan diamonds) from Fig.16.We show the best exponential fits for each population and find the brown dwarf accretion rates decrease at a faster rate compared to the stellar.

Figure 18 .
Figure18.Accretion rate vs mass for the star and brown dwarf population colored by specific ages.The solid lines and bands show the best linear fits and 1σ standard deviation where age and mass are both free parameters for each mass population.The dashed lines are extrapolations to the DBL.

Fig
Fig.20and Table10show the results.Algorithms which include x errors but do not properly take into account upper limits (ODR, BCES) still closely match linmix, while those that do not take into account errorbars do not reproduce linmix.This is most prominent in the substellar regime, where the fits diverge significantly.Those without x errors consistently overestimate

Figure 21 .
Figure 21.Accretion rate vs mass for ρ Ophiuchus (left), Taurus (center), and Upper Scorpius (right).Colors and markers are as in Fig. 15.The best linear fits (solid line) and 1σ scatter (band) are shown for the star (black) and brown dwarf (cyan) populations.The dashed lines show the extrapolation into the planetary regime.

Table 1 .
Literature Reference Star Forming Regions

Table 2 .
CASPAR sections * * See Table 6 for description of individual columns

Table 4 .
Best fit parameters

Table 6 .
Contents of CASPAR and Literature Database

Table 6 continued
Table 6 (continued) He i λ 4471, He i λ 4713, He i λ 5016, He i λ 5876, He i λ 6678, He i λ 7065, He i λ 10830, He ii λ 4686, Ca ii K, Ca ii H, Ca ii λ 8498, Ca ii λ 8542, Ca ii λ 8662, Na i λ 5889, Na i λ 5896, O i λ 8446, C iv λ 1549 Note-Only a portion of this table is shown here to demonstrate its form and content.Machine readable versions of CASPAR (Part 1) and the Literature Database (Part 2) are available.

Table 7 .
Star forming regions and Association ages and distances Jones et al. 2018).If the object had no Gaia measurement, we assume the average distance to the region listed in Table 7 (N = 189 in Appendix B).The uncertainty on each Gaia distance is from either Bailer-Jones et al. (2018) or Bailer-Jones et al. (

Table 8 .
Literature References

Table 9 .
Linear regression algorithms fitting algorithm y err x err upper limits CASPAR Accretion rate vs mass with linear best fits from different linear regression algorithms.linmixtakesintoaccountall parameters, while the others take into account some parameters as listed in Table9.model.The Akaike weights are the relative likelihoods of the models; we assume that w > 0.95 indicates the statistically favored model.Stars: We find similar AIC statistics for the best total fit to CASPAR and star-only fit, with neither model preferentially preferred.This indicates that either fit can be used to model the data.

Table 11 .
Population statistics 2 : coefficient of determination