PROTOPLANETARY DISK MASSES FROM STARS TO BROWN DWARFS

The masses of the primordial disks girdling newborn stars and brown dwarfs (BDs), and the degree of grain growth in these disks, are key to the processes of accretion, planet formation, and migration within them. Much work has gone into inferring disk masses and grain growth around solar-type and higher-mass stars (e.g., Beckwith et al. 1990; Andrews & Williams 2005, 2007a; Ricci et al. 2010a, 2010b, hereafter B90, AW05, AW07a, R10a,b, respectively), and such studies now extend into the substellar domain as well (Klein et al. 2003; Scholz et al. 2006; Schaefer et al. 2009).

Here, we present the results of a JCMT/SCUBA-2 850 μm pilot survey of seven very low mass stars (VLMS) and BDs, spanning 0.02–0.2 M_☉ in mass, and located in the ∼1 Myr old Taurus star-forming region and the ∼10 Myr old TW Hydrae Association (TWA). All the sources are known to be accreting from surrounding primordial disks. The study, undertaken as part of SCUBA-2 first-light observations, is 3–10 times deeper than previous such surveys of young VLMS/BDs, and includes four out of the five confirmed VLMS/BD accretors in the TWA, none of which have been examined at these wavelengths before.⁹

Our goals are twofold. First, by combining our data with previous large submillimeter/millimeter (sub-mm/mm) disk surveys (specifically, those of AW05; AW07a; Scholz et al. 2006; Schaefer et al. 2009), we wish to investigate trends in disk mass and grain growth as a function of stellar mass, and the attendant implications for planet formation and accretion. Our analysis generally follows that of AW05 and AW07a, with a few significant differences: (1) we concentrate solely on Class II objects (or, roughly equivalently, classical T Tauri (cTT) sources), in order to avoid confusion by both envelope contamination in earlier evolutionary types and a possible lack of primordial disks altogether in more evolved sources; (2) we combine the Taurus and ρ Oph samples of AW05 and AW07a, and extend the study to much lower (sub)stellar masses by the addition of objects from our own and other surveys; (3) we take advantage of various more recent investigations of disk radii, surface densities, and temperatures, in order to define a set of realistic fiducial disk parameters that allows us to focus on the primary unknowns of interest here—disk mass and grain growth; and (4) we explicitly frame the analysis in terms of generalized equations that extend the Rayleigh–Jeans (RJ) formalism of B90 to the non-RJ regime (which we show is especially vital for VLMS/BDs).

Second, and equally important, we wish to present a Bayesian framework for some of the above analysis, which makes maximal use of the non-detections/upper limits in the combined data set. The majority of disks around VLMS and BDs, as well as a significant fraction of those around solar-type and higher-mass stars, remain undetected in the sub-mm/mm. While these non-detections clearly have something to say about the underlying disk mass distribution, the question is how to combine them with the detections in order to extract the maximum information encoded in all the data. This is obviously not an issue restricted to studies of disk masses, but one that is central to all surveys that include upper limits. Unfortunately, it has often not been addressed satisfactorily in the young stellar community. Non-detections are frequently simply ignored, or upper limit values—ranging arbitrarily from 2σ to 5σ—used as true detections in order to estimate the sample distribution. A more sophisticated approach has sometimes been to invoke survival analysis, based on the classic work by Feigelson & Nelson (1985). However, as the latter authors explicitly caution, this technique (1) is not appropriate when the upper limits are correlated with the variable under discussion (e.g., if the upper limits on the observed flux are primarily due to the sensitivity of a flux-limited survey), which is clearly the case in many if not most astronomical studies; and (2) does not account for noise in the data. A Bayesian approach, on the other hand, provides an elegant and intuitively simple way of dealing with both upper limits and noise. We outline the general method, and apply it to our combined disk mass data set; we hope that the community adopts the technique more widely for analyzing analogous surveys.

Our sample, observations, and data reduction procedure are summarized in Section 2, and our derivation of stellar parameters is described in Section 3. Section 4 provides a brief overview of the equations we use to model the disk spectral energy distribution (SED), with more details in Appendix A. We summarize our method for analyzing grain growth and disk mass in Section 5, and discuss our choice of fiducial disk parameters for this analysis in Appendix B. Our results are described in Sections 6–8, which cover the validity of the RJ approximation (Section 6), trends in grain growth and disk mass (Section 7), and a Bayesian analysis of the ratio of disk to stellar mass (Section 8, with an outline of the general Bayesian technique in Appendix C). The implications for planet formation and mass accretion rates are discussed in Section 9, and our conclusions presented in Section 10.

2.1. SCUBA-2 Sample and Additional Sources

We were awarded 10 hr of shared-risk (first light) time on JCMT/SCUBA-2 to investigate primordial disks around VLMS/BDs. The observed sample consists of seven young objects: CHFT-BD Tau 12 (M6.5), GM Tau (M6.5), and J044427+2512 (M7.25) in Taurus, and TWA 30A (M5), TWA 30B (M4), 2MASS 1207-3932 (M8), and SSSPM 1102-3431 (M8) in the TWA. All are optically revealed sources (i.e., without surrounding remnant envelopes) known to host accretion disks (from optical/UV accretion signatures and infrared dust emission), i.e., they are Class II cTT analogs. This choice allows us to focus on disk properties, without confusion from either envelope emission (Class 0/I sources) or the absence of a disk altogether (Class III). J044427+2512 had been detected earlier over 450 μm–3.7 mm (Scholz et al. 2006; Bouy et al. 2008), and was selected to test our sensitivity; GM Tau had been observed before at 1.3 mm and 2.6 mm but not detected (Schaefer et al. 2009); and the rest had never been observed earlier in the sub-mm/mm. Our four TWA sources, combined with the previously observed Hen 3-600A and TW Hya (see below), comprised all known VLMS/BD accretors, and six out of the seven known accretors of any stellar mass, in this association at the time of observation. The one confirmed TWA accretor excluded here is TWA 5A (Mohanty et al. 2003), which is a roughly solar-type cTT multiple not yet observed in the sub-mm/mm.¹⁰ Very recently, Schneider et al. (2012) have announced two new VLMS TWA members, TWA 33 and 34. Their relatively weak Hα emission suggests very little accretion (applying the Hα equivalent width criterion devised by Barrado y Navascues & Martin); nevertheless, the mid-infrared excesses in WISE bands indicate surviving primordial disks (Schneider et al. 2012). These two sources have not yet been observed in the sub-mm/mm, and were announced too late for inclusion in our survey.

To our sample observed with SCUBA-2, we add most of the Taurus, ρ Oph, and TWA Class II and/or cTT analogs (i.e., optically revealed accretors) observed in the sub-mm/mm (850 μm and/or 1.3 mm) so far. The Taurus and ρ Oph data are taken from the surveys by AW05, AW07a, Scholz et al. (2006), and Schaefer et al. (2009). From the first two, we select only those classified as Class II (based on power-law fits over 2–60 or 2–25 μm; see AW05, AW07a). From Scholz et al. (2006), we only choose the nine objects classified as cTTs by Mohanty et al. (2005), Muzerolle et al. (2005), or Luhman (2004) based on optical/infrared diagnostics: CFHT-BD Tau 4, J041411+2811, J043814+2611, J043903+2544, J044148+2534, J044427+2512, KPNO Tau 6, KPNO Tau 7, and KPNO Tau 12. These are all in Taurus, and many of them are also known to harbor disks from Spitzer mid-infrared data (Luhman et al. 2010). There is a potential danger that such a selection criterion might discard very weak accretors which nonetheless still harbor disks (which describes a number of Class II sources). In practice, however, we find that this criterion includes all the sources Scholz et al. detect at 1.3 mm, as well as a few which they do not. In other words, the sources we have dropped present no evidence at all of disks, from optical to mm wavelengths, and there is no rationale for including them at this point: while a few may very well harbor disks that are currently completely undetected, the same may be said of objects classified as Class III by AW05 and AW07a, which we have also ignored.

From the Taurus survey by Schaefer et al. (2009), we include all sources with spectral-type M4 or later (i.e., VLMS/BDs) that have been classified as cTTs by Mohanty et al. (2005) or Muzerolle et al. (2003), and/or as Class II by Luhman et al. (2010; based on mid-infrared Spitzer data out to 24 μm). This yields seven objects: CIDA 1, CIDA 14, FN Tau, FP Tau, GM Tau, MHO 5, and V410 Anon 13 (an eighth object fitting these criteria, CIDA 12, has been previously observed by AW05, and we use their value for consistency). We ignore the remaining 13 stars (all earlier than M4) from Schaefer et al.'s survey, for the following reasons. Seven of these have also been observed by AW05, of which six are classified by them as Class II (GO Tau, DN Tau, IQ Tau, CIDA 7, CIDA 8, and CIDA 11); these are already included in our study (with AW05's values). Among the remaining six stars not observed by AW05, various issues arise in some (e.g., large spectral-type uncertainty, or not classified as cTTs or Class II), so we conservatively choose to discard all six. Given the large sample at these spectral types that our study already includes, from AW05 and AW07a, this decision has no discernible impact on our results.

Finally, we include the sub-mm/mm fluxes for the TWA cTTs Hen 3-600A and TW Hya from Zuckerman (2001) and Weintraub et al. (1989), respectively.

Our combined sample consists of 134 objects—48 in ρ Oph, 80 in Taurus, and 6 in the TWA—spanning masses of ∼15 M_J to 4 M_☉. For uniformity, we recalculate the stellar and disk parameters for all our sources taken from the literature, based on the spectral types and disk fluxes cited by the authors and using the methods discussed in Sections 3 and 5. The full final sample is listed in Table 1, with our SCUBA-2 sub-sample marked with asterisks.

Table 1. Sub-mm/mm Fluxes, Binarity, and Derived Properties for Class II/cTTs in ρ Oph, Taurus, and the TWA

IDa	SpT	Teff	M*	F[850]b	1σ Errorb	F[1300]b	1σ Errorb	α850 − 1.3c	Md, νd	1σ Error	Region	Dist	Refe	Multf	Refg
		(K)	(M☉)	(mJy)	(mJy)	(mJy)	(mJy)		(MJup)	(MJup)		(pc)
