Hints for Small Disks around Very Low Mass Stars and Brown Dwarfs^∗

Our sample comprises 11 VLMOs and brown dwarfs from two nearby star-forming regions, Taurus and Chamaeleon I (d = 162 pc; Luhman 2008). Because we use stellar parameters from Rebull et al. (2010) for the majority of our Taurus sources, we adopt their distance of 137 pc for the Taurus star-forming region (Torres et al. 2007, 2009). Ten of the sources were selected by us and observed with the Herschel/PACS spectrograph as part of Proposal ID OT2_ipascucc_2; see Section 2.2 for more details. CIDA-1 is the only other VLMO observed with the same setting (Proposal ID OT1_ascholz_1); hence, we include it in our analysis.

We chose to focus on the Taurus and Chamaeleon I star-forming regions because the low mass end of the stellar population is well characterized and, unlike ρ Oph, suffers from minimal contamination from extended cloud emission. In addition to being VLMOs/brown dwarfs, our targets were selected to posses a dust disk (based on infrared excess emission) and to have a flux density at 24 μm (F₂₄) greater than 50 mJy (Guieu et al. 2007; Luhman et al. 2008; Rebull et al. 2010). CIDA-1 also fits these criteria. The 24 μm flux limit was applied to ensure bright 70 μm fluxes (F₇₀), greater than ∼70 mJy, based on the empirical relation ${\nu }_{70}{F}_{70}\sim 0.5\times {\nu }_{24}{F}_{24}$ from the VLMO disk sample of Harvey et al. (2012). Because of the known positive correlation between the far-infrared continuum and the [O i] 63 μm emission line (e.g., Howard et al. 2013; Keane et al. 2014), our target selection criteria were intended to maximize the detection of the [O i] forbidden line (see also Section 2.2).

Table 1 provides the main properties of our targets. Spectral types are converted to effective temperatures for each object using the relationship given in Luhman et al. (2003). Because this relationship deviates from linearity between spectral types M6 and M9, quadratic interpolation was used to estimate temperatures for non-integer spectral types. Stellar masses were computed by comparing each source position in the Hertzsprung–Russell diagram to the Baraffe et al. (1998) evolutionary tracks. While this is a standard approach and the only one available for our sample, one should keep in mind that model-inferred masses have large uncertainties of 50%–100% below $1\,{M}_{\odot }$ as demonstrated by Stassun et al. (2014) using a sample of eclipsing binary systems.

Table 1.  Source Properties

2MASS	Common	Region	SpTy	References	${A}_{{\rm{v}}}$	log ${L}_{\mathrm{bol}}$	References	${T}_{\mathrm{eff}}$	M*
	Name					(${L}_{\odot }$)		(K)	(${M}_{\odot }$)
J04141760+2806096	CIDA-1	Taurus	M5	L17	4.6	−0.7	R10	3125	0.19
J04193545+2827218	FR Tau	Taurus	M5.25	L17	0.0	−0.9	R10	3090	0.14
J04233539+2503026	FU Tau A	Taurus	M7	L17	2.0	−0.72	L09	2880	0.16
J04295950+2433078	CFHT-20	Taurus	M5	L17	4.6	−0.7	R10	3125	0.19
J04381486+2611399	J04381486	Taurus	M7.25	L17	3.5	−2.3	R10	2846	0.05
J04382134+2609137	GM Tau	Taurus	M5	L17	2.0	−1.15a	H08	3125	0.14
J04393364+2359212	J04393364	Taurus	M5	L17	1.3	−1.0	R10	3125	0.18
J04394748+2601407	CFHT-4	Taurus	M7	L17	4.4	−1.0	R10	2880	0.11
J11062554–7633418	ESO Hα 559	Cha I	M5.25	L07	3.58b	−1.28	L07	3088	0.14
J11071668–7735532	Cha Hα 1	Cha I	M7.75	L07	0.0b	−1.82	L07	2765	0.04
J11105597–7645325	Hn 13	Cha I	M5.75	L07	2.91b	−0.89	L07	3021	0.14

Notes. A distance of 137 and 162 pc is assumed for the Taurus (R10) and Chamaeleon I sources (L07).

^aGM Tau luminosity from D14. ^b L07 report the extinction at J band. From those values we computed the ${A}_{{\rm{V}}}$ in this table using the interstellar extinction law in Cardelli et al. (1989) and an R_v of 3.1.

References. D14: Davies et al. 2014; H08: Herczeg & Hillenbrand 2008; L07: Luhman et al. 2007; L09: Luhman et al. 2009; L17: Luhman et al. 2017; R10: Rebull et al. 2010.

Download table as:  ASCIITypeset image