AS 205	K5	4350	1.000	891	11	450	10	1.61 ± 0.06	59.0070	0.7285	ρ Oph	150.	1, 1	w, sb	1
AS 209	K5	4350	1.000	551	10	300	10	1.43 ± 0.09	36.4903	0.6623	ρ Oph	150.	1, 1	...	...
DoAr 16	K6	4205	0.850	47	8	<50	...	>−0.15	3.1126	0.5298	ρ Oph	150.	1, 1	...	...
DoAr 24	K5	4350	1.000	...	...	<30	...	...	<6.0340	...	ρ Oph	150.	..., 1	...	...
DoAr 24E	K1	5080	2.481	158	6	70	20	1.92 ± 0.68	10.4636	0.3974	ρ Oph	150.	1, 1	w	1
DoAr 25	K5	4350	1.000	461	11	280	10	1.17 ± 0.10	30.5300	0.7285	ρ Oph	150.	1, 1	...	...
DoAr 28	K5	4350	1.000	...	...	<75	...	...	<15.0851	...	ρ Oph	150.	..., 1	...	...
DoAr 32	K6	4205	0.850	...	...	<45	...	...	<9.0511	...	ρ Oph	150.	..., 1	...	...
DoAr 33	K4	4590	1.295	79	7	40	10	1.60 ± 0.62	5.2318	0.4636	ρ Oph	150.	1, 1	...	...
DoAr 44	K3	4730	1.538	181	6	105	11	1.28 ± 0.26	11.9868	0.3974	ρ Oph	150.	1, 1	...	...
DoAr 52	M2	3560	0.579	...	...	<55	...	...	<11.0624	...	ρ Oph	150.	..., 1	...	...
EL 18	K6	4205	0.850	...	...	<10	...	...	<2.0113	...	ρ Oph	150.	..., 1	...	...
EL 24	K6	4205	0.850	838	8	390	10	1.80 ± 0.06	55.4970	0.5298	ρ Oph	150.	1, 1	...	...
EL 26	M0	3850	0.691	<51	...	15	5	<2.88	3.0170	1.0057	ρ Oph	150.	1, 1	...	...
EL 27	K8	3955	0.740	678	10	300	10	1.92 ± 0.09	44.9009	0.6623	ρ Oph	150.	1, 1	...	...
EL 31	M0	3850	0.691	...	...	<10	...	...	<2.0113	...	ρ Oph	150.	..., 1	...	...
EL 32	K7	4060	0.785	...	...	<50	...	...	<10.0567	...	ρ Oph	150.	..., 1	...	...
EL 36	A7	7850	3.639	...	...	<10	...	...	<2.0113	...	ρ Oph	150.	..., 1	...	...
GSS 26	K8	3955	0.740	298	7	125	20	2.04 ± 0.38	19.7352	0.4636	ρ Oph	150.	1, 1	...	...
GY 284	M3	3415	0.478	...	...	130	10	...	26.1475	2.0113	ρ Oph	150.	..., 1	...	...
IRS 2	K3	4730	1.538	...	...	<25	...	...	<5.0284	...	ρ Oph	150.	..., 1	...	...
IRS 37	M4	3270	0.286	93	8	<10	...	<5.25	6.1590	0.5298	ρ Oph	150.	1, 1	...	...
IRS 39	M2	3560	0.579	63	5	<15	...	>3.38	4.1722	0.3311	ρ Oph	150.	1, 1	...	...
IRS 49	K8	3955	0.740	52	5	25	5	1.72 ± 0.52	3.4437	0.3311	ρ Oph	150.	1, 1	...	...
RNO 90	G5	5770	3.220	162	4	25	5	4.40 ± 0.47	10.7285	0.2649	ρ Oph	150.	1, 1	...	...
ROXs 25	M2	3560	0.579	...	...	30	5	...	6.0340	1.0057	ρ Oph	150.	..., 1	...	...
ROXs 42C	K6	4205	0.850	...	...	<30	...	...	<6.0340	...	ρ Oph	150.	..., 1	...	...
ROXs 43A	G0	6030	3.354	...	...	<35	...	...	<7.0397	...	ρ Oph	150.	..., 1	...	...
SR 4	K5	4350	1.000	142	7	31	6	3.58 ± 0.47	9.4040	0.4636	ρ Oph	150.	1, 1	...	...
SR 9	K5	4350	1.000	<25	...	15	5	<1.20	3.0170	1.0057	ρ Oph	150.	1, 1	...	...
SR 10	M2	3560	0.579	...	...	<25	...	...	<5.0284	...	ρ Oph	150.	..., 1	...	...
SR 13	M4	3270	0.286	118	6	60	10	1.59 ± 0.41	7.8146	0.3974	ρ Oph	150.	1, 1	c	1
SR 21	G3	5830	3.251	397	6	95	15	3.37 ± 0.37	26.2916	0.3974	ρ Oph	150.	1, 1	w	2
SR 22	M4	3270	0.286	31	3	<20	...	1.03	2.0530	0.1987	ρ Oph	150.	1, 1	...	...
Wa-Oph 4	K4	4590	1.295	...	...	<13	...	...	<2.6148	...	ρ Oph	150.	..., 1	...	...
Wa-Oph 5	M2	3560	0.579	...	...	<25	...	...	<5.0284	...	ρ Oph	150.	..., 1	...	...
Wa-Oph 6	K6	4205	0.850	379	7	130	10	2.52 ± 0.19	25.0995	0.4636	ρ Oph	150.	1, 1	...	...
WL 10	K8	3955	0.740	...	...	<60	...	...	<12.0681	...	ρ Oph	150.	..., 1	...	...
WL 14	M4	3270	0.286	...	...	30	10	...	6.0340	2.0113	ρ Oph	150.	..., 1	...	...
WL 18	K7	4060	0.785	...	...	85	10	...	17.0965	2.0113	ρ Oph	150.	..., 1	w	1
WSB 19	M3	3415	0.478	<78	...	...	...	...	<5.1656	...	ρ Oph	150.	1, 1	...	...
WSB 37	M5	3125	0.150	...	...	<25	...	...	<5.0284	...	ρ Oph	150.	..., 1	...	...
WSB 46	M2	3560	0.579	...	...	<20	...	...	<4.0227	...	ρ Oph	150.	..., 1	...	...
WSB 49	M4	3270	0.286	...	...	<25	...	...	<5.0284	...	ρ Oph	150.	..., 1	...	...
WSB 52	K5	4350	1.000	...	...	51	10	...	10.2579	2.0113	ρ Oph	150.	..., 1	...	...
WSB 60	M4	3270	0.286	149	7	89	2	1.21 ± 0.12	9.8676	0.4636	ρ Oph	150.	1, 1	...	...
WSB 63	M2	3560	0.579	...	...	<25	...	...	<5.0284	...	ρ Oph	150.	..., 1	...	...
YLW 16C	M1	3705	0.596	...	...	60	5	...	12.0681	1.0057	ρ Oph	150.	..., 1	...	...
04113+2758	M2	3560	0.579	...	...	410	40.	...	71.8364	7.0084	Taurus	140.	..., 2	w	3
04278+2253	F1	7050	3.567	36.	7.	...	...	...	2.0768	0.4038	Taurus	140.	2, ...	...	...
AA Tau	K7	4060	0.785	144.	5.	88	9.	1.16 ± 0.25	8.3073	0.2884	Taurus	140.	2, 2	...	...
AB Aur	A0	9520	3.791	359.	67.	103	18.	2.94 ± 0.60	20.7107	3.8652	Taurus	140.	2, 2	...	...
BP Tau	K7	4060	0.785	130.	7.	47	0.7	2.39 ± 0.13	7.4997	0.4038	Taurus	140.	2, 2	...	...
CFHT-BDTau 4	M7	2880	0.057	10.8	1.8	2.38	0.75	3.56 ± 0.84	0.6230	0.1038	Taurus	140.	3, 4	...	...
*CFHT-BDTau 12	M6.5	2935	0.074	4.06	1.32	...	...	...	0.2342	0.0762	Taurus	140.	5, ...	...	...
CIDA 1	M5.5	3058	0.125	...	...	13.5	2.8	...	2.3653	0.4906	Taurus	140.	..., 6	...	...
CIDA 3	M2	3560	0.579	<9.	...	...	...	...	<0.5192	...	Taurus	140.	2, ...	...	...
CIDA 8	M4	3270	0.286	27.	3.	...	...	...	1.5576	0.1731	Taurus	140.	2, ...	...	...
CIDA 9	M0	3850	0.691	71.	7.	...	...	...	4.0960	0.4038	Taurus	140.	2, ...	w	4
CIDA 11	M3	3415	0.478	<8.	...	...	...	...	<0.4615	...	Taurus	140.	2, ...	c	4
CIDA 12	M4	3270	0.286	<7.	...	...	...	...	<0.4038	...	Taurus	140.	2, ...	...	...
CIDA 14	M5	3125	0.150	...	...	<4.5	...	...	<0.7884	...	Taurus	140.	..., 6	...	...
CI Tau	K7	4060	0.785	324.	6.	190	17.	1.26 ± 0.22	18.6915	0.3461	Taurus	140.	2, 2	...	...
CoKu Tau/1	M0	3850	0.691	35.	7.	<12	...	>2.52	2.0191	0.4038	Taurus	140.	2, 2	...	...
CoKu Tau/3	M1	3705	0.596	<8.	...	<16	...	...	<0.4615	...	Taurus	140.	2, 2	w	4
CoKu Tau/4	M2	3560	0.579	9.0	2.9	<15	...	>−1.20	0.5192	0.1673	Taurus	140.	2, 2	c	4
CW Tau	K2	4900	1.971	66.	6.	96	8.	−0.88 ± 0.29	3.8075	0.3461	Taurus	140.	2, 2	...	...
CX Tau	M3	3415	0.478	25.	6.	<40	...	>−1.11	1.4422	0.3461	Taurus	140.	2, 2	...	...
CZ Tau	M2	3560	0.579	<9.	...	<30	...	...	<0.5192	...	Taurus	140.	2, 2	c	4
DD Tau	M1	3705	0.596	<42.	...	17	4.	<2.13	2.9786	0.7008	Taurus	140.	2, 2	c	4
DE Tau	M2	3560	0.579	90.	7.	36	5.	2.16 ± 0.37	5.1921	0.4038	Taurus	140.	2, 2	...	...
DF Tau	M1	3705	0.596	8.8	1.9	<25	...	>−2.46	0.5077	0.1096	Taurus	140.	2, 2	c	4
‡ DH Tau	M1	3705	0.596	57.	9.	<57	...	>0.00	3.2883	0.5192	Taurus	140.	2, 2	w	2
DK Tau	K7	4060	0.785	80.	10.	35	7.	1.95 ± 0.56	4.6152	0.5769	Taurus	140.	2, 2	w	4
DL Tau	K7	4060	0.785	440.	40.	230	14.	1.53 ± 0.26	25.3835	2.3076	Taurus	140.	2, 2	...	...
DM Tau	M1	3705	0.596	237.	12.	109	13.	1.83 ± 0.30	13.6725	0.6923	Taurus	140.	2, 2	...	...
DN Tau	M0	3850	0.691	201.	7.	84	13.	2.05 ± 0.37	11.5957	0.4038	Taurus	140.	2, 2	...	...
DO Tau	M0	3850	0.691	258.	42.	136	11.	1.51 ± 0.43	14.8840	2.4230	Taurus	140.	2, 2	...	...
DP Tau	M1	3705	0.596	<10.	...	<27	...	...	<0.5769	...	Taurus	140.	2, 2	c	4
DQ Tau	M0	3850	0.691	208.	8.	91	9.	1.95 ± 0.25	11.9995	0.4615	Taurus	140.	2, 2	sb	3
DR Tau	K5	4350	1.000	533.	7.	159	11.	2.85 ± 0.17	30.7487	0.4038	Taurus	140.	2, 2	...	...
FM Tau	M0	3850	0.691	32.	8.	<36	...	>−0.28	1.8461	0.4615	Taurus	140.	2, 2	...	...
FN Tau	M5	3125	0.150	...	...	<17.5	...	...	<3.0662	...	Taurus	140.	..., 6	...	...
FO Tau	M2	3560	0.579	13.	3.	<14	...	>−0.17	0.7500	0.1731	Taurus	140.	2, 2	c	4
FP Tau	M4	3270	0.286	...	...	<9.3	...	...	<1.6295	...	Taurus	140.	..., 6	...	...
‡ FV Tau	K5	4350	1.000	48.	5.	15	4.	2.74 ± 0.67	2.7691	0.2884	Taurus	140.	2, 2	w	4
‡ FV Tau/c	M4	3270	0.286	<25.	...	<16	...	...	<1.4422	...	Taurus	140.	2, 2	w	4
FX Tau	M1	3705	0.596	17.	3.	<30	...	>−1.34	0.9807	0.1731	Taurus	140.	2, 2	w	4
‡ FY Tau	K7	4060	0.785	<27.	...	16	5.	<1.23	2.8034	0.8761	Taurus	140.	2, 2	...	...
‡ GG Tau A	K7	4060	0.785	1255.	57.	593	53.	1.76 ± 0.24	72.4007	3.2883	Taurus	140.	2, 2	c	4
‡ GH Tau	M2	3560	0.579	15.	3.	<30	...	>−1.63	0.8653	0.1731	Taurus	140.	2, 2	c	4
‡ GK Tau	K7	4060	0.785	33.	7.	<21	...	>1.06	1.9038	0.4038	Taurus	140.	2, 2	...	...
GM Aur	K3	4730	1.538	...	...	253	12.	...	44.3283	2.1025	Taurus	140.	..., 2	...	...
*GM Tau	M6.5	2935	0.074	0.86	0.87	<4.8	...	...	0.0496	0.0502	Taurus	140.	5, 6	...	...
GN Tau	M2.5	3488	0.571	12.	3.	<50	...	>−3.36	0.6923	0.1731	Taurus	140.	2, 2	c	4
GO Tau	M0	3850	0.691	173.	7.	83	12.	1.73 ± 0.35	9.9803	0.4038	Taurus	140.	2, 2	...	...
Haro 6-37	K6	4205	0.850	245.	7.	<88	...	>2.41	14.1340	0.4038	Taurus	140.	2, 2	m	4
HN Tau	K5	4350	1.000	29.	3.	<15	...	>1.55	1.6730	0.1731	Taurus	140.	2, 2	w	4
HO Tau	M1	3705	0.596	44.	6.	<30	...	>0.90	2.5384	0.3461	Taurus	140.	2, 2	w	3
IQ Tau	M1	3705	0.596	178.	3.	87	11.	1.68 ± 0.30	10.2688	0.1731	Taurus	140.	2, 2	...	...
IS Tau	K7	4060	0.785	30.	3.	<20	...	>0.95	1.7307	0.1731	Taurus	140.	2, 2	c	4
IT Tau	K2	4900	1.971	22.	3.	<33	...	>−0.95	1.2692	0.1731	Taurus	140.	2, 2	w	4
J041411+2811	M6.25	2963	0.086	...	...	0.91	0.65	...	0.1594	0.1139	Taurus	140.	..., 4	...	...
J043814+2611	M7.25	2838	0.049	...	...	2.29	0.75	...	0.4012	0.1314	Taurus	140.	..., 4	...	...
J043903+2544	M7.25	2838	0.049	...	...	2.86	0.76	...	0.5011	0.1332	Taurus	140.	..., 4	...	...
J044148+2534	M7.75	2753	0.035	...	...	2.64	0.64	...	0.4626	0.1121	Taurus	140.	..., 4	...	...
*J044427+2512	M7.25	2838	0.049	9.85	0.76	7.55	0.89	0.63 ± 0.33	0.5682	0.0438	Taurus	140.	5, 4	...	...
JH 112	K6	4205	0.850	30.	10.	<18	...	>1.20	1.7307	0.5769	Taurus	140.	2, 2	m	4
JH 223	M2	3560	0.579	<7.	...	<19	...	...	<0.4038	...	Taurus	140.	2, 2	w	4
KPNO Tau 6	M8.5	2555	0.020	...	...	−0.66	0.79	...	−0.1156	0.1384	Taurus	140.	..., 4	...	...
KPNO Tau 7	M8.25	2633	0.025	...	...	0.70	0.88	...	0.1226	0.1542	Taurus	140.	..., 4	...	...
KPNO Tau 12	M9	2400	0.014	...	...	−0.92	0.70	...	−0.1612	0.1226	Taurus	140.	..., 4	...	...
LkCa-15	K5	4350	1.000	428.	11.	167	6.	2.22 ± 0.10	24.6912	0.6346	Taurus	140.	2, 2	...	...
MHO 5	M6	2990	0.096	...	...	<9.0	...	...	<1.5769	...	Taurus	140.	..., 6	...	...
RW Aur	K3	4730	1.538	79.	4.	42	5.	1.49 ± 0.30	4.5575	0.2308	Taurus	140.	2, 2	w	4
RY Tau	K1	5080	2.481	560.	30.	229	17.	2.10 ± 0.22	32.3063	1.7307	Taurus	140.	2, 2	...	...
St 34	M3	3415	0.478	<11.	...	<15	...	...	<0.6346	...	Taurus	140.	2, 2	sb, w	3, 4
SU Aur	G2	5860	3.267	74.	3.	<30	...	>2.12	4.2690	0.1731	Taurus	140.	2, 2	...	...
T Tau	K0	5250	2.847	628.	17.	280	9.	1.90 ± 0.10	36.2292	0.9807	Taurus	140.	2, 2	m	3
UX Tau	K2	4900	1.971	173.	3.	63	10.	2.38 ± 0.38	9.9803	0.1731	Taurus	140.	2, 2	m	4
UY Aur	K7	4060	0.785	102.	6.	29	6.	2.96 ± 0.51	5.8844	0.3461	Taurus	140.	2, 2	w	4
UZ Tau	M1	3705	0.596	560.	7.	172	15.	2.78 ± 0.21	32.3063	0.4038	Taurus	140.	2, 2	m	3, 4
V410 Anon 13	M5.75	3024	0.109	...	...	<9.0	...	...	<1.5769	...	Taurus	140.	..., 6	...	...
‡ V710 Tau	M1	3705	0.596	152.	6.	60	7.	2.19 ± 0.29	8.7689	0.3461	Taurus	140.	2, 2	w	3, 4
V836 Tau	K7	4060	0.785	74.	3.	37	6.	1.63 ± 0.39	4.2690	0.1731	Taurus	140.	2, 2	...	...
V892 Tau	A0	9520	3.791	638.	54.	234	19.	2.36 ± 0.28	36.8061	3.1153	Taurus	140.	2, 2	w	3
‡ V955 Tau	M0	3850	0.691	14.	2.	<19	...	>−0.72	0.8077	0.1154	Taurus	140.	2, 2	c	4
VY Tau	M0	3850	0.691	<10.	...	<17	...	...	<0.5769	...	Taurus	140.	2, 2	c	4
*2M1207-3932	M8	2710	0.035	0.07	1.75	...	...	...	0.0006	0.0141	TWA	52.4	5, ...	c	5
Hen3-600A	M4	3270	0.218	65.	5.	...	...	...	0.4783	0.0368	TWA	50.	7, ...	c, sb	6
*SSPM1102-3431	M8	2710	0.035	−0.94	1.57	...	...	...	−0.0084	0.0141	TWA	55.2	5, ...	...	...
*TWA 30A	M5	3125	0.114	1.70	1.63	...	...	...	0.0088	0.0085	TWA	42.	5, ...	w	7
*TWA 30B	M4	3270	0.218	0.013	1.86	...	...	...	6.7e–5	0.0097	TWA	42.	5, ...	w	7
TW Hya	K6	4205	0.949	1450.	310.	874	54.	1.19 ± 0.52	13.3840	2.8614	TWA	56.	8, 8	...	...