Two of our sources, FU Tau and Hn 13, are multiple. FU Tau has an M9.25 companion at a separation of $5\buildrel{\prime\prime}\over{.} 7$ and may have an additional spectroscopic binary companion (Luhman et al. 2009). Hn 13 has a $0\buildrel{\prime\prime}\over{.} 13$ near-equal mass companion (imaged by Ahmic et al. 2007).

Our sample was observed in 2013 February and March with the ESA Herschel Space Observatory (Pilbratt et al. 2010) Spectrograph (PACS; Poglitsch et al. 2010), which provides an integral field unit (IFU) with a 5 × 5 array of spaxels, each with a pixel size of $9\buildrel{\prime\prime}\over{.} 4$. The sources were observed under two different programs (OT2_ipascucc_2 10 sources and OT1_ascholz_1 for CIDA-1). Both programs used the standard point-source spectroscopy line scan mode with chop/nod, the chopper throw, set to small and the standard faint line mode selected (high grating sampling density). Spectra were centered on the [O i] 63 μm line and covered from 62.93 to 63.43 μm. The red channel observations were a bonus and covered from 188.77 to 190.35 μm.

We set the exposure times for the 10 OT2_ipascucc_2 sources by extrapolating the 63 μm continuum-line relationship from Howard et al. (2013) down to a minimum flux density of 70 mJy (see Section 2.1) and found a corresponding [O i] 63 μm line flux of 4 × 10⁻¹⁸ W m⁻². The exposure time was set to 16,420 s to achieve a 3σ detection on this line, i.e., a sensitivity of 1.4 × 10⁻¹⁸ W m⁻², about 3 times better than that achieved for TT stars with the GASPS survey (Dent et al. 2013). CIDA-1, the brightest source in our sample and observed under OT1_ascholz_1, has a shorter exposure time of 9868 s. The particulars of the observations, including the Herschel observation identification numbers (ObsIDs), are listed in Table 2.

The data were reduced from level 0 using the Herschel Interactive Processing Environment (HIPE; Ott 2010) version 14.0.1 with calibration set product 72. Because the targets are all faint point sources, the "ChopNodBackgroundNormalization" (version 1.38.4.3) pipeline script was used for reduction. This pipeline flags bad data (e.g., mechanism movements, saturation, overly noisy or bad pixels), applies signal corrections (e.g., nonlinearities, cross-talk, transients, capacitance ratios), performs wavelength and flux calibration, and corrects data for spacecraft velocity. For the final cube rebinning, the wavelength grid was sampled using oversample = 2 and upsample = 4. Proper flat-fielding was verified interactively from within HIPE for each target.

Because telescope jitter and pointing errors can result in flux beyond the central spaxel, the appropriate spectra output from the task extractCentralSpectrum must be chosen as described in the PACS Data Reduction Guide: Spectroscopy (Section 8.4.1). The central spaxel with the application of the c1-to-total point-source correction was used to extract the final spectrum. To further test this choice, the spectra of neighboring spaxels were examined for extended emission using the "comboplot" output. All neighboring spaxels were confirmed to have no signal above the rms. For all of our sources the central spaxel (c1) spectra had the best signal-to-noise ratio (S/N) and a signal that was comparable to those extracted from the 3 × 3 corrected spectra. The absolute flux calibration uncertainty for the PACS spectrograph at ∼63 μm is estimated to be 10% (PACS Observer Manual; HERSCHEL-HSC-DOC-0832, Version 2.5.1). The final spectra for all of our targets are shown in Figure 1.

Zoom In Zoom Out Reset image size

To identify [O i] 63 μm detections, we smooth the spectrum using a uniform filter (width of three resolution elements) before fitting each spectrum within ±0.1 μm of the line using a Levenberg–Marquardt algorithm assuming a Gaussian for the line profile and a first-order polynomial for the continuum. The 1σ uncertainties on the line fluxes are evaluated from the standard deviation of the pixels in the spectrum minus the best-fit model. We consider a line to be detected when its flux is greater than 3 times the 1σ uncertainty. In case of nondetections, we fit the same spectral range with a first-order polynomial, and we provide in Table 3 the 3σ upper limits from the rms in the baseline-subtracted spectrum using a Gaussian with a width equal to that of an unresolved line (FWHM of 98 km s^–1 at 63.18 μm). Sources are considered detected in the continuum if the S/N at 63.2 μm is greater than 3.

Following this approach, we detect the [O i] 63 μm line of two sources in our sample (toward FU Tau A and Hn 13) and the continuum of three bright disks (CIDA-1, J04393364, and ESO Hα 559). Two of these bright disks (CIDA-1 and ESO Hα 559) are also detected in the continuum at 189 μm. FU Tau A has a known molecular outflow (Monin et al. 2013), and it is possible that the outflow contributes to the [O i] 63 μm line, in which case the flux we report in Table 3 would be an upper limit for the disk emission.