Notes. ^a"*" marks sources observed in this paper, at 850 μm on SCUBA-2. "‡" marks sources with more binarity information in footnote 21 of the main text. ^bMeasured fluxes and 1σ errors, or 3σ upper limits (for sources from the literature where the measured value is not cited). ^cM_{d, [850]} if either a measured value of the 850 μm flux is available, or the upper limit on M_{d, [850]} is smaller than the upper limit on M_{d, [1300]}; otherwise M_{d, [1300]}. Computed using Equation (9). These are the values used in our Bayesian analysis. ^dMeasured values of α with 1σ error bars, or 3σ upper or lower limits. See Equation (8) and discussion in Section 5. ^eReferences for fluxes at 850 μm (first reference) and 1.3 mm (second reference). (1) AW07a, and references therein; (2) AW05, and references therein; (3) Klein et al. 2003; (4) Scholz et al. 2006; (5) This study (these sources are also marked with an asterisk); (6) Schaefer et al. 2009; (7) Zuckerman 2001; and (8) Weintraub et al. 1989. ^fc = close binary (sep < 100 AU), w = wide binary (sep ⩾ 100 AU), sb = spectroscopic binary, and m = multiple. ^gReferences for multiplicity information. (1) R10a, and references therein; (2) Reipurth & Zinnecker 1993; (3) AW05, and references therein; (4) Kraus et al. 2011, and references therein; (5) Chauvin et al. 2005; (6) Andrews et al. 2010, and references therein; and (7) Looper et al. 2010a. See also footnote 21 in the main text, for sources IDs marked with ‡. ^hFor TW Hya, this is the flux (and 1σ error) at 1.1 mm (Weintraub et al. 1989); we have cited it here for completeness, since close to 1.3 mm.

Download table as:  ASCIITypeset images: 1 2 3 4

2.2. SCUBA-2 Observations and Fluxes

We calculate stellar masses, radii, and effective temperatures for all our sources using previously determined spectral types, a single assumed age for all sources in a given star-forming region or association, and the predictions of theoretical evolutionary tracks. Spectral types are converted to effective temperatures using the conversion scheme in Kenyon & Hartmann (1995) for types M0 and earlier, and the scheme devised by Luhman et al. (2003) for later types. Objects in Taurus and ρ Oph are assumed to have an age of ∼1 Myr, and those in the TWA are assumed to have an age of ∼10 Myr. The derived temperatures are then compared to solar-metallicity evolutionary model predictions for the assumed ages to infer stellar masses and radii. We use (1) the models of Siess et al. (2000) for stars with mass >1.4 M_☉; (2) the models of Baraffe et al. (1998) for stellar masses 0.08–1.4 M_☉, specifically, models with a mixing length to pressure scale height ratio of α_mix = 1.0 for 0.08–<0.6 M_☉ and models with α_mix = 1.9 (value required to fit the Sun) for 0.6–1.4 M_☉ (see discussion in Baraffe et al. 2002); and (3) the "Dusty" models of Chabrier et al. (2000) for masses <0.08 M_☉ (which are mid- to late-M types at ages of a few mega-years, the spectral-type regime where photospheric dust starts to become important).

We divide our sample into four stellar mass bins: intermediate-mass stars (>1–4 M_☉), solar-type stars (0.3–1 M_☉), VLMS (0.075–<0.3 M_☉), and BDs (≲0.075 M_☉). Within each, we adopt the median stellar mass as a representative of its class: 2.5 M_☉ for intermediate-mass stars (⇒ R_* ≈ 4 R_☉ and T_* ≈ 5000 K at 1 Myr); 0.75 M_☉ for solar types (2 R_☉, 4000 K; the usual parameters adopted for solar-type cTTs in the literature); 0.2 M_☉ for VLMS (1.5 R_☉, 3200 K); and 0.05 M_☉ for BDs (0.55 R_☉, 2850 K). These fiducial masses will serve to illustrate the trends in disk mass and grain growth implied by our disk models for each mass bin. In Taurus, we have 11 intermediate-mass stars, 49 solar types, 9 VLMS, and 11 BDs; in ρ Oph, we have 9, 32, 7, and 0 in the same bins; and in the TWA, we have 0, 1, 3, and 2, respectively.

The theory of disk SEDs is outlined in Appendix A. The main assumptions and resulting equations relevant to our analysis are as follows.

We assume that the disk surface density and temperature have power-law radial profiles¹¹:

$$\begin{equation} \Sigma (r) = \Sigma _0\left(\frac{r}{r_0}\right)^{-p}, \quad T(r) = T_0\left(\frac{r}{r_0}\right)^{-q}, \end{equation} \tag{ 1 }$$

where Σ₀ and T₀ are the values at the disk inner edge r₀. We further assume that the opacity is a power law in frequency:

$$\begin{equation} \kappa _\nu = \kappa _f \,\left(\frac{\nu }{\nu _f}\right)^\beta , \end{equation} \tag{ 2 }$$

where κ_f is the opacity at some appropriate fiducial frequency ν_f. The spectral index of the emission is defined as

$$\begin{equation} \alpha \equiv \frac{d({\rm ln}F_\nu)}{d({\rm ln}\,\nu)}, \end{equation} \tag{ 3 }$$

where F_ν is the disk flux density at frequency ν measured by an observer.

In the RJ limit, the flux density is given by a polynomial expression (B90), derived in Appendix A. In this case, the spectral index reduces to

$$\begin{equation} \alpha \approx 2 + \frac{\beta }{1+\Delta }, \end{equation} \tag{ 4 }$$

where Δ is the ratio of optically thick to optically thin emission (see Appendix A).¹² For grains much larger than the wavelength observed, the opacity κ_ν is independent of the frequency (i.e., β ∼ 0), yielding α ≈ 2 regardless of the disk optical thickness. Conversely, optically thick emission (Δ → ∞) also implies α ≈ 2, independent of β and hence grain size. In the optically thin limit (Δ → 0), on the other hand, α ≈ 2 + β.

If the RJ approximation is not valid, then the general expression for the flux density is

$$\begin{eqnarray} F_\nu & \approx & \nu ^3 \left(\frac{4\pi h}{c^2}\right) f_0 \left(\frac{{\rm cos}\,i}{D^2}\right) \left[ \frac{(2-p){\bar{\tau }}_\nu }{2} \,R_d^p \right]\nonumber\\ &&\times\left[\int _{r_1}^{R_d} \frac{r^{1-p}}{{\rm exp}[h\nu /kT(r)]-1}\,dr \right] (1+\Delta) \end{eqnarray} \tag{ 5 }$$

for a disk with inner radius r₀, outer radius R_d, situated at a distance D from the observer, and inclined at an angle i relative to the line of sight (so that i = 90° for an edge-on orientation). r₁ is the radius at which the emission at ν changes from optically thick to thin, ${\bar{\tau }}_\nu$ is the average optical depth in the disk ($\equiv (\kappa _\nu M_d)/(\pi R_d^2\,{\rm cos}\,i)$, where M_d is the disk mass), and f₀ ∼ 0.8 is a correction factor (see Appendix A). The product of all the terms outside the last parentheses is the optically thin contribution to the flux density; the ratio Δ of the optically thick to thin contributions is given by

$$\begin{equation} \Delta \equiv \left(\frac{2}{(2-p){\bar{\tau }}_\nu \, R_d^p} \right) \left[\frac{\int _{r_0}^{r_1} \frac{r}{{\rm exp}[h\nu /kT(r)]-1}\,dr}{\int _{r_1}^{R_d} \frac{r^{1-p}}{{\rm exp}[h\nu /kT(r)]-1}\,dr}\right]. \end{equation} \tag{ 6 }$$

Equations (5) and (6) must be evaluated numerically, and the spectral index α computed directly from its definition, Equation (3). In particular, even if the opacity remains a power law in frequency, α does not reduce to the simple form of Equation (4), and will generally be smaller in the sub-mm/mm (because the spectrum is flatter) than the RJ asymptotic values of 2 and 2 + β for optically thick and thin disks, respectively. We show in Section 6 that the RJ limit is not ideal at 850 μm and 1.3 mm for our sample, and instead use the generalized forms of F_ν and α to estimate disk masses and grain growth in Section 7, via the technique described below.

For sources observed at both 850 μm and 1.3 mm, we have an estimate of F_ν as well as α. Using the foregoing equations, we wish to investigate their disk masses and grain properties. The latter, however, are specified by three unknown parameters—M_d, κ_f, and β—while F_ν and α represent only two independent observables. As such, we can only derive β and the product κ_fM_d from the observed SED (e.g., Natta et al. 2004; R10a, b, note that only this product enters the expression for the observed flux, Equation (5), through the quantity ${\bar{\tau }}_\nu$).

Even this, of course, still requires specifying the other six parameters that F_ν and α depend on the disk inclination i, the disk inner and outer radii r₀ and R_d, the surface density power-law index p, the temperature power-law index q, and the temperature at the disk inner edge T₀. A rigorous determination of these entails spatially resolved photometry and spectroscopy from optical to mm wavelengths, which have not been obtained for the vast majority of sources. Instead, we adopt fiducial parameters guided by theory and current observations, which suffices to estimate the broad trends in κ_fM_d and β. The detailed rationale behind our choices is presented in Appendix B; the final adopted values are

$$\begin{eqnarray} i &=& 60^\circ , \,\,r_0=5\ R_{\ast } , \,\,R_d=100\,{\rm AU}, \,\,p=1, \nonumber\\ q &=& 0.58, \,\,T_0=880(T_{\ast } /4000\,{\rm K})\,{\rm K}. \end{eqnarray} \tag{ 7 }$$