Riviere-Marichalar et al. (2016) analyzed 362 Herschel sources, which included our 11 targets, using HIPE 12.0. Our results are consistent with Riviere-Marichalar et al. (2016) except for FU Tau A, where we report an $\sim 5\sigma$ [O i] 63 μm detection ($5.17\pm 1.04\times {10}^{-18}\,{\rm{W}}\,{{\rm{m}}}^{-2}$), while Riviere-Marichalar et al. (2016) report a $3\sigma$ upper limit of ($\lt 3\times {10}^{-18}\,{\rm{W}}\,{{\rm{m}}}^{-2}$).

All of the eight Taurus sources were also observed with the PACS photometer by Bulger et al. (2014). Literature 70 μm flux densities are reported in the last column of Table 3. Two of the sources (FR Tau and GM Tau) are fainter than the 70 mJy flux density estimated from the 24 μm photometry available at the time of the proposal submission. Except for CIDA-1, our continuum values and upper limits are consistent with the literature values within 1σ of the uncertainties we quote. In the case of CIDA-1 the 63 and 70 μm flux densities are within 2σ. It is possible that source variability further contributes to the flux discrepancy for this source, as 70 μm variability can be close to ∼20% even for non-embedded, young, low mass stars (Class II SEDs; see Billot et al. 2012). Given the better sensitivity of the PACS photometer and smaller uncertainty of the 70 μm photometric values, we will use them, when available, in the analysis that follows.

We combine continuum observations at multiple wavelengths and carry out continuum RT modeling to constrain some of the main disk properties, with a focus on the outer disk radius, as it might affect the [O i] 63 μm emission (e.g., Kamp et al. 2011). Each source's SED is shown in Figure 2, and individual fluxes and references are given in Table 4 and are available online.

Zoom In Zoom Out Reset image size

We model the SEDs with the 3D axisymmetric Monte Carlo RT code MCMax (Min et al. 2009). MCMax self-consistently calculates the disk vertical structure for a given dust surface density profile and gas-to-dust ratio, assuming hydrostatic equilibrium and ${T}_{\mathrm{gas}}={T}_{\mathrm{dust}}$. The vertical structure of the dust is calculated from an equilibrium between dust settling and vertical mixing as described in Dullemond & Dominik (2004). Calculation of the disk temperature and density structure are iterated until a self-consistent solution is reached. For a given grain size distribution, the only free parameter controlling the vertical structure is the turbulent mixing strength, which was not found to be different in modeling the median SEDs of Herbig stars, TT stars, and brown dwarfs by Mulders & Dominik (2012).

We begin by fitting the stellar parameters and disk inclinations of our sources using a genetic algorithm. The stellar luminosity and temperature are taken from Table 1 and allowed to vary only within their typical uncertainties, 30% for the luminosity and ±100 K for the temperature (Luhman et al. 2007). Higher luminosities were required for CFHT-4 (2 times higher) and J04381486 (10 times higher) in order to achieve good fits. This is likely due to the high inclination of these disks. Ricci et al. (2014) estimate a disk inclination of $77^\circ$ from a resolved millimeter continuum image of CFHT-4, and Luhman et al. (2007) estimate an inclination of 67°–71° for J04381486 using a Hubble Space Telescope scattered light observation. GM Tau required a ${T}_{\mathrm{eff}}$ that was 189 K lower than the literature value in order to fit the optical/near-IR portion of the SED. These quantities are used to compute stellar radii from the Stefan–Boltzmann relation. The source distance is fixed to that reported in Table 1. The literature extinction in Table 1 provides a good fit for all sources. Eight equally spaced (in cosine) inclinations are sampled to constrain the viewing angle of the disk.

The general setup for all modeling follows the one described in Mulders & Dominik (2012) to fit the median SEDs of TT stars and brown dwarfs. Given the rather coarse SED sampling at long wavelengths and degeneracies in SED modeling (e.g., Woitke et al. 2016), we limit the number of free parameters. The surface density is a power law of the form ${\rm{\Sigma }}\propto {r}^{-1}$, while the minimum and maximum grain sizes are $a=0.1\,\mu {\rm{m}}$ and 1 mm with a particle size distribution proportional to ${a}^{-3.5}$, following Mulders & Dominik (2012). With this approach, we can explore the impact of the outer radius on the SED under the explicit assumption that the surface density slope and grain size distribution are stellar mass independent.

To estimate confidence intervals, we run a grid of 4400 models for each source that explores the disk structure (inner and outer radius) and the disk mass. For the disk dust mass, 20 logarithmically spaced steps between 10⁻² and ${10}^{-9}{M}_{\odot }$ are sampled. The disk inner radius is explored using 11 logarithmically spaced steps between 10⁻⁴ and 10¹ au. For outer radius, 20 logarithmically spaced steps from 1 to 250 au are used. These ranges are chosen such that there is a peak in the Bayesian probability distribution. For sources where a clear peak is not observed, physically unrealistic values justify the boundaries of our grids (e.g., dust masses >25% stellar mass, ${R}_{\mathrm{in}}$ $\ll$ silicate sublimation radius).

We also determine the likelihood of each model in order to establish confidence intervals for the parameter space of our model grids. The likelihood of each model is calculated by comparing photometric data with model flux densities using ${\rm{\exp }}(-{\chi }_{R}^{2}/2)$, where ${\chi }_{R}^{2}$ is the reduced chi-squared metric. For photometric detections, we assume a typical uncertainty of 20% for each flux measurement. Upper limits are included in our fitting by using the modified ${\chi }^{2}$ statistic described by Sawicki (2012).

Relative probabilities are then calculated for each parameter by summing the likelihood of common parameter values normalized by the sum of all likelihoods. The 68.3% confidence intervals for the parameters ${R}_{\mathrm{out}}$, ${R}_{\mathrm{in}}$, and ${M}_{\mathrm{dust}}$ are estimated from Bayesian probability distributions and reported in Table 5. Figure 3 shows the probability distribution and confidence interval for ${R}_{\mathrm{out}},$ for each source. See Appendix B for additional information on the confidence intervals and for examples of probability distributions for ${R}_{\mathrm{in}}$ and ${M}_{\mathrm{dust}}$.

Zoom In Zoom Out Reset image size

Table 5.  Maximum Likelihood Parameters for Source Models

	i	${R}_{\mathrm{in}}$		${R}_{\mathrm{out}}$		Log ${M}_{\mathrm{dust}}$		$\langle {T}_{{\rm{d}}}\rangle$
Source	(deg)	(au)		(au)		(${M}_{\odot }$)		(K)
CIDA-1	35	0.32	${}_{-0.31}^{+1.5}$	7.6	${}_{-3.9}^{+60}$	−4.9	${}_{-0.92}^{+0.55}$	48
FR Tau	77	0.32	${}_{-0.32}^{+1.5}$	3.2	${}_{-2.3}^{+25}$	−6.1	${}_{-1.3}^{+0.92}$	46
FU Tau A	56	0.00032	${}_{-0.00026}^{+0.017}$	1.3	${}_{-0.47}^{+36}$	−7.2	${}_{-0.92}^{+1.3}$	140
CFHT-20	56	0.0001	${}_{-4.4\times {10}^{-5}}^{+0.018}$	4.3	${}_{-2.7}^{+46}$	−5.7	${}_{-1.3}^{+0.92}$	64
J04381486	70	0.032	${}_{-0.031}^{+0.025}$	10	${}_{-3.6}^{+5.6}$	−4.9	${}_{-0.18}^{+0.55}$	33
GM Tau	35	0.0032	${}_{-0.0031}^{+0.015}$	1.8	${}_{-0.92}^{+3.2}$	−6.8	${}_{-0.55}^{+3.1}$	73
J04393364	35	0.01	${}_{-0.0094}^{+0.55}$	4.3	${}_{-3.4}^{+4.6}$	−5.7	${}_{-1.3}^{+2.4}$	42
CFHT-4	80	0.01	${}_{-0.0099}^{+0.0078}$	78	${}_{-66}^{+43}$	−5.3	${}_{-0.92}^{+0.18}$	15
ESO Hα 559	35	0.032	${}_{-0.03}^{+0.53}$	33	${}_{-17}^{+18}$	−4.2	${}_{-0.18}^{+1.7}$	18
Cha Hα 1	35	0.1	${}_{-0.099}^{+0.46}$	33	${}_{-24}^{+88}$	−3.1	${}_{-0.55}^{+1.3}$	16
Hn 13	35	0.01	${}_{-0.0099}^{+0.0078}$	3.2	${}_{-2.3}^{+13}$	−6.1	${}_{-1.3}^{+0.92}$	64

Note. The maximum likelihood values are shown for the parameters ${R}_{\mathrm{in}},\ {R}_{\mathrm{out}},\ {\rm{and}}\ {M}_{\mathrm{dust}}$, along with their corresponding 1σ confidence intervals. $\langle {T}_{{\rm{d}}}\rangle$ is calculated from the best-fit model selected for the source as described in Section 3.2.

Download table as:  ASCIITypeset image