Additionally, we stipulate that the disk temperature cannot fall below a minimum value T_min ≡ 10 K, the temperature declines with radius as a power law until this value is reached, and remains at a constant 10 K at larger radii. This is justified in Appendix A, on the basis of heating by the interstellar radiation field (ISRF). Note that such a broken power law does not do any violence to our generalized Equations (5) and (6), where the precise form of T(r) is unspecified. Finally, we also test the consequences of varying R_d over the range 10–300 AU and p over the range 0.5–1.5, for different stellar mass bins. With these six parameters specified, we determine κ_fM_d and β by comparing the observed α and F_ν to the predictions of the generalized equations in Section 4, as follows.

With only two observed wavelengths, 850 μm and 1.3 mm, the spectral index defined by Equation (3) reduces to

$$\begin{equation} \alpha \equiv \frac{{\rm ln}\,F_{[850]} - {\rm ln}\,F_{[1300]}}{{\rm ln}\,1300 - {\rm ln}\,850}. \end{equation} \tag{ 8 }$$

We use this expression to calculate α for both the data and the SED models they are compared to. Above, and henceforth, the subscript [λ] denotes a quantity evaluated at a frequency ν = c/λ. For sources detected at 1.3 mm but not at 850 μm, we compute the 3σ upper limits on α by replacing F_[850] above by its 3σ upper limit. Similarly, for sources detected at 850 μm but not at 1.3 mm, we find 3σ lower limits on α by replacing F_[1300] by its 3σ upper limit. Sources with only upper limits at both 850 μm and 1.3 mm, or those not observed at all at one or the other wavelength, are excluded from this analysis.¹³

Next, following common practice, we choose to normalize our opacity power law, Equation (2), at a fiducial frequency ν_f ≡ 2.3 × 10¹¹ Hz, corresponding to λ_f = 1300 μm (1.3 mm). We denote κ_f, the (a priori unknown) opacity at this frequency, explicitly as κ_[1300]. We further denote the commonly adopted value of this fiducial opacity by $\tilde{\kappa }_{[1300]} \equiv 0.023$ cm² g⁻¹ (e.g., B90). Note that the latter value assumes an interstellar medium (ISM) like gas-to-dust mass ratio of 100:1.

Finally, we define the apparent disk mass M_{d, ν} as the mass derived from the observed flux at some frequency ν, assuming that the emission (1) is optically thin, (2) arises from an isothermal region of the disk with known temperature $\tilde{T}$, and (3) is due to material with a known opacity $\tilde{\kappa }_\nu$. Thus,

$$\begin{equation} M_{d,\nu } \equiv \frac{F_{\nu } D^2}{\tilde{\kappa }_{\nu } B_{\nu }(\tilde{T})}. \end{equation} \tag{ 9 }$$

We adopt $\tilde{T} = 20$ K and $\tilde{\kappa }_\nu = \tilde{\kappa }_{[1300]}$ (ν/2.3 × 10¹¹ Hz)^β, with $\tilde{\kappa }_{[1300]}$ as defined above and β = 1. We will write M_{d, ν} explicitly as M_{d, [850]} or M_{d, [1300]} when F_ν is the flux density at 850 μm or 1.3 mm. With $\tilde{T}$ and $\tilde{\kappa }_\nu$ fixed, M_{d, ν} is simply a proxy for the specific disk luminosity F_νD². The advantage of using this formulation is that M_{d, ν} (evaluated at 850 μm or 1.3 mm, with $\tilde{T}$ and $\tilde{\kappa }_\nu$ similar or identical to our adopted values) is widely used as a simplistic estimate of the disk mass (e.g., AW05; AW07a; Klein et al. 2003; Scholz et al. 2006). This allows us to directly compare our results to literature values for disk masses.

Turning now to the theoretical models, we define the opacity-normalized disk mass as

$$\begin{equation} M_d^{\kappa} \,\equiv \, \frac{\kappa _{[1300]} M_d}{\tilde{\kappa }_{[1300]}} \,=\, \left(\frac{\kappa _{[1300]}}{0.023\ {\rm cm}^2\ {\rm g}^{-1}}\right)M_d, \end{equation} \tag{ 10 }$$

$M_d^{\kappa}$ is merely a scaled proxy for the real variable κ_[1300]M_d, but is useful to employ in lieu of the latter as a more intuitive quantity (e.g., Natta et al. 2004), specifically, $M_d^{\kappa}$ equals the real disk mass M_d if the true opacity κ_[1300] equals the fiducial value $\tilde{\kappa }_{[1300]}$. Without knowing the true value of the opacity, $M_d^{\kappa}$ is the closest we can get to the real disk mass.

Our analysis now is straightforward. We first show, in Section 6, that the RJ approximation is not ideal at the disk temperatures expected in our sample. We thus use the generalized equations (5), (6), and (3) (with the latter approximated by Equation (8)). From these, we calculate the F_ν, Δ, and α predicted by our six fixed disk parameters and a physically plausible range of $M_d^{\kappa}$ and β; we further convert the predicted F_ν to a predicted M_{d, ν} using Equation (9). Comparing these theoretical M_{d, ν} and α values to the values derived from the observed fluxes enables us to (1) gauge what fraction of the emission is optically thin; and (2) estimate the true $M_d^{\kappa}$ and β when the emission is predominantly optically thin, or put lower limits on $M_d^{\kappa}$ (and no constraints on β) when it is optically thick.

Basically, all we are doing is calculating (within the context of our fiducial disk parameters) what the emitted flux and α should be for any specified $M_d^{\kappa }$ and β, and comparing these to the observed flux and α to infer what the $M_d^{\kappa }$ and β for a source really are; i.e., a simple inversion. The only potential confusion for the reader is that, instead of directly using flux in these comparisons, we use its scaled proxy M_{d, ν}, a quantity which is often cited in the literature and thus useful to have at hand. However, the various assumptions that go into calculating M_{d, ν} (fixed T = 20 K, fixed β = 1) have no effect on our results: the predicted and observed fluxes are converted to predicted and observed M_{d, ν} using exactly the same multiplicative factors, so the end result is precisely the same as simply using flux instead.

Figure 1 shows the α predicted by the generalized equations for our fiducial disk model, Equation (7), for the median stellar temperature and radius within each of our four stellar mass bins: intermediate-mass stars, solar-type stars, VLMS, and BDs. The abscissae of the plots show the corresponding minimum temperature in the disk (i.e., the temperature at the outer disk radius R_d = 100 AU); in VLMS and BDs, this value levels out at our fixed lower limit of 10 K. Note that the plotted α refers to emission integrated over the entire disk, not just from R_d. The opacity-normalized disk mass, $M_d^{\kappa}$, is either fixed at a very large value, so that the entire disk is optically thick (i.e., the transition radius r₁ from Equation (A4) formally satisfies r₁ ≫ R_d; top panel of Figure 1), or fixed at a very small value, so that the entire disk is optically thin (r₁ ≪ r₀; bottom panel). β is varied from 0 to 2.

Zoom In Zoom Out Reset image size

The top panel shows that, as expected, α in the optically thick limit is independent of the opacity index β. However, the asymptotic RJ value in this case, α = 2, is not achieved for any of our stars. For the relatively hot disks around intermediate-mass stars, the deviation from RJ is quite small, ∼0.25. For the cooler disks around solar-type stars to BDs, on the other hand, the SED is markedly flatter, yielding an α significantly smaller—by 0.5–0.8—than the RJ expectation. Similarly, the bottom panel shows that while α does depend on β in the optically thin limit, as expected, the asymptotic RJ value in this limit, α = 2 + β, is not reached. Again, the deviation is relatively small for intermediate-mass stars (∼0.2), but appreciable in VLMS and BDs (∼0.4–0.6).

It is often assumed that the RJ limit applies at ∼mm wavelengths for disk temperatures ≳15 K. Figure 1 shows that this is not very accurate, with a discernible α deviation of ∼0.2 even for intermediate-mass stars, where the minimum disk temperature is ∼20 K (and most of the disk is considerably hotter still). More importantly, disks significantly larger than 100 AU around solar-type and intermediate-mass stars, as are often observed, as well as even 100 AU disks around VLMS/BDs (like those plotted in Figure 1), may have temperatures down to ∼10 K in their outer regions, where most of the mass resides. Figure 1 shows that the RJ approximation is severely strained under these circumstances.¹⁴ The danger in assuming RJ values is that a low observed α from an optically thin disk will lead one to infer a spuriously low β (i.e., too much grain growth), and/or a spuriously high optically thick contribution. Equally importantly, a number of disks appear to have α < 2 (as we shall shortly see); these cannot be explained at all under the RJ assumption.

Even in the general non-RJ case, however, we see that α remains very sensitive to β for optically thin emission, with a roughly linear relationship between the two when all other disk parameters are fixed (i.e., analogous to the RJ case, except that now the precise value of α depends on the disk temperature as well, in both optically thin and thick limits). Thus, α from the general equations can still be used to probe grain growth, as we do next.

7.1. Spectral Slopes and Apparent Disk Masses from the Observed Fluxes

We first discuss some general trends in the spectral slopes (α) and apparent disk masses (M_{d, ν}) that we derive from the observed fluxes, before comparing to model predictions.

Figure 2 shows the α for our sample as a function of stellar mass, for sources observed at both 850 μm and 1.3 mm and detected at at least one of these two wavelengths. Sources detected at both are shown as filled circles with 1σ error bars,¹⁵ while 3σ upper and lower limits are plotted as downward and upward triangles, respectively. Note that the plotted sources are mainly at M_* ≳ 0.5 M_☉, since most lower mass stars and BDs have either not been observed at both wavelengths or are undetected at both (see Table 1). The lower limits are also concentrated at the lower end of this stellar mass range, where the disks are still bright enough to be detected at 850 μm but too faint for 1.3 mm; this reflects the empirical fact that disks become fainter with diminishing M_*.

Zoom In Zoom Out Reset image size

For the sub-sample detected at both wavelengths, we find a mean spectral index of 〈α〉 = 1.98 ± 0.06 (1σ error; mean and associated 3σ error plotted in Figure 2), in agreement with the average α derived by AW05.¹⁶ While our mean is driven predominantly by stars ≳ 0.5 M_☉, the four VLMS/BDs with measured α are consistent with this value. Also, while we have ignored upper and lower limits in this calculation, Figure 2 shows that all but two of the lower limits and two out of the four upper limits are compatible, within the errors, with this mean; including the four deviant points negligibly affects the outcome.

In the upper panels of Figure 3, we plot M_{d, [850]} and M_{d, [1300]} for our sample (derived via Equation (9), using the distances D discussed in Section 3), as a function of stellar mass. Going from BDs to solar-type stars, there is a clear trend of M_{d, ν} increasing on average with higher M_*; we return to this in Section 8. For now, note that this simply reflects an increase in the average emitted specific disk luminosity (F_νD²) with stellar mass. Conversely, there is a decline in the upper envelope of M_{d, ν} from the solar-type to intermediate-mass stars, going approximately as $[M_{d,\nu }]_{{\rm max}} \propto M_{\ast } ^{-1/2}$ (shown by a dotted line in Figure 3). While this relationship is defined by only a handful of stars, we can at least state that there is no evidence of [M_{d, ν}]_max increasing with M_* for these stars; we discuss this in Section 7.3.

Zoom In Zoom Out Reset image size

In the lower panels of Figure 3, we plot the corresponding M_{d, [850]}/M_* and M_{d, [1300]}/M_* as a function of stellar mass. Excluding the four upper limits among the VLMS/BDs in TWA (green downward triangles), the mean M_{d, ν}/M_* appears roughly constant from BDs to solar types, reflecting the increase in M_{d, ν} with M_* noted above. Conversely, the decline in the upper envelope from solar-type to intermediate-mass stars is now more pronounced, as expected: $[M_{d,\nu }]_{{\rm max}} \propto M_{\ast }^{-1/2}$ in the upper panel translates to a steeper slope $[M_{d,\nu }/M_{\ast }]_{{\rm max}} \propto M_{\ast }^{-3/2}$ in the lower one. The limiting mass ratio above which gravitational instabilities set in, M_d/M_* ∼ 0.1 (e.g., Lodato et al. 2005), is also marked for later reference.

Finally, since the M_{d, ν} results at 1.3 mm are nearly identical to those at 850 μm, we focus henceforth on the latter for clarity. To estimate the true grain growth and opacity-normalized disk masses, we must now compare these M_{d, [850]} and α to the model predictions.

7.2. Model Predictions

Our model calculations are carried out for four fiducial M_* (with corresponding R_* and T_* at 1 Myr), representing median values within the four stellar mass bins in our sample (see Section 3): 2.5 M_☉ (intermediate-mass stars), 0.75 M_☉ (solar types), 0.2 M_☉ (VLMS), and 0.05 M_☉ (BDs). For each of these M_*, models are constructed for β = 0–2 (lowest possible to standard ISM) in steps of 0.1, and $M_d^{\kappa }$ = 10⁻⁶–1 M_☉ in steps of 1 dex (sufficient to capture the observed range in M_{d, [850]}, as illustrated shortly). Furthermore, if $M_d^{\kappa }$ represents the true disk mass M_d (i.e., if the true opacity κ_[1300] equals the fiducial value 0.023 cm² g⁻¹), then the disk would nominally become gravitationally unstable above $M_{d,GI}^{\kappa }\equiv 0.1\,M_{\ast }$. We explicitly include this value in our models for each stellar bin, for reasons explained further below. We examine a range of R_d over 10–300 AU, and also vary the surface density exponent p over 0.5–1.5 (instead of fixing it at 1) in some cases. The rest of the parameters are as specified in Equation (7). Inserting these into Equations (5), (6), and (8), we compute the theoretical F_ν, Δ, and α, and further convert the predicted F_[850] to a predicted M_{d, [850]} with Equation (9).