Along with the resulting parameters of our best-fit models, we also report the mass-averaged disk temperature ($\langle {T}_{{\rm{d}}}\rangle$) in Table 5. Here, we take $\langle {T}_{{\rm{d}}}\rangle$ to be the average of the dust temperature at each model grid point weighted by the dust mass at that point. In order to demonstrate the relationship between this value and the disk outer radius, we choose the model with the maximum likelihood from the set of models containing the most probable outer radius value as our fiducial model with which to calculate $\langle {T}_{{\rm{d}}}\rangle$.

Because we find maximum likelihood values of ${R}_{\mathrm{out}}$ that are substantially smaller than those of disks around TT stars (∼10 times smaller on average), we also run a second grid of models keeping ${R}_{\mathrm{out}}$ fixed to 200 au, a typical value used for modeling TT disks. These models use the same 20 and 11 steps of ${M}_{\mathrm{dust}}$ and ${R}_{\mathrm{in}}$, respectively, outlined above. We use these grids to similarly derive a $\langle {T}_{{\rm{d}}}\rangle$ with which to compare our first grid of variable ${R}_{\mathrm{out}}$ against (see Figure 5 below).

Because heating of the dust in the outer disks of VLMOs by interstellar radiation may be significant, all of our models include an isotropic interstellar radiation field (approximated by a diluted 20,000 K blackbody) in combination with a cosmic microwave background temperature of 2.7 K to properly compute the temperature in the outer disk (see Woitke et al. 2009).

In addition to modeling individual sources, we carry out two sets of RT models to assess the impact of the disk outer radius (${R}_{\mathrm{out}}$) on the derived disk mass and $\langle {T}_{{\rm{d}}}\rangle$. Both parameter studies are calculated for a range of stellar masses, but in one case ${R}_{\mathrm{out}}$ is held fixed at 200 au (hereafter, fixed ${R}_{\mathrm{out}}$ study), while in the other case (hereafter, variable ${R}_{\mathrm{out}}$ study) ${R}_{\mathrm{out}}$ scales linearly with stellar mass (see also Table 6). For both studies, the scaling of disk mass with stellar mass is chosen to reproduce the best-fit relation to the 887 μm flux densities versus stellar masses in the Chameleon I star-forming region (Pascucci et al. 2016).

Table 6.  Parameter Study Model Inputs

Parameter	Variable ${R}_{\mathrm{out}}$	Fixed ${R}_{\mathrm{out}}$
${R}_{\mathrm{out}}$(au)	$200({M}_{* }/{M}_{\odot })$	200
${M}_{\mathrm{dust}}$	$\propto {M}_{* }^{2.0}$	$\propto {M}_{* }^{1.2}$
M* (M⊙)	0.05–2.239
${R}_{\mathrm{in}}$(au)	$0.05\times {L}_{* }^{0.5}$

Download table as:  ASCIITypeset image

To cover a large range in stellar mass, we combine the evolutionary tracks from Baraffe et al. (2015) and Feiden (2016) as in Pascucci et al. (2016). To find corresponding luminosities and temperatures for each stellar mass, we use the 2 Myr isochrone because this age is a good match to the age of the star-forming regions we consider in this study. We set 15 equally log-spaced stellar masses from 0.05 to 2.239 M_⊙, which cover from 0.019 to 6.33 L_⊙.

All models adopt the same dust properties, turbulent mixing strength, and surface density power-law profiles as those to fit individual SEDs (Section 3.2). Inner radii follow a scaling relationship consistent with the temperature at which silicate sublimation would occur. All parameters are summarized in Table 6.

One important difference between the two parameter studies is the scaling of disk mass with stellar mass. The fixed ${R}_{\mathrm{out}}$ study requires a shallower disk mass–stellar mass relation to reproduce the best-fit 887 μm flux and yields a different temperature–luminosity relation. Implications of this parametric study will be discussed in Section 4.3.

As mentioned in Section 3.1, we detect the [O i] 63 μm line only toward 2 out of the 11 VLMO disks in the sample. Figure 4 shows the [O i] line flux versus the continuum of our sample (red and black symbols) in the context of literature values for more luminous/massive objects (gray symbols). The nondetections are somewhat surprising given that most sources should have been detected if they followed the line luminosity–continuum relation for TT stars with no known jets from Howard et al. (2013) (dashed line). This suggests that VLMO disks may be underluminous in the [O i] 63 μm line.

Zoom In Zoom Out Reset image size

In order to consider the possibility that the trend in the line flux–continuum relation does not extend to the VLMO regime and that the VLMO [O i] 63 μm emission is underluminous, we have refit the same data and included the TT sources with upper limits that Howard et al. omitted from their fit. By doing this, we are able to then refit the full TT sample (detections and upper limits) with our VLMO nondetections⁸ (see Figure 4).