The model predictions are plotted in Figures 4–7. Each figure is for a fixed M_*; the top and bottom sets of plots in each figure correspond to two illustrative values of R_d. The four panels in each plot show (from left to right and top to bottom) the following.

• 1.  
  The predicted F_[850] (arbitrarily scaled to the distance of Taurus for illustration) versus β, for the range of β and $M_d^{\kappa}$ examined. For optically thin emission, the flux density increases with both $M_d^{\kappa}$ and β. As the optically thick contribution rises with $M_d^{\kappa}$, F_[850] becomes less sensitive to both parameters, until finally, in the thick limit, the flux density is independent of β and $M_d^{\kappa}$, and saturates at a level set by the disk temperature and radial extent.
• 2.  
  The predicted α versus M_{d, [850]} (the latter calculated from the predicted F_[850]; independent of distance), for our grid of β and $M_d^{\kappa}$. As discussed in Section 6, α increases roughly linearly with β for optically thin disks, and saturates to a fixed value when the emission becomes optically thick. Unlike in the RJ limit, however, the value of α in both cases depends on the disk temperature (and thus on T_*, or equivalently M_*, as well as on R_d).Overplotted in this panel for comparison are the α and M_{d, [850]} derived from our observed fluxes, for sources in the relevant stellar mass bin (objects with upper/lower limits in α are excluded for clarity). The horizontal dashed line marks the mean index for the full sample, 〈α〉 ≈ 2. The vertical dashed line denotes the maximum observed M_{d, [850]} (equivalently, maximum observed 850 μm flux) in our data for the relevant stellar mass bin. This observed upper limit sets a lower boundary on the size R_d of the brightest disks, as follows. For a given surface density and temperature profile, a disk of a specified size has a maximum allowed flux, which is its optically thick limit, and this maximal value increases with disk size. Thus, a minimum disk size is required to explain the maximum observed flux (or equivalently, maximum observed M_{d, [850]}).¹⁷ This minimum R_d changes as a function of the adopted density and temperature profiles.An additional constraint may be set by the disk mass. The optically thick limit for a disk of a given size is reached above a threshold disk mass, and this threshold value increases with disk size. Thus, the minimum R_d derived above also corresponds to a minimum disk mass. However, one expects the gravitational instability limit, M_d ∼ 0.1 M_*, to set a reasonable upper limit to the allowed disk mass. Thus, the minimum mass corresponding to the derived lower bound on R_d must not greatly exceed the instability limit. If it does, then some other parameter must be changed (e.g., p or q) to resolve the discrepancy.This analysis is complicated by the fact that we do not know the real opacity κ_[1300], so we cannot associate our disk models with a true mass M_d, but only with its opacity-normalized counterpart $M_d^{\kappa }$. Nevertheless, the mass constraint described above can still be implemented by setting plausible limits on how far the true opacity may deviate from our fiducial value $\tilde{\kappa }_{[1300]} = 0.023$ cm² g⁻¹. Using detailed self-consistent dust models, R10b show that the fiducial opacity may be an underestimation by up to a factor of 10 for β ∼ 2, and an overestimation by the same factor for β ∼ 0. If the true opacity equaled the fiducial one, then $M_{d,GI}^{\kappa }$ (≡ 0.1 M_*, as defined earlier) would represent a disk that was actually at the unstable limit; the R10b analysis suggests that in fact, the instability limit may lie anywhere in the range 0.1–10 $M_{d,GI}^{\kappa }$. As such, in determining the minimum disk radius, we must verify that the corresponding minimum $M_d^{\kappa }$ does not exceed at least the upper limit 10 $M_{d,GI}^{\kappa }$.
• 3.  
  Predicted Δ versus M_{d, [850]}, explicitly showing the relative contributions of optically thin and thick emission for each $M_d^{\kappa}$ and β. The emission flux, and thus the corresponding M_{d, [850]}, saturates to a fixed value in the optically thick limit.
• 4.  
  Predicted ratio of M_{d, [850]} to $M_d^{\kappa}$, versus M_{d, [850]}. For hot stars with optically thin disks, the assumption of isothermal emission at $\tilde{T} = 20$ K in deriving M_{d, [850]} can cause the latter to exceed the true value $M_d^{\kappa}$, since a good fraction of the emission is from radii hotter than 20 K for these sources. For cooler stars/BDs in the optically thin limit, the same assumption can make M_{d, [850]} an underestimation of $M_d^{\kappa}$, due to emission from radii cooler than 20 K. In optically thick disks, M_{d, [850]} saturates and always represents a lower limit on the true $M_d^{\kappa}$.The dashed line in this panel plots the [$M_d^{\kappa}$, β] locus, derived from point (2) above, corresponding to the mean spectral index of the sample 〈α〉 ≈ 2. This reveals the average trend in $M_{d,[850]}/M_d^{\kappa}$ among our sources in each stellar mass bin.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

7.3. Implications of Model Comparisons to Data

In the bottom two panels of Figure 3, the mean values of both M_{d, [850]}/M_* and M_{d, [1300]}/M_* appear roughly constant with M_*, among the approximately coeval Taurus and ρ Oph populations. This apparently flat distribution of M_{d, ν}/M_* has been commented on in previous work as well; it suggests that on average, M_{d, ν}/M_* ∼ 10⁻² (e.g., Scholz et al. 2006). However, the presence of a large number of upper limits, especially among the VLMS and BDs, makes the veracity of this claim hard to judge by eye alone. Instead, we use a Bayesian analysis to test this. The technique is described in Appendix C, and the results are discussed below.

Zoom In Zoom Out Reset image size

We emphasize that we only analyze the distribution of the apparent disk mass M_{d, ν}, and not of the true disk mass M_d, or even of the opacity-normalized mass $M_d^{\kappa}$. As we have discussed, translating the first into either of the latter two quantities requires knowledge of a number of disk parameters, which are unknown for most of our sample. As such, our discussion above of various broad trends in $M_d^{\kappa}$ and M_d suggested by the data is the best we can do; precise determination of these two quantities for all the individual stars, required for a statistical investigation of the underlying distribution, is not currently possible. Nevertheless, the statistics of M_{d, ν} alone are still valuable as an initial indicator of the possible behavior of the true disk mass. Equally importantly, the analysis serves to illustrate the Bayesian techniques that can be applied to the true M_d distribution when it is derived in the future, as well as to any other distribution that is both noisy and plagued by upper limits.

To begin with, we combine our M_{d, [850]} and M_{d, [1300]} estimates to get the largest possible sample of apparent disk masses. Specifically, we use M_{d, [850]} if available, otherwise M_{d, [1300]} (see Table 1). Furthermore, our sample includes a number of known binaries and higher-order multiples (Table 1). Within our Taurus and ρ Oph sub-samples, most such systems have not been resolved in the sub-mm/mm data presented here.²¹ As such, the M_{d, ν} we derive corresponds to the total apparent disk mass in these systems; we do not know how this is partitioned between the components.²² In all these cases, we assume that the disk material is present solely around the primary, i.e., we adopt M_{d, ν}/M_* ≡ M_{d, ν}/M_primary. This is because binarity studies are very much incomplete for our sample (especially for the ρ Oph sources and the VLMS/BDs), many of the objects assumed to be single here have not been subject to as wide a companion parameter search as the identified binaries, and many have not been examined for multiplicity at all. For unidentified binaries in our sample, the M_{d, ν}/M_* we calculate implicitly corresponds to assigning the total disk mass entirely to the primary²³; doing the same for the known binaries/multiples is thus necessary for uniformity.

The final combined sample is plotted in Figure 9. The full sample is shown in the top panel, and the "single" objects (i.e., sample with known binaries/multiples removed) in the bottom panel; some of the latter may have as yet unknown companions. Note that most of the VLMS/BDs plotted as 3σ upper limits have actual measured values at <3σ significance (Table 1). Our plotting convention is simply to facilitate visual comparison to objects from the literature for which only 3σ upper limits in flux have been published. The actual measurements, where available, are used in our subsequent Bayesian analysis.

Zoom In Zoom Out Reset image size

For the analysis, we first divide our full sample into four populations: (1) Taurus solar-type stars (49 sources); (2) ρ Oph solar-type stars (32 sources); (3) all (Taurus + ρ Oph) VLMS/BDs (27 sources); and (4) all (Taurus + ρ Oph) intermediate-mass stars (20 sources). Taurus and ρ Oph objects, which are roughly coeval, are lumped together in the VLMS/BD and intermediate-mass bins to increase the sample sizes therein. The small and much older set of six TWA objects is excluded from this analysis. The population of Taurus solar-type stars, which is the best constrained (in that it includes the most data and the least number of upper limits), serves as a baseline against which the other populations are compared.

The Taurus solar types have also been subject to more thorough binarity surveys than the other populations. To evaluate the effects of multiplicity, therefore, we also compare the full sample of Taurus solar types to the single Taurus solar types (18 sources); the Taurus solar-type close binaries (projected component separation <100 AU; 13 sources) to the singles; and the Taurus solar-type wide binaries (⩾100 AU; 12 sources) to the singles.

For each population, we assume that M_{d, ν}/M_* is described by an underlying lognormal distribution specific to that population. Our Bayesian analysis then reveals the probability distributions for the mean (μ) and standard deviation (s) of this lognormal in each case.

Finally, we include only the Gaussian photon noise in the observed fluxes in our error analysis, and not the systematic calibration uncertainties, nor the uncertainties in the assumed values for disk parameters used to calculate M_{d, ν} ($\tilde{T}$ = 20 K, β = 1, $\tilde{\kappa }_{[1.3]}$ = 0.023 cm² g⁻¹, and gas-to-dust ratio = 100:1; see Section 5), nor the uncertainties in the M_* inferred in Section 2. Including the additional calibration uncertainty would mostly only strengthen our results, as pointed out at appropriate junctures below. Moreover, our lack of knowledge about the true disk parameters for a large fraction of our sources prevents us from deriving M_d or $M_d^{\kappa}$ for individual stars in the first place, which is why we concentrate here only on the apparent disk mass (we do point out the effect of using β < 1, as seems appropriate for many of our stars, on our results). Similarly, while uncertainties and systematic errors are undoubtedly present in the evolutionary tracks used to calculate M_*, these are difficult to quantify precisely with our present knowledge. We therefore consider only the Gaussian noise in the flux here.

Figures 10–15 show the outcomes of our analysis. In Figures 10–12, we inter-compare the results for various sub-samples of the Taurus solar-type population (full sample, singles, close binaries, and wide binaries). In Figures 13–15, we compare the result for the full sample of Taurus solar types to that for each of the other populations (ρ Oph solar types, all intermediate masses, and all VLMS/BDs). In each figure, the three panels show (clockwise from top right): (1) the full two-dimensional (2D) probability distribution of the lognormal parameters μ (mean) and s (standard deviation) for a given population of stars, with the contours enclosing 68.27% (1σ), 95.45% (2σ), and 99.73% (3σ) of the distribution; (2) the 1D probability distribution of the mean μ, marginalized (integrated) over all standard deviations s (i.e., the distribution of μ independent of the precise value of s), with contours again at 1σ–3σ; and similarly (3) the 1D probability distribution of s marginalized over μ. Note that the μ and s of the actual lognormal distribution, Equation (C9), are in natural log units (ln[M_d/M_*]); in these figures, we plot them in base-10 (log₁₀[M_d/M_*]) units instead (i.e., we plot μ log₁₀e and s log₁₀e instead of μ and s), for greater intuition. We obtain four main results as following.

• 1.  
  M_{d, ν}/M_* for Taurus solar-type stars. Figure 10 shows that the lognormal for the full sample of Taurus solar-type stars has a most likely mean of μ log₁₀e = $-2.4^{+0.3}_{-0.4}$, and a most likely standard deviation of s log₁₀e = $0.7^{+0.4}_{-0.2}$ (where the ranges are the ±3σ spread in probable values around the peak of the 1D μ and s distributions). Note that β for many of the solar-type stars appears to lie in the range 0–1 (see Section 7.3.2), lower than our adopted β = 1. Hence, the mean of the true disk-to-stellar mass ratio M_d/M_* (modulo the absolute value of the opacity κ) may be a factor of ∼3 higher for these stars (as noted in Section 7.3.2), corresponding to μ log₁₀e ∼−2.
• 2.  
  Effects of binarity in Taurus solar-type stars. Figure 10 also shows that the 1D distribution of the mean for the Taurus single solar types deviates by nearly 2σ from that for the full Taurus sample of these stars, with the most probable mean for the singles occurring at μ log₁₀e ≈−2 (2.5 times higher than in the full sample). Since the full sample comprises both singles and multiples, this result suggests that multiplicity has some effect on the disk mass distribution. To investigate this further, we follow AW05 in dividing the Taurus solar-type binaries into close and wide systems (projected separations of <100 AU and ⩾100 AU, respectively), and compare the two sub-samples separately to the Taurus solar-type singles. To keep the analysis clean, we exclude the handful of higher-order multiples (which can have both close and wide components) and spectroscopic binaries (whose circumbinary disks may be as bright as the disks in wide binaries and around single stars; Harris et al. 2012).Figure 11 shows that the 1D distribution of the mean for the Taurus solar-type close binaries deviates by ∼3σ from that for the single stars, with the most probable mean for the close systems occurring at μ log₁₀e ≈−3 (i.e., 10 times smaller than for the singles). The deviation between the wide binaries and singles is significantly smaller: Figure 12 shows that the 1D distribution of the mean for the wides is within 2σ of that for the singles, and peaks at μ log₁₀e =2.5 (i.e., three times smaller than in the singles). We note in passing that the 1D distribution of the standard deviation shows much less variation between the singles and close and wide binaries: as Figures 11 and 12 reveal, the s distributions in all three cases lie within ∼1σ of each other.These results are consistent with those of AW05 and Harris et al. (2012), who found (in samples consisting predominantly of solar-type stars, in line with the sample tested here) that the total disk mass in wide binaries is similar to (or modestly smaller than) around singles, while it is substantially smaller in close binaries. It is also worth recalling that the means cited above refer to log₁₀[M_{d, ν}/M_primary], i.e., are calculated under the assumption that all the disk material is associated with the primary. Since, in reality, the material is often distributed between both binary components (albeit possibly with disproportionately more around the primary; see Harris et al. 2012), the disk mass per individual component stellar mass in binaries is likely to deviate even more from the disk-to-stellar mass ratio in singles, than the numbers above indicate. Overall, these effects probably arise from a combination of tidal truncation of disks in binaries and the binary formation mechanism itself; we refer the reader to Harris et al. (2012) for a more detailed discussion.Our goal here is to examine how the disk-to-stellar mass ratio changes as a function of stellar mass. Ideally, therefore, we should conduct the same fine analysis—distinguishing between singles, close binaries, and wides—when comparing the Taurus solar-type population to the ρ Oph sample, and to the other stellar mass bins. However, as noted earlier, binarity surveys at all stellar masses in ρ Oph, and among VLMS/BDs in both Taurus and ρ Oph, are far less complete than among the Taurus solar types, precluding such a detailed investigation at present. Instead, we conservatively compare the full sample of Taurus solar types, regardless of binarity, to the other populations (whose binarity fraction is largely unknown), with the stipulation (as before) that in binaries and higher-order multiples our M_{d, ν}/M_* refers to the total disk mass per primary stellar mass.
• 3.  
  M_{d, ν}/M_* for different stellar mass bins.

  • (a)  
    The ρ Oph solar-type sample, and the [Taurus + ρ Oph] VLMS/BD and intermediate-mass samples, are all somewhat less constrained than the Taurus solar-type population, due to a combination of fewer data points and more upper limits. Nevertheless, the μ and s of their lognormal models are all consistent with that of the Taurus solar types to within 1σ (both in the full 2D and marginalized 1D parameter spaces; Figures 13–15). We thus conclude that the current data are consistent with a constant log₁₀[M_{d, ν}/M_*] ≈−2.4 all the way from intermediate-mass stars to low-mass BDs. Note that including the flux calibration uncertainties would only broaden each of the probability distributions, increasing the overlap with the Taurus solar types. Whether the binary and single star populations in these different mass bins are also individually similar to the corresponding populations among the Taurus solar types remains to be clarified in the future, with better binary statistics.
• 4.  
  Validity of a lognormal distribution. We have so far assumed that the underlying distribution of M_d/M_* for any stellar mass bin is a lognormal. Since M_d/M_* starts out relatively large (with the disk making up most of the star+disk system in the earliest evolutionary phases) and evolves eventually to zero (as the material either accretes on the star, or forms planetary bodies, or is ejected from the system), it is reasonable to suppose that in roughly coeval objects known to have disks but be more evolved than the earliest (Class 0/I) phases, M_d/M_* will be clustered about some mean, with fewer systems remaining very far above or having evolved very far below this value. The most natural such distribution is the lognormal, hence our choice.²⁴Nevertheless, our analysis thus far does not provide evidence that the distribution is a lognormal, but only supplies the most probable model parameters if it is one. There are two ways to investigate the validity of the lognormal assumption itself. In the Bayesian context, one would undertake a model comparison to derive the relative merits of various possibilities (e.g., lognormal versus power law versus exponential²⁵). We defer such rigorous analysis to the future, when more data become available. The other option is to compare the best-fit lognormals we have derived to the results of a survival analysis, where the latter yields an estimate of the underlying distribution of the data independent of any particular model, under the assumptions that the upper limits are randomly distributed and measurement errors are negligible. The general invalidity of these assumptions is what led us to a Bayesian inference in the first place. The comparison is still useful, however, both as a consistency check when upper limits do not dominate the data and measurement noise is relatively small, and as a means of underlining the limitations of survival analysis when this is not true.Figure 16 illustrates the comparison between our Bayesian results and the outcome of a survival analysis based on the KM estimator²⁶ (Feigelson & Nelson 1985; AW05; AW07a). For each of the four stellar mass bins we have considered, we plot the cumulative distribution corresponding to the lognormal with the most probable mean and standard deviation inferred from the Bayesian 1D marginalizations (see Figures 13–15). This is compared to the underlying cumulative distribution predicted by the Kaplan–Meier (KM) estimator. Note that the error bars on the latter only reflect the uncertainty in binning the (assumed random) upper limits in the survival analysis, and do not include any actual measurement noise in the data (which, as mentioned before, cannot be treated in survival analysis).We see that, for the Taurus solar-type stars and the Taurus + ρ Oph intermediate-mass stars—the two populations with the lowest fraction of upper limits and the least noise—our lognormals are in excellent agreement with the KM estimator. For the ρ Oph solar-type stars—comprising more upper limits, but still dominated by detections—our lognormal deviates slightly more from the KM estimator, but is still in good overall agreement with it. Finally, for the Taurus + ρ Oph VLMS/BDs, in which upper limits outnumber detections and the noise is greatest, our most-probable lognormal departs significantly from the survival analysis prediction. These results argue that (1) a lognormal distribution is indeed valid for the solar-type and intermediate-mass stars; and (2) by Occam's razor, it is also appropriate for VLMS/BDs, with the mismatch here between our lognormal and the KM estimator arising from the breakdown of the assumptions underlying survival analysis itself. Rigorous evidence for the validity of a lognormal for the VLMS/BDs must await more detections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Payne & Lodato (2007) have carried out simulations of planet formation via core accretion around a range of stellar/substellar masses. They find that low-mass stars with M_* ∼ 0.3 M_☉ and a mean M_d/M_* ∼ 5% can form terrestrial mass planets efficiently, as well as intermediate-mass and giant planets (though recovering the observed semimajor axis distribution of the latter two classes requires some fine tuning of Type I migration). 0.05 M_☉ BDs with the same mean M_d/M_* (translating to disk masses of a few Jupiters), on the other hand, form terrestrial mass planets less frequently, with M_planet ≳ 0.3 M_⊕ arising in ∼10% of cases, while giant planets are completely absent in this domain (timescales for collisionally driven planetesimal growth are longer around such low masses, and the disk gas dissipates before rocky cores become large enough to initiate runaway growth).

Given that the upper envelope of our observed distribution of M_{d, ν}/M_* for VLMS/BDs at ∼1 Myr lies at about 1%–2%, and that our most probable mean M_{d, ν}/M_* for these objects is only about ∼0.5%, it would appear that terrestrial planets may be hard to form around VLMS, and be extremely rare around BDs (for 0.05 M_☉ BDs with mean M_d a fraction of a Jupiter mass; Payne & Lodato find M_planet ≳ 0.3 M_⊕ in only 0.035% of cases). However, for the few VLMS/BDs with measured continuum slopes α, we showed that M_{d, ν} may underestimate the true opacity-normalized disk-mass $M_d^{\kappa}$ by a factor of ≳3–5; in general, underestimates $M^\kappa_d$ by a factor of at least a few in VLMS/BDs. M_{d, ν}, if their average continuum slope is comparable to that in more massive stars, 〈α_{850 − 1.3}〉 ∼ 2 (Section 7.3.1). Thus, disk masses of a few percent of the (sub)stellar mass may in fact be more appropriate than the naive estimate M_{d, ν}. Moreover, Payne & Lodato's simulations refer to the initial M_d (at ∼10⁴ yr); the values we observe at ∼1 Myr must be somewhat lower due to subsequent accretion/outflow. Putting the two effects together, an initial mean M_d/M_* ∼ 5% is plausible, implying that terrestrial mass planets may indeed form copiously around VLMS, and albeit less frequently, around BDs as well. This must be tested through future multi-wavelength observations that constrain the spectral slope and degree of grain growth in VLMS/BD disks.

Our inferred disk masses have implications for accretion too. Empirically, the accretion rate $\dot{M}$ onto the central star/BD seems to falloff with stellar mass roughly as $\dot{M} \propto M_{\ast } ^2$ (though with considerable scatter), all the way from intermediate-mass stars to very low mass BDs (Mohanty et al. 2005; Muzerolle et al. 2005; Natta et al. 2006). Padoan et al. (2005) proposed that this is linked to Bondi–Hoyle accretion onto the disk from the surrounding molecular cloud as the star+disk system moves through the cloud. However, there are severe problems with this hypothesis (see discussion by Mohanty et al. 2005; Hartmann et al. 2006), and it is unlikely to be the main reason for the observed relationship. Dullemond et al. (2006) and Vorobyov & Basu (2008, 2009) have presented two alternate theories. According to the former authors, the dependence of $\dot{M}$ on M_* ultimately results from the distribution of rotation rates of the parent molecular cloud cores; according to the latter, it derives from gravitational instability-driven accretion. In both cases, $\dot{M} \propto M_d$; insofar as the $\dot{M} \propto M_{\ast } ^2$ trend is real, therefore, both require $M_d \propto M_{\ast } ^2$. Our results, on the other hand, indicate a significantly shallower trend, M_{d, ν}∝M_*, all the way from intermediate-mass stars to very low mass BDs, and therefore do not support either of these two theories. Note that we have assumed a constant characteristic disk temperature $\tilde{T} = 20$ K independent of stellar mass, for estimating M_{d, ν}; this is likely to be systematically slightly in error, with the appropriate temperature being somewhat higher for the hot stars and lower for BDs, as discussed in Sections 7.3.2 and 9. Correcting for this, however, leads to disks systematically less massive with increasing stellar mass than our naive M_{d, ν} estimate, making the true disk mass an even shallower function of stellar mass than we find, which is opposite to the required effect. As such, theories requiring $M_d \propto M_{\ast } ^2$ to explain the observed $\dot{M}$ appear unsupported by data.

Alternatively, Alexander & Armitage (2006) show that the $\dot{M} \propto M_{\ast } ^2$ relationship can be reproduced even if M_d∝M_* (as our M_{d, ν} results suggest), if the initial disk radius increases with decreasing (sub)stellar mass. They note that viscous spreading means that stellar disks should rapidly expand to BD disk sizes, and moreover that external factors such as binary truncation can strongly influence the disk size, so that the inverse scaling of the initial disk radius with M_* may not be easily observable. However, their model does predict that few if any objects should have $\dot{M} < 10^{-12}$ M_☉ yr⁻¹, and that the observed scatter in accretion rates (assumed to be age-related) should be smaller in BDs than in higher-mass stars. Sensitive UV measurements by Herczeg et al. (2009) have revealed $\dot{M} \sim 10^{-12}$–10⁻¹³ M_☉ yr⁻¹ in three BDs so far, nominally at odds with this hypothesis, but more data are required to clarify whether such low rates are indeed common in BDs or a rarity as this scenario predicts.

However, the trend $\dot{M} \propto M_{\ast } ^2$ itself remains open to question (e.g., Clarke & Pringle 2006); for instance, it is unclear whether it applies to the average $\dot{M}$ at any given mass or only to the upper envelope, and also whether the exponent is indeed 2 everywhere or smaller (e.g., Vorobyov & Basu 2009 argue that $\dot{M} \propto M_{\ast } ^{1.3}$ among solar-type stars). A more careful consideration of the selection effects and upper limits in the data (as suggested by Clarke & Pringle 2006), as well as more sensitive $\dot{M}$ observations, are required to resolve this question.

We have obtained SCUBA-2 850 μm data for seven accreting VLMS and BDs in Taurus and the TWA, and combined our observations with other recent sub-mm/mm surveys of Taurus, ρ Ophiuchus, and the TWA to investigate the trends in disk mass and grain growth during the cTT phase. Assuming a disk gas-to-dust ratio of 100:1 and fiducial surface density and temperature profiles guided by current observations (Σ∝r⁻¹ and T∝r^−0.6), and using generalized equations for the disk emission, we find that:

• 1.  
  The RJ approximation is less than ideal for disks around solar-type and intermediate-mass stars, and even worse VLMS and BD disks.
• 2.  
  The minimum disk outer radius required to explain the upper envelope of sub-mm/mm fluxes is ∼100 AU for intermediate-mass stars, solar types, and VLMS, and ∼20 AU for BDs.
• 3.  
  There is a marked flattening/decrease in the upper envelope of observed fluxes going from solar-type to intermediate-mass stars, given roughly by $F_{\nu } \propto M_{\ast } ^{-1/2}$; this may be due to greater photoevaporation with increasing stellar mass. Grain growth and accretion, while present, appear unlikely to cause this trend.
• 4.  
  Many of the (likely) optically thin disks around Taurus and ρ Oph intermediate-mass and solar-type stars evince an opacity power-law index of β ∼ 0–1 (as found in previous studies as well), suggesting substantial grain growth. The situation for the VLMS/BDs in these regions is not clear: most have been observed at only one wavelength, and the few with a small measured α are consistent with either grain growth or optically thick disks.
• 5.  
  We have observed the TWA VLMS TWA30A and B only at 850 μm, but a comparison of their derived apparent disk masses and their observed accretion rates and age suggests substantial grain growth, similar to that in the TWA accretors Hen 3-600A and TW Hya.
• 6.  
  For the TWA BD 2M1207A, an analogous analysis suggests M_d/M_* ∼ 10⁻³, comparable to the lower end of measured values among Taurus and ρ Oph VLMS/BDs; the disk mass for the similar TWA BD SSSPM1102 appears to be a factor of 10 smaller still (of order 1 M_⊕). The available data put no strong constraints on grain growth in the two disks.
• 7.  
  The observed spread of ≳2 dex in apparent disk masses, in all stellar mass bins, may reflect a spread in initial disk masses, or in accretion and/or photoevaporation rates. There is no strong evidence at present of this spread being caused by a reduction in disk flux due to grain growth, but more data on faint disks are required to test this proposition.We have also examined the relationship between apparent disk mass and stellar mass through a Bayesian analysis (with the caveat that the trend in the real disk mass may be different, depending on grain properties and the absolute value of the opacity, which are unknown for a large fraction/most of the sources). We find the following.
• 8.  
  Among Taurus solar-type stars, the disk mass in close binary systems (projected separation <100 AU) is ∼10× smaller than around single stars, while the disk mass in wide binaries (⩾100 AU) is closer to that around singles (∼3× smaller), in line with previous studies.
• 9.  
  The apparent disk-to-stellar mass ratio is consistent with a lognormal distribution with a mean of log₁₀[M_disk/M_*] ≈ −2.4, all the way from intermediate-mass stars to VLMS/BDs, in agreement with previous qualitative suggestions that the ratio is roughly 1% throughout the stellar/substellar domain.
• 10.  
  We caution against the application of survival analysis techniques to astrophysical data sets in which upper limits are not random, but dependent on survey sensitivities.Finally, we have examined the implications of our results for planet formation around VLMS/BDs, and accretion rates as a function of stellar mass. We find that:
• 11.  
  While the apparent disk masses, M_{d, ν}, suggest that there may not be enough material around VLMS/BDs to efficiently form terrestrial mass planets, our detailed analysis suggests that these estimates might be systematically too low by factors of a few to 10. If so, terrestrial mass planets may form copiously around VLMS, and in the most massive BD disks as well.
• 12.  
  The observations do not support theories which require $M_d \propto M_{\ast } ^2$ in order to explain the apparent empirical correlation $\dot{M} \propto M_{\ast } ^2$ from intermediate-mass stars to BDs. Either the mechanism behind this relationship is different (e.g., related to a mass-dependent variation in disk size), or the relationship itself is more nuanced (either not as universal as thought, or driven by observational selection effects).

We thank the anonymous referee for a very close reading of the paper and extremely useful comments. S.M. is very grateful to the JCMT for awarding shared-risk SCUBA-2 time to this project, and to the International Summer Institute for Modeling in Astrophysics (ISIMA) for affording him the time and research environment necessary to take this work forward; he also acknowledges the funding support of the STFC Grant ST/H00307X/1. Part of the work by A.S. was funded by the Science Foundation Ireland through grant 10/RFP/AST2780. G.L. acknowledges financial support from PRIN MIUR 2010-2011, project "The Chemical and Dynamical Evolution of the Milky Way and Local Group Galaxies," prot. 2010LY5N2T.

We outline the theory of disk SEDs, initially following the exposition by B90. The disk flux density at a frequency ν, measured by an observer a distance D away, is

$$\begin{equation} F_\nu = \frac{{\rm cos}\,i}{D^2}\int _{r_0}^{R_d}B_\nu (T) (1 - e^{-\tau _\nu })\,2\pi r dr, \end{equation} \tag{ A1 }$$

where i is the inclination of the disk relative to the observer (with i = 90° for edge-on), B_ν(T) is the Planckian specific intensity emitted by a locally blackbody disk with temperature T(r) at a radial distance r from the central star, τ_ν(r) is the line-of-sight optical depth of the emitting material, and r₀ and R_d are the disk inner and outer edge radii, respectively. This formula is only valid if the source function (=B_ν(T)) is roughly constant with optical depth through the disk, and thus inapplicable for very large viewing angles (i → 90°); we assume that the disk is sufficiently far from edge-on to satisfy this constraint. The optical depth is then τ_ν(r) = κ_νΣ(r)/cos i, where Σ(r) is the surface density of the disk at r (so that the total disk mass is $M_d = \int _{r_0}^{R_d} \Sigma (r) \,2\pi r dr$), and κ_ν is the total opacity (gas+dust) of the emitting material (i.e., the actual opacity of the emitting dust grains scaled by the gas-to-dust ratio). We assume the surface density and temperature are power laws in radius, Σ(r) = Σ₀(r/r₀)^−p and T(r) = T₀(r/r₀)^−q, where Σ₀ and T₀ are the values at the disk inner edge r₀.

Since the surface density declines with increasing radius, the inner regions of the disk are more optically thick than the outer parts at any frequency, by the definition of τ_ν. For a given ν, the change from optically thick (τ_ν ≫ 1) to optically thin (τ_ν ≪ 1) conditions occurs at a radius r₁, found by setting τ_ν(r₁) ∼ 1:

$$\begin{equation} r_1 \,=\, \tau _0^{1/p} r_0 \,\approx \, \left[\frac{(2-p){\bar{\tau }}_\nu }{2}\right]^{1/p}R_d. \end{equation} \tag{ A2 }$$

Here, τ₀ is the optical depth at the inner edge r₀, and ${\bar{\tau }}_\nu$ is the average optical depth in the disk at ν, ${\bar{\tau }}_\nu \equiv \kappa _\nu \,M_d/(\pi R_d^2\,{\rm cos}\,i)$. The latter uses the definition of M_d in terms of the integral of Σ(r), and assumes that R_d ≫ r₀ (i.e., we have a disk, not a narrow annulus).

Lastly, for a disk with some global minimum temperature T_min, the flux densities emitted from all radii lie in the RJ limit for all frequencies ν ≪ kT_min/h. If the observed frequencies are in this regime, then the Planck spectrum reduces to B_ν(T)_RJ ≈ (2ν²/c²)kT.

Expressing Equation (A1) as the sum of two integrals representing the flux densities from optically thick radii (r₀ ⩽ r ⩽ r₁) and optically thin ones (r₁ ⩽ r ⩽ R_d), noting that the attenuation factor $(1-e^{-\tau _\nu })$ becomes ∼1 for τ_ν ≫ 1 and ∼τ_ν for τ_ν ≪ 1, and stipulating that we are in the RJ regime, finally yields

$$\begin{eqnarray} F_\nu & \approx & \nu ^2 \left(\frac{4\pi k}{c^2}\right) f_0 \left(\frac{{\rm cos}\,i}{D^2}\right) \left[T_0 r_0^q \frac{(2-p){\bar{\tau }}_\nu }{2}\,R_d^p\right]\nonumber\\ &&\times\left[\frac{R_d^{2-p-q}-r_1^{2-p-q}}{2-p-q}\right] (1+\Delta), \end{eqnarray} \tag{ A3 }$$

where the dimensionless factor f₀ corrects for our having neglected the curvature of $(1-e^{-\tau _\nu })$ over the transition from optically thick to thin regimes; f₀ ∼ 0.8 yields good agreement between Equation (A3) and numerical integration of Equation (A1) in the RJ limit (B90). The product of all the terms outside the last parentheses is the optically thin contribution to the flux density, while the term Δ is the ratio of the optically thick to thin contributions:

$$\begin{equation} \Delta \equiv \left[\frac{2}{(2-p){\bar{\tau }}_\nu \,R_d^p}\right] \left(\frac{2-p-q}{2-q}\right) \left(\frac{r_1^{2-q} - r_0^{2-q}}{R_d^{2-p-q} - r_1^{2-p-q}}\right). \end{equation} \tag{ A4 }$$

Note that B90 further simplify the product of the two square brackets in Equation (A3); their expression implicitly assumes that the r₁ derived from Equation (A2) satisfies r₀ ⩽ r₁ ⩽ R_d. If, however, Equation (A2) formally yields r₁ < r₀, then the entire disk is optically thin, and we must set r₁ ≡ r₀ for a physical solution to Equation (A3); similarly, if formally r₁ > R_d, then all radii are optically thick, and we must set r₁ ≡ R_d. In either case, B90's contraction does not hold, but our explicit Equation (A3) remains valid.

If the opacity is moreover a power law in frequency, κ_ν = κ_f(ν/ν_f)^β, where κ_f is the opacity at some fiducial frequency ν_f, and the spectral index is defined as α ≡ d(lnF_ν)/d(ln ν), then Equation (A3) can be manipulated to find α ≈ 2 + (β/1 + Δ) in the RJ limit.

We now depart from B90's discussion. If the RJ approximation is not valid, then Equations (A3) and (A4) fail. General expressions for F_ν and Δ can still be derived by dividing the disk into optically thick and thin parts, as above, but now using the general Planck formula instead of its RJ limit. The results are cited in Equations (5) and (6) in Section 4.

We assume i ≈ 60°, corresponding to the mean value of cos(i) for a random distribution of orientations (note that for optically thin disks, F_ν becomes independent of i). We also adopt r₀ ≈ 5 R_*, the usual value for a magnetospherically truncated inner edge (the precise location of the hot inner edge has negligible impact on the long wavelength emission). The location of the outer edge R_d is harder to constrain; we adopt a fiducial value of R_d = 100 AU, but also examine a large range of R_d = 10–300 AU, consistent with most sources in resolved observations of stars (e.g., AW07a; Isella et al. 2009; Andrews et al. 2010), and with the few available constraints for BDs (Scholz et al. 2006; Luhman et al. 2007).

For the minimum mass solar nebula, Hayashi (1981) derived a surface density exponent of p = 1.5. More recent spatially resolved observations of disks all imply shallower radial profiles, mostly in the range 0.5–<1.5 (e.g., Wilner et al. 2000; Isella et al. 2009; Andrews & Williams 2007b; Andrews et al. 2009, 2010). Specifically, Andrews et al. (2010) find that the surface densities are consistent with a uniform exponent of ∼1 independent of stellar mass, for solar-type to intermediate-mass stars (∼0.3–4 M_☉). We thus adopt a fiducial p = 1, and discuss the effects of varying this over the range 0.5–1.5 at appropriate junctures.

Finally, both the normalization of the temperature (T₀ here) and its radial exponent q can be estimated from SED fits to optically thick emission in the mid- to far-infrared (see B90). AW05 do this for sources in their sample with infrared data, to find median values of q = 0.58 and T(1 AU) = 148 K (they choose to normalize at 1 AU instead of at the inner edge r₀). This q is intermediate between that for flat disks, q = 3/4, and fully flared ones, q = 3/7 (Chiang & Goldreich 1997, hereafter CG97), suggesting some degree of dust settling. While the majority of stars AW05 examine in this regard are solar types, a significant fraction of the small number of intermediate-mass stars they analyze also evince intermediate q. Similarly, VLMS/BD disks exhibit mid-infrared signatures of disk flattening as well (Apai et al. 2005; Scholz et al. 2006). We thus adopt AW05's median q = 0.58 as a reasonable fiducial value.

Furthermore, we scale AW05's temperature normalization, applicable mainly to solar-type stars, to the wide range of stellar/BD masses in our sample as follows. We assume that the disks are passive, i.e., heated predominantly by stellar irradiation instead of by accretion (a good assumption for average stellar/BD accretion rates in the Class II phase). The general expression for the disk effective temperature at a radius r is then (CG97), T(r) = (α/2)^1/4(R_*/r)^1/2 T_*, where α is the grazing angle at which stellar photons impinge on the disk at r. Assuming now that α is roughly constant at the disk inner edge for all stellar masses (justified below), and recalling that we adopt r₀/R_* ≈ 5 in all cases, then yields T₀∝T_*, with the constant of proportionality independent of stellar parameters. Moreover, for solar-type stars with fiducial values T_* ∼ 4000 K and R_* ∼ 2 R_☉ (Section 3), AW05's median T(1 AU) and q (derived mainly from solar types) imply T₀(solar  types) = 148 (5 R_*/1 AU)^−0.58 K ≈ 880 K. The disk temperature normalization at r₀ for any stellar temperature T_* is then