To do this, we use the Bayesian method of regression fitting described in Kelly (2007),⁹ which takes into account nondetections and errors bars.

We exclude FU Tau A from the fitting because it is a known outflow source (Monin et al. 2013). Table 7 gives the fit coefficients.

The blue line in Figure 4 shows the best fit with the addition of our sources. VLMO disks may be underluminous in the [O i] 63 μm line, as suggested by the brightest (in continuum) half of the VLMO sample with no [O i] detection, but the upper limits are not stringent enough to conclude that they are actually underluminous at a confidence greater than 89%. On average our samples appear to be underluminous by a factor of 1.8.

Our RT modeling suggests greater likelihoods for outer disk radii that are smaller (maximum likelihoods $\leqslant 78\,\mathrm{au}$ with a mean value of 15 au) than typical values (>150 au) for disks around TT stars (e.g., Andrews et al. 2009; Isella et al. 2009; Guilloteau et al. 2011). We also note that for every source, the confidence intervals (see Figure 3) lie entirely well below 200 au. In fact, with the exception of CFHT-4, the confidence intervals lie entirely below 100 au, and the Bayesian probability distributions decrease toward larger outer radii.

Although these best-fit models may not be unique solutions and stellar mass dependencies in other disk parameters might yield equally good fits (e.g., Woitke et al. 2016), this approach shows that smaller disks around substellar objects are consistent with the SEDs.

Disk sizes are mostly unknown in the brown dwarf regime. For our sample spatially resolved millimeter images exist only for CIDA-1, CFHT-4, and ESO Hα 559 (Ricci et al. 2014; Pascucci et al. 2016). Detailed modeling has been carried out for the first two with estimated disk radii of $66{\pm }_{12}^{14}$ au for CIDA-1 and greater than 80 au for CFHT-4 (Ricci et al. 2014). ESO Hα 559 was resolved at $887\,\mu {\rm{m}}$, and a fitted elliptical Gaussian resulted in an FWHM of 023 × 015 (37 au × 24 au; Pascucci et al. 2016). The maximum likelihood outer radius values found by our modeling of CFHT-4 (78 au) and ESO Hα 559 (33 au) are in good agreement with these findings.

Ricci et al. (2014) find an outer radius for CIDA-1 that differs significantly from our value of 7.6 au; however, our modeling does not reject large sizes, and the $1\sigma$ confidence interval of 60 au is consistent with the estimate of Ricci et al. These three disks are all larger than the median VLMO disks in our sample, but they are also three of the four brightest disks. This suggests that smaller disks are indeed more common in the brown dwarf regime (see also Testi et al. 2016). Additionally, Luhman et al. (2007) found that their models of J04381486 with outer radii of 20 and 40 au agreed "reasonably well" with Hubble Wide Field Planetary Camera observations. These values are larger than our finding of ${R}_{\mathrm{out}}={10}_{-3.6}^{+5.6}$ au, in line with observations showing that submicron grains and gas disks, as probed via scattered light images, are typically larger than dust disks traced at millimeter wavelengths (e.g., de Gregorio-Monsalvo et al. 2013).

If VLMO disks are indeed smaller than those around TT stars, an important implication is that they are hotter than we would estimate if they were the same size as disks around TT stars, i.e., their $\langle {T}_{{\rm{d}}}\rangle$ is higher.

The left panel of Figure 5 illustrates this point: circles are our best-fit values color-coded by disk outer radius, while the dotted and solid lines are the $\langle {T}_{{\rm{d}}}\rangle \mbox{--}{L}_{\mathrm{bol}}$ relations by Andrews et al. (2013) and by our fixed ${R}_{\mathrm{out}}$ models (Section 3.3). In both models the dust disk radius is 200 au regardless of stellar luminosity/mass, but our fixed ${R}_{\mathrm{out}}$ models are significantly cooler because we compute the disk vertical structure self-consistently, which results in less vertically extended disks. Regardless of these differences, it is clear that our best-fit RT models point to disk radii smaller than 200 au and higher $\langle {T}_{{\rm{d}}}\rangle$ for VLMO disks than typically assumed. The SED fits with fixed ${R}_{\mathrm{out}}$ have temperatures in between those predicted by the Andrews et al. (2013) and our fixed ${R}_{\mathrm{out}}$ grid relations (right panel of Figure 5).

Zoom In Zoom Out Reset image size

Note that the Herschel models with fixed 200 au outer radius fall above (and not on) the trend line produced by our 200 au parameter models for two reasons. First, the Herschel sources were selected to be brighter, and hence hotter, disks, whereas the model grid was fit to the median disk mass in Chameleon. Second, the parameter models use disk masses and inner disk radii that scale with stellar mass and bolometric luminosity, respectively, whereas our best-fit SED models allow these values to vary for a best fit with photometry.

Interestingly, our variable ${R}_{\mathrm{out}}$ grid VLMO models are found to be optically thick at $63\,\mu {\rm{m}}$ for all of our sources with the exception of CFHT-4 (the source with the largest ${R}_{\mathrm{out}}$) and at $850\,\mu {\rm{m}}$ for all but two (CFHT-4 and FR Tau). This is different than what is proposed in Harvey et al., who, however, targeted a fainter $70\,\mu {\rm{m}}$ sample of brown dwarf disks. Measuring dust disk sizes will be important to establish if, and by how much, brown dwarf disks are optically thick at these wavelengths.

Dust disk masses are typically measured by using a sub-mm/mm data point, source distance, dust opacity, and average dust temperature as inputs to the Planck function (e.g., Beckwith et al. 1990) while assuming optically thin emission. This same approach has been recently used to estimate dust disk masses for hundreds of stars in nearby star-forming regions and associations. As discussed in Pascucci et al. (2016), such estimates, as well as the disk mass–stellar mass scaling relation, are very sensitive to the assumed $\langle {T}_{{\rm{d}}}\rangle$ and how it scales with stellar luminosity.

In the previous section we showed that the SEDs of VLMOs can be well reproduced with disks that have smaller radii than TT disks and, as such, do not fall on the temperature–luminosity relation derived by Andrews et al. (2013). Instead, they are hotter, plotting above this relation (see Figure 6 and Table 5). We use the two RT grids discussed in Section 3.3 to quantify the difference in the $\langle {T}_{{\rm{d}}}\rangle \mbox{--}{L}_{\mathrm{bol}}$ relation for fixed outer radii disks and radii scaling with stellar mass (see Figure 6). The fixed ${R}_{\mathrm{out}}$ models result in mass-averaged disk temperatures that decrease with stellar luminosity and become lower than 10 K for brown dwarfs. This is somewhat surprising given that dust grains in giant molecular clouds, heated by the interstellar radiation field, stabilize at ∼10 K (Mathis et al. 1983). We remind the reader that we have included interstellar radiation in our modeling (see Section 3.2), but, as also found in van der Plas et al. (2016), this extra heating does not change appreciably the dust disk temperature (see in particular their Figure 5(d)).

Zoom In Zoom Out Reset image size

In the very different assumption of outer radii scaling with stellar mass, $\langle {T}_{{\rm{d}}}\rangle$ remains confined within ∼10–20 K over four orders of magnitude in luminosity and shows an opposite behavior by slightly increasing toward the lower-luminosity/lower-mass objects. For example, the difference in $\langle {T}_{{\rm{d}}}\rangle$ for the two assumptions when ${L}_{* }=0.1{L}_{\odot }$ is 10.4 K, which results in mass estimates that differ by a factor of 9.

Disk mass estimates using a $\langle {T}_{{\rm{d}}}\rangle$ taken from a ${L}_{\mathrm{bol}}$ relationship with a fixed outer radius will ultimately result in temperatures that are too low for low-luminosity objects and temperatures that are too high for high-luminosity objects (Figure 6). Consequently, this will result in overpredicting and underpredicting the dust disk masses of low- and high-luminosity sources, respectively.

$\langle {T}_{{\rm{d}}}\rangle$ can vary by a factor of 3 to 5 times for luminosities between 0.01 and 100 ${L}_{\odot }$. If we consider VLMOs to be represented by ${L}_{\mathrm{bol}}\in [0.005,0.2]{L}_{\odot }$ (the range of our Herschel sample), we find $\langle {T}_{{\rm{d}}}\rangle$ to be lower by a factor of 8 to 2 with a mean difference (in logarithmic space) of 5 times. The disagreement in the two treatments is much less for TT objects where considering typical ${L}_{\mathrm{bol}}\in [1.5,8]{L}_{\odot }$ (comparable to the TT disks considered by Howard et al. 2013) results in an overestimate of $\langle {T}_{{\rm{d}}}\rangle$ by a factor of ∼1.5. In addition, the fixed outer radius relationship gives too low of a $\langle {T}_{{\rm{d}}}\rangle$ in the VLMO regime. Luminosities below $0.053{L}_{\odot }$ (which would include 4 of our 11 sources) result in $\langle {T}_{{\rm{d}}}\rangle$ under 5 K.