$$\begin{equation} T_0 \approx 880 \,(T_{\ast } /4000\,{\rm K})\,{\rm K}. \end{equation} \tag{ B1 }$$

In fact, the T₀ derived this way turn out to be very close to those expected for a flat disk at r₀ (T₀ ≈ 0.2^1/4(R_*/r₀)^3/4 T_*; CG97). In the latter case, α ≈ 0.4 R_*/r₀, which is constant for all stars with our choice of r₀. Since passive disks should indeed be nearly flat very close to the star (where the flaring term in α is negligible; CG97), we see that our assumption of a constant α at our chosen r₀, and hence adoption of Equation (B1), is self-consistent.

Note that for the range of T_* covering most of our sample, ∼2500–5500 K, the T₀ we derive are lower than the dust sublimation temperature ∼1500 K; our adoption of r₀ = 5 R_* (set by magnetospheric truncation) as the inner edge of the dust disk is thus formally self-consistent (in the sense that dust can exist at these temperatures). In reality, the disk may be significantly hotter here due to accretion and/or frontal heating of the disk inner wall; the real inner edge of the dust disk may then be set by sublimation instead, and lie at somewhat larger radii (e.g., Eisner et al. 2005). However, the hot emission from regions close to r₀ contributes negligibly to F_ν in the sub-mm/mm, so the precise temperature and location of r₀ is not crucial in itself to our calculations. The true importance of T₀ for us is as a scaling factor for the temperatures at larger (cooler) radii, which do contribute to the sub-mm/mm emission. Heating at these radii is expected to be dominated by stellar irradiation instead of accretion (CG97; D'Alessio et al. 1999), and not subject to frontal irradiation effects either, so the temperatures here can be self-consistently derived by scaling under the assumption of a passive disk throughout (tapering to a small grazing angle at r₀), as we have done.

Lastly, the disk temperature cannot continue decreasing as a power law indefinitely; eventually, other heating sources, such as cosmic rays, radionuclides, and the ISRF must dominate over the stellar irradiation. Cosmic rays and radionuclides do not appear important in this context (e.g., for a disk around a solar-type star, heating from these two sources becomes important only when the temperature due to stellar irradiation has formally dropped to <3 K, i.e., below the cosmic microwave background temperature; see D'Alessio et al. 1998). The ISRF, however, is critical. Specifically, Mathis et al. (1983) calculate the temperature of grains in both the ISM and in Giant Molecular Clouds (GMCs) due to ISRF heating. Given the very high visual extinction to the midplane of accretion disks, the GMC case is more germane to our analysis. For their calculated ISRF near the galactocentric distance of the Sun, they find a grain temperature of ∼8–15 K, for silicate and graphite grains, respectively. We therefore adopt a minimum disk temperature of 10 K due to the ISRF. Note that the gas does not influence this grain temperature at all; when the density is sufficiently high (as it is in accretion disks), the gas simply reaches thermal equilibrium with the dust at the same temperature.

Thus, our final adopted temperature profile is one that declines as a power law in radius, due to stellar irradiation, until the 10 K threshold is reached, and then remains constant at this temperature, due to ISRF heating, at all larger radii.

Consider a specified model $\mathcal {M}_i$ invoked to explain a set of data $\mathcal {D}$. Then, the fundamental equation of Bayesian statistics states that

$$\begin{equation} P({\mathcal {M}}_i|{\mathcal {D}}) = \frac{P({\mathcal {M}}_i)\,P({\mathcal {D}}|{\mathcal {M}}_i)}{P({\mathcal {D}})}, \end{equation} \tag{ C1 }$$

where $P({\mathcal {M}}_i|{\mathcal {D}})$ is the posterior probability that the model is correct given the data (posterior for short); $P({\mathcal {M}}_i)$ is the prior probability that we assign to the model's veracity (before comparing to the data), given all our pre-existing information/biases about the world (prior for short); $P({\mathcal {D}}|{\mathcal {M}}_i)$ is the probability of obtaining the data given this model (known as the likelihood); and $P({\mathcal {D}}) \equiv \sum _i{P({\mathcal {M}}_i)\,P({\mathcal {D}}|{\mathcal {M}}_i)}$—where the summation is over all possible models—is a normalization factor (known as the evidence) that ensures $\sum _i{P({\mathcal {M}}_i|{\mathcal {D}})} = 1$.

Equation (C1) elegantly expresses the intuitive notion that the probability of a given model being correct, in light of the data, is proportional to both our a priori preference for the model and the likelihood of actually obtaining the observed data if the model were true. Thus, a model which fits the data well but is judged outlandish for whatever other reason must still be considered improbable, and so must a model which seems reasonable to start with but fits the data very poorly. The normalization factor simply expresses the condition that the data must be explicable by some subset of the available models.

This technique can be used for either model comparison (evaluating the posterior probabilities of different models) or parameter estimation (evaluating the posteriors for different parameter values within the context of a single model). In this work, we are concerned with the latter; we will consider one model (in particular a lognormal, discussed further below), specified by N parameters. Our task is to estimate the probability of these parameters taking on any particular set of values {θ_{n = 1..N}} (in our case, a specific mean and standard deviation), given a set of actual measurements, upper limits, and uncertainties. Thus, consider data with N_d measured values {${\hat{m}}_{d=1..N_d}$} and associated errors {σ_d}, and N_u upper limits {${\hat{m}}_{{\rm lim},u=1..N_u}$} with errors {σ_u}. Then, from Equation (C1), the posterior probability of the model parameters taking on the values {θ_n} may be explicitly written as

$$\begin{equation} P(\lbrace \theta _n\rbrace \,|\, \lbrace {\hat{m}}_d \rbrace ,\lbrace {\hat{m}}_{{\rm lim},u} \rbrace) = \frac{P(\lbrace \theta _n\rbrace)\,\, P(\lbrace {\hat{m}}_d \rbrace ,\lbrace {\hat{m}}_{{\rm lim},u} \rbrace \,|\,\lbrace \theta _n\rbrace)}{\int {P(\lbrace \theta _n\rbrace) P(\lbrace {\hat{m}}_d \rbrace ,\lbrace {\hat{m}}_{{\rm lim},u} \rbrace \,|\, \lbrace \theta _n\rbrace)}}, \end{equation} \tag{ C2 }$$

where the evidence in the denominator is expressed as an integral, instead of a summation, for a continuous instead of discrete range for each model parameter.

The evidence integral is immaterial for inferring the relative probabilities of model parameters. We will also assume here that the prior is a uniform distribution, expressing our lack of any strong preference for one set of model parameters over another. Then, the posterior is simply proportional to the likelihood:

$$\begin{equation} P(\lbrace \theta _n\rbrace \,|\, \lbrace {\hat{m}}_d \rbrace ,\lbrace {\hat{m}}_{{\rm lim},u} \rbrace) \propto P(\lbrace {\hat{m}}_d \rbrace ,\lbrace {\hat{m}}_{{\rm lim},u} \rbrace \,| \,\lbrace \theta _n\rbrace). \end{equation} \tag{ C3 }$$

The likelihood is given by the product of the likelihood of each measured value or upper limit, i.e.,

$$\begin{eqnarray} P(\lbrace {\hat{m}}_d \rbrace ,\lbrace {\hat{m}}_{{\rm lim},u} \rbrace \,| \,\lbrace \theta _n\rbrace) &=& \left[\prod _{d=1}^{N_d}P({\hat{m}}_d \,|\, \lbrace \theta _n\rbrace)\right]\nonumber\\ &&\times \left[\prod _{u=1}^{N_u}P({\hat{m}}_{{\rm lim},u} \,|\, \lbrace \theta _n\rbrace)\right].\quad\quad \end{eqnarray} \tag{ C4 }$$

The one subtlety in evaluating the above is that the true value of each data point (as opposed to the measured or upper limit value, which is the true value scattered by the measurement noise) must be taken into account (i.e., marginalized over), since we want the underlying distribution of the population. This is accomplished by making use of the known uncertainties (noise) as follows.

C.1. Measured Values

We will assume that the true value m_d is non-negative, and that the noise is additive and Gaussian. Then, for objects with measurements, the individual likelihoods are given by

$$\begin{eqnarray} P({\hat{m}}_d \,|\, \lbrace \theta _n\rbrace) &=& \int _0^\infty P(m_d \,|\, \lbrace \theta _n\rbrace) \,\frac{1}{(2\pi)^{1/2}{{\sigma }}_d }\nonumber\\ &&\times{\rm exp}\left[-\frac{1}{2}\left(\frac{{\hat{m}}_d - m_d}{{{\sigma }}_d }\right)^2\right] \,{d}m_d.\quad\quad \end{eqnarray} \tag{ C5 }$$

In general, the integral must be evaluated numerically once the model is specified, which is nevertheless usually straightforward since the integrand is only significant in the approximate range ${\rm MAX}(0,{\hat{m}}_d - 3{{\sigma }}_d) \lesssim m_d \lesssim ({\hat{m}}_d + 3{{\sigma }}_d)$. Note that the measured value ${\hat{m}}_d$ may well be negative; it is only the true value m_d that is required to be non-negative.

C.2. Upper Limits

For objects without reported measured values, the individual likelihood is complicated by the fact that the upper limit is consistent with any case in which the unreported measurement ${\hat{m}}_u$ (i.e., the true value m_u scattered by the measurement noise σ_u) is less than or comparable to the reported limit ${\hat{m}}_{{\rm lim},u}$. Hence, it is necessary to marginalize not only over the unknown true value m_u, but also over the unknown measurement ${\hat{m}}_u$. Assuming again that the noise is additive and Gaussian, the likelihood for each object with a reported upper limit becomes

$$\begin{eqnarray} &&P({\hat{m}}_{{\rm lim},u} \,|\, \lbrace \theta _n\rbrace) = \int _0^\infty P(m_u \,|\, \lbrace \theta _n\rbrace) \,\left\lbrace \int _{-\infty }^{{\hat{m}}_{{\rm lim},u} }\frac{1}{(2\pi)^{1/2}{{\sigma }}_{u} }\right. \nonumber\\ &&\left. \quad\quad\times{\rm exp}\left[-\frac{1}{2}\left(\frac{{\hat{m}}_u - m_u}{{{\sigma }}_{u} }\right)^2\right] {d}{\hat{m}}_u \right\rbrace \,{d}m_u. \end{eqnarray} \tag{ C6 }$$

The inner integral is given by

$$\begin{eqnarray} &&\int _{-\infty }^{{\hat{m}}_{{\rm lim},u} }\frac{1}{(2\pi)^{1/2}{{\sigma }}_{u} } {\rm exp}\left[-\frac{1}{2}\left(\frac{{\hat{m}}_u - m_u}{{{\sigma }}_{u} }\right)^2\right] {d}{\hat{m}}_u\nonumber\\ &&\quad\quad= \frac{1}{2}\,{\rm erfc}\left(\frac{{\hat{m}}_{{\rm lim},u} - m_u}{{{\sigma }}_{u} }\right), \end{eqnarray} \tag{ C7 }$$

where erfc is the complementary error function, so that the likelihood simplifies to

$$\begin{equation} P({\hat{m}}_{{\rm lim},u} |\lbrace \theta _n\rbrace) = \int _0^\infty\! P(m_u |\lbrace \theta _n\rbrace) \frac{1}{2}\,{\rm erfc}\left(\frac{{\hat{m}}_{{\rm lim},u} - m_u}{{{\sigma }}_{u} }\right) {d}m_u. \end{equation} \tag{ C8 }$$

This integral must be evaluated numerically, but this is again straightforward, with the integrand appreciable only over the range $0 \le m_u \lesssim ({\hat{m}}_{{\rm lim},u} + 3{{\sigma }}_{u})$.

C.3. Model

We assume that the underlying distribution of M_{d, ν}/M_* is a lognormal (justified at the end of Section 8), specified by two parameters: θ₁ ≡ μ (the mean) and θ₂ ≡ s (the standard deviation). For any given μ and s, the probability of any particular positive value m (≡ M_{d, ν}/M_*) in this distribution is

$$\begin{equation} P(m \,|\, \lbrace \mu , s\rbrace) \equiv \frac{1}{(2\pi)^{1/2}\,m\,s}{\rm exp}\left[-\frac{1}{2}\left(\frac{{\rm ln}(m) - \mu }{s}\right)^2\right]. \end{equation} \tag{ C9 }$$

Substituting this into Equation (C5) with m ≡ m_d, and into Equation (C8) with m ≡ m_u, then gives us the individual likelihoods for any measured value ${\hat{m}}_d$ or upper limit ${\hat{m}}_{{\rm lim},u}$. The product of these, as specified in Equation (C4), yields the total likelihood of the data given any particular parameter set {μ, s}, which leads to the posterior probability of this set via Equation (C3). Repeating this procedure over the whole range of model parameters then gives the posterior probability distribution of the model parameters (in reality, we carry this out over a fine but discrete and finite mesh covering the range of plausible μ and s). This method of calculating posterior probabilities is graphically illustrated in Figure 17.

Zoom In Zoom Out Reset image size