Ultimately these disagreements in $\langle {T}_{{\rm{d}}}\rangle$ are significant because they result in an even greater discrepancy in mass predictions. For example, using the fixed $\langle {T}_{{\rm{d}}}\rangle \mbox{--}{L}_{\mathrm{bol}}$ relationships and applying them to the source J04381486 ($0.005{L}_{\odot };$ ${F}_{890\mu {\rm{m}}}=6$ mJy), a total disk mass (assuming a gas-to-dust radio of 100) of 5 and $0.07{M}_{* }$ is estimated (depending on whether the vertical structure is self-consistently calculated or not), due to the low $\langle {T}_{{\rm{d}}}\rangle$ of 6.7 and 2.5 K, respectively. Table 8 compares the effects of model assumption on estimated disk masses for the five sources in our sample with $\sim 887\,\mu {\rm{m}}$ detections.

Table 8.  Effect of Model Assumptions on Total Disk Mass (Gas + Dust) Estimates

	Fixed Radius Estimates						Scaled Radius Estimates
	This Work			Andrews et al. (2013)			This Work
	Disk Mass		$\langle {T}_{{\rm{d}}}\rangle$	Disk Mass		$\langle {T}_{{\rm{d}}}\rangle$	Disk Mass		$\langle {T}_{{\rm{d}}}\rangle$
Object	(${M}_{\odot }$)	(${M}_{* }$)	(K)	(${M}_{\odot }$)	(${M}_{* }$)	(K)	(${M}_{\odot }$)	(${M}_{* }$)	(K)
CIDA-1	0.017	0.09	7.40	0.0036	0.019	16.71	0.0040	0.021	15.67
J04381486	0.274	5.05	2.45	0.0039	0.071	6.65	0.0005	0.008	20.29
GM Tau	0.008	0.11	3.25	0.0003	0.004	8.42	0.0001	0.001	18.99
CFHT-4	0.004	0.03	6.01	0.0006	0.005	14.06	0.0004	0.004	16.45
ESO Hα 559	0.069	0.48	4.96	0.0078	0.055	11.97	0.0043	0.030	17.21

Download table as:  ASCIITypeset image

Because the dust disk size likely depends on many factors, including disk mass, height, turbulent mixing, age, grain growth, and drift rates, perhaps even more so than it depends on host luminosity, characterizing the dust disk outer radius becomes critical in understanding disk masses.

For our 11 observed sources we had only two [O i] 63 μm line detections. Figure 4 shows the fit (dashed line) to TT disks found by Howard et al. (2013), which we used to set the sensitivity of our observations. Based on this, we expected more detections, but over two-thirds of our sources resulted in upper limits that fall below their fit.

In order to assess how feasible it is that these disks are underluminous (and to pinpoint a likely origin of this underluminosity), we compare our observations to the thermochemical models by Greenwood et al. (2017). These models use 2D RT and a complex chemical network on top of a parameterized disk structure, allowing investigations into the effects of disk geometry on the [O i] 63 μm line flux by computing a grid of models of varying disk masses and radii.

Figure 7 shows that for a given mass, a decreasing of the disk size leads to a decrease in [O i] emission while continuum emission increases or stays constant. This is consistent with the [O i] underluminosity of VLMOs being caused by smaller disk sizes.

Zoom In Zoom Out Reset image size

The reason for the reduced [O i] emission is due to the radii at which the emission originates. Figure 8 shows the radius within which 70% of the [O i] and continuum emissions originate normalized to the taper radius,¹⁰ a measure of the disk size. For all disk sizes, the [O i] emission originates from larger radii than the continuum. For large disks, the majority of the emission is well within the taper radius for both the [O i] 63 μm line flux and the continuum flux density. However, as the disks become smaller (going from blue to red in the figure), the area of [O i] 63 μm line emission falls beyond the taper radius, while the continuum emission remains well within. Thus, one possible contribution to the observed underluminosity is that the line-emitting area of [O i] 63 μm is disproportionately small in VLMO disks.

Zoom In Zoom Out Reset image size

Another contributing factor may be that VLMO disks tend to be less flared than their higher-mass counterparts, for which there is some observational evidence (Szűcs et al. 2010; Liu et al. 2015; however, see Harvey et al. 2012). The models by Greenwood et al. (2017) show that changing the flaring index (β) from 1.2 to 1.15, 1.10, and 1.05 decreases the [O i] 63 μm emission in a VLMO disk to 73%, 40%, and 13% of the $\beta =1.2$ levels, respectively, while the continuum emission decreases to 87%, 40%, and 19% of the $\beta =1.2$ levels (A. Greenwood 2017, personal communication). These results suggest that the [O i] 63 μm line flux is more sensitive to disk flaring than the continuum flux, causing a sample of weakly flared disks to be underluminous in their [O i] 63 μm emission.