<figure><img src="image_1.png" /></figure>

HAL
open science

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

Nathanial P. Hendler, Gijs D. Mulders, Ilaria Pascucci, Aaron Greenwood, Inga Kamp, Thomas Henning, François Ménard, William R. F. Dent, Neal J. Evans

► To cite this version:

Nathanial P. Hendler, Gijs D. Mulders, Ilaria Pascucci, Aaron Greenwood, Inga Kamp, et al. Hints for Small Disks around Very Low Mass Stars and Brown Dwarfs. The Astrophysical Journal, 2017, 841, 10.3847/1538-4357/aa71b8. insu-03692508

HAL Id: insu-03692508

https://insu.hal.science/insu-03692508v1

Submitted on 3 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Hints for Small Disks around Very Low Mass Stars and Brown Dwarfs*

Nathanial P. Hendler<fnref n="1" />, Gijs D. Mulders<fnref n="1" />,<fnref n="7" />, Ilaria Pascucci<fnref n="1" />,<fnref n="7" />, Aaron Greenwood<fnref n="2" />, Inga Kamp<fnref n="2" />, Thomas Henning<fnref n="3" />, François Ménard<fnref n="4" />, William R. F. Dent<fnref n="5" />, and Neal J. Evans II<fnref n="6" />

Received 2017 February 19; revised 2017 May 2; accepted 2017 May 2; published 2017 May 31

## Abstract

The properties of disks around brown dwarfs and very low mass stars (hereafter VLMOs) provide important boundary conditions on the process of planet formation and inform us about the numbers and masses of planets than can form in this regime. We use the Herschel Space Observatory PACS spectrometer to measure the continuum and [O I] 63 µm line emission toward 11 VLMOs with known disks in the Taurus and Chamaeleon I star-forming regions. We fit radiative transfer models to the spectral energy distributions of these sources. Additionally, we carry out a grid of radiative transfer models run in a regime that connects the luminosity of our sources with brighter T Tauri stars. We find that VLMO disks with sizes 1.3–78 au, smaller than typical T Tauri disks, fit well the spectral energy distributions assuming that disk geometry and dust properties are stellar mass independent. Reducing the disk size increases the disk temperature, and we show that VLMOs do not follow previously derived disk temperature–stellar luminosity relationships if the disk outer radius scales with stellar mass. Only 2 out of 11 sources are detected in [O I] despite a better sensitivity than was achieved for T Tauri stars, suggesting that VLMO disks are underluminous. Using thermochemical models, we show that smaller disks can lead to the unexpected [O I] 63 µm nondetections in our sample. The disk outer radius is an important factor in determining the gas and dust observables. Hence, spatially resolved observations with ALMA—to establish if and how disk radii scale with stellar mass—should be pursued further.

**Key words:** brown dwarfs – circumstellar matter – protoplanetary disks

**Supporting material:** machine-readable table

## 1. Introduction

Brown dwarfs are substellar objects that represent the low mass end of the star formation process. Similar to their stellar counterparts, T Tauri (TT) stars, young (1–2 Myr) brown dwarfs possess protoplanetary disks of dust and gas (e.g., Luhman et al. 2007) that are the environment of planet formation (e.g., Apai et al. 2005). Indeed, the recent discovery of the TRAPPIST-1 exoplanet system at the substellar boundary (Gillon et al. 2016, 2017) suggests that brown dwarfs may be capable of forming planetary systems. Observing protoplanetary disks in this regime allows us to observe the lower-boundary conditions of planet formation.

The conditions in protoplanetary disks are strongly dependent on stellar mass. In particular, scaling laws have been identified between host star mass and protoplanetary dust disk mass (e.g., Pascucci et al. 2016). Planet formation models typically adopt these scaling laws to predict the planets that may form around low mass stars (e.g., Raymond et al. 2007; Mordasini et al. 2012). However, disks around the lowest-mass stars may not be scaled-down versions of their higher-mass counterparts. Disks around brown dwarfs show differences in their dust and gas evolution (Pascucci et al. 2009). There are

also indications that disks around brown dwarfs have a different disk geometry (e.g., Szűcs et al. 2010; Daemgen et al. 2016), though this is not always found (Harvey et al. 2010; Mulders & Dominik 2012).

Because brown dwarfs are significantly fainter than TT stars, only limited observations are present to directly constrain their dust masses (e.g., Mohanty et al. 2013; Pascucci et al. 2016, and references therein). Several studies have used Herschel far-infrared photometry (e.g., Harvey et al. 2012; Daemgen et al. 2016; van der Plas et al. 2016) and ALMA (e.g., van der Plas et al. 2016) to estimate dust masses. One parameter that has remained largely unexplored in these studies is the dust disk outer radius. From a theoretical perspective, disks around brown dwarfs are predicted to form with smaller gas disks (Bate 2012) and experience more efficient radial drift, leading to smaller dust disks (Pinilla et al. 2013). Observations with ALMA of a handful of spatially resolved brown dwarf dust and gas disks range in size from 20 to 140 au (Ricci et al. 2013, 2014; Testi et al. 2016), smaller than TT disks (22–440 au with a mean size of 165 au; Isella et al. 2009; Andrews et al. 2010; Guilloteau et al. 2011).

In this work, we explore the impact of smaller disk radii around brown dwarfs and very low mass stars (hereafter VLMOs) on observables related to the dust and the gas in those disks. We present a Herschel survey of [O I] 63 µm for 11 VLMOs probing the disk gas in Sections 2 and 3.1. Detailed thermochemical models of protoplanetary disks have shown

<footnote n="1">¹Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA; equant@lpl.arizona.edu</footnote>

<footnote n="2">²Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="3">³Max Planck Institute for Astronomy, Konigstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="4">⁴Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="5">⁵Department of Engineering, Atacama Large Millimeter/submillimeter Array (ALMA) Santiago Central Offices, Alonso de Córdova 3107, Vitacura, Casilla 763 0355, Santiago, Chile</footnote>

<footnote n="6">⁶Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA</footnote>

<footnote>* Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>

<footnote n="7">⁷Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science.</footnote>
<table><caption>**Table 1**
Source Properties</caption><thead><tr><th>2MASS</th><th>Common Name</th><th>Region</th><th>SpTy</th><th>References</th><th>A<sub>v</sub></th><th>log L<sub>bol</sub> (L<sub>☉</sub>)</th><th>References</th><th>T<sub>eff</sub> (K)</th><th>M<sub>☉</sub> (M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>J04141760+2806096</td><td>CIDA-1</td><td>Taurus</td><td>M5</td><td>L17</td><td>4.6</td><td>-0.7</td><td>R10</td><td>3125</td><td>0.19</td></tr><tr><td>J04193545+2827218</td><td>FR Tau</td><td>Taurus</td><td>M5.25</td><td>L17</td><td>0.0</td><td>-0.9</td><td>R10</td><td>3090</td><td>0.14</td></tr><tr><td>J04233539+2503026</td><td>FU Tau A</td><td>Taurus</td><td>M7</td><td>L17</td><td>2.0</td><td>-0.72</td><td>L09</td><td>2880</td><td>0.16</td></tr><tr><td>J04295950+2433078</td><td>CFHT-20</td><td>Taurus</td><td>M5</td><td>L17</td><td>4.6</td><td>-0.7</td><td>R10</td><td>3125</td><td>0.19</td></tr><tr><td>J04381486+2611399</td><td>J04381486</td><td>Taurus</td><td>M7.25</td><td>L17</td><td>3.5</td><td>-2.3</td><td>R10</td><td>2846</td><td>0.05</td></tr><tr><td>J04382134+2609137</td><td>GM Tau</td><td>Taurus</td><td>M5</td><td>L17</td><td>2.0</td><td>-1.15<sup>a</sup></td><td>H08</td><td>3125</td><td>0.14</td></tr><tr><td>J04393364+2359212</td><td>J04393364</td><td>Taurus</td><td>M5</td><td>L17</td><td>1.3</td><td>-1.0</td><td>R10</td><td>3125</td><td>0.18</td></tr><tr><td>J04394748+2601407</td><td>CFHT-4</td><td>Taurus</td><td>M7</td><td>L17</td><td>4.4</td><td>-1.0</td><td>R10</td><td>2880</td><td>0.11</td></tr><tr><td>J11062554-7633418</td><td>ESO Hα 559</td><td>Cha I</td><td>M5.25</td><td>L07</td><td>3.58<sup>b</sup></td><td>-1.28</td><td>L07</td><td>3088</td><td>0.14</td></tr><tr><td>J11071668-7735532</td><td>Cha Hα 1</td><td>Cha I</td><td>M7.75</td><td>L07</td><td>0.0<sup>b</sup></td><td>-1.82</td><td>L07</td><td>2765</td><td>0.04</td></tr><tr><td>J11105597-7645325</td><td>Hn 13</td><td>Cha I</td><td>M5.75</td><td>L07</td><td>2.91<sup>b</sup></td><td>-0.89</td><td>L07</td><td>3021</td><td>0.14</td></tr></tbody></table>

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

a GM Tau luminosity from D14.

b L07 report the extinction at J band. From those values we computed the A<sub>v</sub> in this table using the interstellar extinction law in Cardelli et al. (1989) and an R<sub>v</sub> 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.

that this transition can be used as an order-of-magnitude disk mass estimator (e.g., Woitke et al. 2010); hence, a large survey toward TT stars (Dent et al. 2013) was carried out. We perform radiative transfer (RT) modeling of the spectral energy distribution (SED) of these sources in Section 3.2, starting from the zero-order assumption that the disk geometry and dust opacity are stellar mass independent. We also run two RT grids with and without a stellar-mass-dependent disk radius to quantify the difference in disk temperature in Section 3.3. We examine the underluminosity of [O I] 63 µm in our VLMO disks in Section 4.1, while in Section 4.2 we compare the likelihoods of small and large disk models by testing model grids against the SEDs of our VLMO sample. We also show how smaller disk radii lead to higher disk temperatures than previously assumed and impact estimates of the dust disk mass.

## 2. Observations and Data Reduction

### 2.1. Sample

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 possess a dust disk (based on infrared excess emission) and to have a flux density at 24 µm (F<sub>24</sub>) 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<sub>70</sub>), 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 Hertz-sprung-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<sub>☉</sub> as demonstrated by Stassun et al. (2014) using a sample of eclipsing binary systems.

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

### 2.2. Herschel/PACS Spectroscopy

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"/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
<table><caption>**Table 2**
Herschel/PACS Sample and Observations</caption><thead><tr><th>Target Name</th><th>ObsID</th><th>Obs. Date</th><th>Duration (s)</th></tr></thead><tbody><tr><td>CIDA-1</td><td>1342268646</td><td>2013 Mar 25</td><td>9868.0</td></tr><tr><td>FR Tau</td><td>1342263516</td><td>2013 Feb 12</td><td>16420.0</td></tr><tr><td>FU Tau A</td><td>1342264241</td><td>2013 Feb 25</td><td>16420.0</td></tr><tr><td>CFHT-20</td><td>1342265469</td><td>2013 Feb 16</td><td>16420.0</td></tr><tr><td>J04381486</td><td>1342265470</td><td>2013 Feb 16</td><td>16420.0</td></tr><tr><td>GM Tau</td><td>1342264239</td><td>2013 Feb 25</td><td>16420.0</td></tr><tr><td>J04393364</td><td>1342263934</td><td>2013 Feb 19</td><td>16420.0</td></tr><tr><td>CFHT-4</td><td>1342264240</td><td>2013 Feb 25</td><td>16420.0</td></tr><tr><td>ESO Hα 559</td><td>1342263489</td><td>2013 Feb 11</td><td>16420.0</td></tr><tr><td>Cha Hα 1</td><td>1342263459</td><td>2013 Feb 10</td><td>16420.0</td></tr><tr><td>Hn 13</td><td>1342263492</td><td>2013 Feb 11</td><td>16420.0</td></tr></tbody></table>

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 \times 10^{-18}$ W m$^{-2}$. The exposure time was set to 16,420 s to achieve a $3\sigma$ detection on this line, i.e., a sensitivity of $1.4 \times 10^{-18}$ W m$^{-2}$, 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.

### 3. Analysis

Our analysis is separated into three subsections. In Section 3.1 we present our findings for the 63 µm [O I] line

emission and continuum and 189 µm continuum for 11 Herschel Space Observatory observations. The [O I] 63 µm emission informs us about the gas content, especially in the disk atmosphere (e.g., Woitke et al. 2010), while the continuum emission informs us about the dust in the disk.

This is followed by two separate sets of RT modeling schemes. The first, Section 3.2, investigates VLMO dust disk properties by fitting our source SEDs. Section 3.3 describes the grid of models we used to analyze the relationship between the disk outer radius and the derived disk mass.

#### 3.1. [O I] Line and Continuum Detections

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\sigma$ 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\sigma$ uncertainty. In case of nondetections, we fit the same spectral range with a first-order polynomial, and we provide in Table 3 the $3\sigma$ 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 ~5σ [O I] 63 µm detection ($5.17 \pm 1.04 \times 10^{-18}$ W m$^{-2}$), while Riviere-Marichalar et al. (2016) report a $3\sigma$ upper limit of ($<3 \times 10^{-18}$ W 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\sigma$ of the uncertainties we quote. In the case of CIDA-1 the 63 and 70 µm flux densities are within $2\sigma$. 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.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Herschel/PACS spectroscopy measurements of the wavelength region around the [O I] 63 µm emission line (solid gray line). Greater than 3σ Gaussian fits (detections) are overplotted as solid green lines. Hypothetical 3σ upper limits (FWHM of unresolved line) are shown as red dashed lines.</figcaption></figure>

## 3.2. Continuum Radiative Transfer Modeling of Individual Sources

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.

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_{gas} = T_{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
<table><caption>**Table 3**
Observed Continuum Flux Densities and Line Fluxes</caption><thead><tr><th>Source</th><th>F<sub>[OII]</sub><br>(10<sup>-18</sup> W m<sup>-2</sup>)</th><th>F<sub>63</sub><br>(mJy)</th><th>F<sub>189</sub><br>(mJy)</th><th>F<sub>70</sub><br>(mJy)</th><th>References<br>(F<sub>70</sub>)</th></tr></thead><tbody><tr><td>CIDA-1</td><td>&le;3.01</td><td>358.7 &plusmn; 50.9</td><td>270.4 &plusmn; 83.3</td><td>266 &plusmn; 2</td><td>B14</td></tr><tr><td>FR Tau</td><td>&le;2.62</td><td>&le;120.7</td><td>&le;145.4</td><td>46 &plusmn; 3</td><td>B14</td></tr><tr><td>FU Tau A</td><td>5.17 &plusmn; 1.04</td><td>&le;98.9</td><td>&le;131.3</td><td>86 &plusmn; 17</td><td>R10</td></tr><tr><td>CFHT-20</td><td>&le;2.89</td><td>&le;127.6</td><td>&le;154.5</td><td>128 &plusmn; 4</td><td>B14</td></tr><tr><td>J04381486</td><td>&le;2.91</td><td>&le;132.0</td><td>&le;152.5</td><td>95 &plusmn; 2</td><td>B14</td></tr><tr><td>GM Tau</td><td>&le;2.54</td><td>&le;121.5</td><td>&le;66.8</td><td>36 &plusmn; 2</td><td>B14</td></tr><tr><td>J04393364</td><td>&le;2.28</td><td>100.3 &plusmn; 32.4</td><td>&le;140.0</td><td>70 &plusmn; 1</td><td>B14</td></tr><tr><td>CFHT-4</td><td>&le;2.73</td><td>&le;110.6</td><td>&le;161.9</td><td>109 &plusmn; 5</td><td>B14</td></tr><tr><td>ESO H&alpha; 559</td><td>&le;3.30</td><td>174.0 &plusmn; 45.6</td><td>317.7 &plusmn; 42.2</td><td>...</td><td>...</td></tr><tr><td>Cha H&alpha; 1</td><td>&le;2.55</td><td>&le;112.5</td><td>172.6 &plusmn; 40.0</td><td>...</td><td>...</td></tr><tr><td>Hn 13</td><td>4.35 &plusmn; 1.30</td><td>&le;114.0</td><td>&le;76.7</td><td>...</td><td>...</td></tr></tbody></table>

Note. 3σ upper limits are reported for nondetections.

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° 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_{\text{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 $\Sigma \propto r^{-1}$, while the minimum and maximum grain sizes are $a = 0.1 \ \mu\text{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^{-2}$ and $10^{-9} M_{\odot}$ are sampled. The disk inner radius is explored using 11 logarithmically spaced steps between $10^{-4}$ and $10^1$ 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_{\text{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 $\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_{\text{out}}$, $R_{\text{in}}$, and $M_{\text{dust}}$ are estimated from Bayesian probability distributions and reported in Table 5. Figure 3 shows the probability distribution and confidence interval for $R_{\text{out}}$, for each source. See Appendix B for additional information on the confidence intervals and for examples of probability distributions for $R_{\text{in}}$ and $M_{\text{dust}}$.

Along with the resulting parameters of our best-fit models, we also report the mass-averaged disk temperature ($\langle T_d \rangle$) in Table 5. Here, we take $\langle T_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_d \rangle$.

Because we find maximum likelihood values of $R_{\text{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_{\text{out}}$ fixed to 200 au, a typical value used for modeling TT disks. These models use the same 20 and 11 steps of $M_{\text{dust}}$ and $R_{\text{in}}$, respectively, outlined above. We use these grids to similarly derive a $\langle T_d \rangle$ with which to compare our first grid of variable $R_{\text{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
<figure><img src="image_5.png" /><figcaption>**Figure 2.** SEDs of our Herschel sources. Orange: photometric fluxes from literature; triangles: upper limits. The best-fit RT models for each source are displayed for each target (black line).</figcaption></figure>

compute the temperature in the outer disk (see Woitke et al. 2009).

### 3.3. Outer Radius Parameter Study

In addition to modeling individual sources, we carry out two sets of RT models to assess the impact of the disk outer radius ($R_{\text{out}}$) on the derived disk mass and $\langle T_{\text{d}} \rangle$. Both parameter studies are calculated for a range of stellar masses, but in one case $R_{\text{out}}$ is held fixed at 200 au (hereafter, fixed $R_{\text{out}}$ study), while in the other case (hereafter, variable $R_{\text{out}}$ study) $R_{\text{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).

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_{\odot}$, which cover from 0.019 to 6.33 $L_{\odot}$.

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_{\text{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.

## 4. Results

### 4.1. Faint [O I] 63 µm Emission from Very Low Mass Star and Brown Dwarf Disks

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.

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<fnref n="8" /> (see Figure 4).

<footnote n="8">⁸ We consider only nondetections from our sample. Of the two VLMOs with detected [O I] lines, one is excluded from our analysis because it has a known jet, and the other is excluded because it has only an upper limit for its 63 µm continuum.</footnote>
<table><caption>**Table 4**
Source Photometry</caption><thead><tr><th>Source</th><th>Wavelength (μm)</th><th>Flux (mJy)</th><th>References</th></tr></thead><tbody><tr><td>CIDA-1</td><td>0.88, 1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 70.00, 100.00, 160.00, 189.00, 890.00, 1300.00, 2600.00</td><td>0.0047, 0.033, 0.06, 0.075, 0.075, 0.12, 0.12, 0.16, 0.2, 0.28, 0.27, 0.29, 0.25, 2.7e+02, 0.035, 0.013, 0.0078</td><td>L00, C03, W10, L10, R10, B14, D16, H17, A13, S09</td></tr><tr><td>FR Tau</td><td>1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 71.42, 160.00, 189.00, 890.00, 1300.00</td><td>0.066, 0.073, 0.068, 0.048, 0.052, 0.056, 0.081, 0.09, 0.083, 0.043, 0.033, 1.5e+02, 0.039, 0.0015</td><td>C03, W10, L10, R10, B14, H17, A13</td></tr><tr><td>FU Tau A</td><td>1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 71.42, 100.00, 160.00, 189.00, 1300.00</td><td>0.078, 0.11, 0.12, 0.11, 0.14, 0.13, 0.15, 0.12, 0.1, 0.086, 0.039, 0.25, 1.3e+02, 0.002</td><td>C03, W10, R10, B14, H17, A13</td></tr><tr><td>CFHT-20</td><td>0.88, 1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 71.42, 160.00, 189.00, 890.00, 1300.00</td><td>0.0023, 0.034, 0.063, 0.079, 0.073, 0.068, 0.055, 0.049, 0.057, 0.078, 0.12, 0.091, 1.5e+02, 0.0089, 0.0034</td><td>G06, C03, W10, L10, R10, B14, H17, A13</td></tr><tr><td>J04381486</td><td>0.68, 0.88, 1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 70.00, 160.00, 189.00, 890.00, 1300.00</td><td>2e-05, 0.00019, 0.0014, 0.0023, 0.0043, 0.011, 0.016, 0.017, 0.018, 0.022, 0.073, 0.095, 0.067, 1.5e+02, 0.006, 0.0023</td><td>G07, C03, W10, L10, R10, B14, H17, A13</td></tr><tr><td>GM Tau</td><td>1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 100.00, 160.00, 189.00, 860.00, 1300.00, 2600.00</td><td>0.012, 0.024, 0.037, 0.049, 0.052, 0.045, 0.045, 0.046, 0.053, 0.036, 0.031, 67, 0.00087, 0.0048, 0.0024</td><td>C03, W10, L10, R10, D16, H17, M13, A13, S09</td></tr><tr><td>J04393364</td><td>1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 70.00, 160.00, 189.00, 890.00, 1300.00</td><td>0.038, 0.049, 0.052, 0.038, 0.039, 0.034, 0.044, 0.051, 0.059, 0.07, 0.044, 1.4e+02, 0.009, 0.0034</td><td>C03, W10, L10, R10, B14, H17, A13</td></tr><tr><td>CFHT-4</td><td>0.68, 0.88, 1.25, 1.75, 2.50, 3.35, 4.49, 5.73, 7.87, 11.56, 23.67, 70.00, 160.00, 189.00, 850.00, 890.00, 1300.00</td><td>0.0001, 0.0013, 0.022, 0.04, 0.049, 0.037, 0.043, 0.042, 0.048, 0.048, 0.075, 0.11, 0.15, 1.6e+02, 0.011, 0.0043, 0.002</td><td>M01, L00, C03, W10, L10, R10, B14, H17, K03, R14, A13</td></tr><tr><td>ESO Hα 559</td><td>1.25, 1.75, 2.50, 3.40, 4.50, 5.80, 8.00, 12.00, 24.00, 100.00, 160.00, 189.00, 850.00, 890.00</td><td>0.025, 0.029, 0.024, 0.015, 0.014, 0.013, 0.013, 0.016, 0.055, 0.23, 0.28, 3.2e+02, 47, 0.044</td><td>C03, C13, L08, B14, H17, P16, B11</td></tr><tr><td>Cha Hα 1</td><td>1.25, 1.75, 2.50, 3.40, 4.50, 5.80, 8.00, 12.00, 24.00, 189.00</td><td>0.0073, 0.0088, 0.009, 0.0061, 0.0061, 0.006, 0.0081, 0.013, 0.028, 1.7e+02</td><td>C03, C13, L08, H17</td></tr><tr><td>Hn 13</td><td>1.25, 1.75, 2.50, 3.40, 4.50, 5.80, 8.00, 12.00, 24.00, 189.00, 850.00</td><td>0.067, 0.079, 0.079, 0.052, 0.049, 0.047, 0.044, 0.039, 0.07, 77, 2.2</td><td>C03, C13, L08, H17, P16</td></tr></tbody></table>

Note. Photometric observations used to create SEDs for model fitting.

References. B14: Bulger et al. 2014; G06: Guieu et al. 2006; L10: Luhman et al. 2010; M01: Martín et al. 2001; R10: Rebull et al. 2010; C03: Cutri et al. 2003; W10: Wright et al. 2010; A13: Andrews et al. 2013; M13: Mohanty et al. 2013; S09: Schaefer et al. 2009; K03: Klein et al. 2003; G07: Guieu et al. 2007; L00: Luhman 2000; P16: Pascucci et al. 2016; H17: this work.

(This table is available in its entirety in machine-readable form.)

To do this, we use the Bayesian method of regression fitting described in Kelly (2007),<fnref n="9" /> 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.

## 4.2. Small Dust Disks

Our RT modeling suggests greater likelihoods for outer disk radii that are smaller (maximum likelihoods ≤78 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\pm14_{12}$ au for CIDA-1 and greater than 80 au for CFHT-4 (Ricci et al. 2014). ESO Hα 559 was resolved at 887 µm, and a fitted elliptical Gaussian resulted in an FWHM of $0.23 \times 0.15$ (37 au $\times$ 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

<footnote n="9">⁹ Implemented by Josh Meyers in Python (2016 January 16 commit).</footnote>
<table><caption>**Table 5**
Maximum Likelihood Parameters for Source Models</caption><thead><tr><th rowspan="2">Source</th><th rowspan="2"><i>i</i><br>(deg)</th><th colspan="2"><i>R</i><sub>in</sub></th><th colspan="2"><i>R</i><sub>out</sub></th><th colspan="2">Log <i>M</i><sub>dust</sub></th><th rowspan="2">⟨<i>T</i><sub>d</sub>⟩<br>(K)</th></tr><tr><th>(au)</th><th></th><th>(au)</th><th></th><th>(M<sub>☉</sub>)</th><th></th></tr></thead><tbody><tr><td>CIDA-1</td><td>35</td><td>0.32</td><td>+1.5<br>-0.31</td><td>7.6</td><td>+60<br>-3.9</td><td>-4.9</td><td>+0.55<br>-0.92</td><td>48</td></tr><tr><td>FR Tau</td><td>77</td><td>0.32</td><td>+1.5<br>-0.32</td><td>3.2</td><td>+25<br>-2.3</td><td>-6.1</td><td>+0.92<br>-1.3</td><td>46</td></tr><tr><td>FU Tau A</td><td>56</td><td>0.00032</td><td>+0.017<br>-0.00026</td><td>1.3</td><td>+36<br>-0.47</td><td>-7.2</td><td>+1.3<br>-0.92</td><td>140</td></tr><tr><td>CFHT-20</td><td>56</td><td>0.0001</td><td>+0.018<br>-4.4×10<sup>-5</sup></td><td>4.3</td><td>+46<br>-2.7</td><td>-5.7</td><td>+0.92<br>-1.3</td><td>64</td></tr><tr><td>J04381486</td><td>70</td><td>0.032</td><td>+0.025<br>-0.031</td><td>10</td><td>+5.6<br>-3.6</td><td>-4.9</td><td>+0.55<br>-0.18</td><td>33</td></tr><tr><td>GM Tau</td><td>35</td><td>0.0032</td><td>+0.015<br>-0.0031</td><td>1.8</td><td>+3.2<br>-0.92</td><td>-6.8</td><td>+3.1<br>-0.55</td><td>73</td></tr><tr><td>J04393364</td><td>35</td><td>0.01</td><td>+0.55<br>-0.0094</td><td>4.3</td><td>+4.6<br>-3.4</td><td>-5.7</td><td>+2.4<br>-1.3</td><td>42</td></tr><tr><td>CFHT-4</td><td>80</td><td>0.01</td><td>+0.0078<br>-0.0099</td><td>78</td><td>+43<br>-66</td><td>-5.3</td><td>+0.18<br>-0.92</td><td>15</td></tr><tr><td>ESO Hα 559</td><td>35</td><td>0.032</td><td>+0.53<br>-0.03</td><td>33</td><td>+18<br>-17</td><td>-4.2</td><td>+1.7<br>-0.18</td><td>18</td></tr><tr><td>Cha Hα 1</td><td>35</td><td>0.1</td><td>+0.46<br>-0.099</td><td>33</td><td>+88<br>-24</td><td>-3.1</td><td>+1.3<br>-0.55</td><td>16</td></tr><tr><td>Hn 13</td><td>35</td><td>0.01</td><td>+0.0078<br>-0.0099</td><td>3.2</td><td>+13<br>-2.3</td><td>-6.1</td><td>+0.92<br>-1.3</td><td>64</td></tr></tbody></table>

Note. The maximum likelihood values are shown for the parameters $R_{in}$, $R_{out}$, and $M_{dust}$, along with their corresponding $1\sigma$ confidence intervals. $\langle T_d \rangle$ is calculated from the best-fit model selected for the source as described in Section 3.2.

well” with *Hubble* Wide Field Planetary Camera observations. These values are larger than our finding of $R_{out} = 10^{+5.6}_{-3.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_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_d \rangle - L_{bol}$ relations by Andrews et al. (2013) and by our fixed $R_{out}$ models (Section 3.3). In both models the dust disk radius is 200 au regardless of stellar luminosity/mass, but our fixed $R_{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_d \rangle$ for VLMO disks than typically assumed. The SED fits with fixed $R_{out}$ have temperatures in between those predicted by the Andrews et al. (2013) and our fixed $R_{out}$ grid relations (right panel of Figure 5).

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_{out}$ grid VLMO models are found to be optically thick at 63 µm for all of our sources with the exception of CFHT-4 (the source with the largest $R_{out}$) and at 850 µ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 µ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.

### 4.3. Disk Temperature–Stellar Luminosity Relation

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_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_d \rangle - L_{bol}$ relation for fixed outer radii disks and radii scaling with stellar mass (see Figure 6). The fixed $R_{out}$ models result in mass-averaged disk tempera- tures 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)).

In the very different assumption of outer radii scaling with stellar mass, $\langle T_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_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.

## 5. Discussion

In this paper we start from the zero-order assumptions that disk geometry and dust opacity are stellar mass independent. With these assumptions, our RT modeling of each source
<figure><img src="image_6.png" /><figcaption>Figure 3. Bayesian probability distribution of R<sub>out</sub>, values. Vertical blue lines are 1σ confidence intervals reported in Table 5.</figcaption></figure>

<table><caption>**Table 6**
Parameter Study Model Inputs</caption><thead><tr><th>Parameter</th><th>Variable R<sub>out</sub></th><th>Fixed R<sub>out</sub></th></tr></thead><tbody><tr><td>R<sub>out</sub>(au)</td><td>200(M<sub>*</sub>/M<sub>☉</sub>)</td><td>200</td></tr><tr><td>M<sub>dust</sub></td><td>~M<sub>*</sub><sup>2.0</sup></td><td>~M<sub>*</sub><sup>1.2</sup></td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>0.05–2.239</td><td></td></tr><tr><td>R<sub>in</sub>(au)</td><td>0.05 × L<sub>*</sub><sup>0.5</sup></td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 4. Fluxes are scaled to a distance of 140 pc as in Howard et al. (2013). Red points are [O I] detections from this work. Black points are this work, with arrows indicating upper limits in continuum and/or line flux observations. Gray points are full disks from Howard et al. (2013) (squares are detections, and triangles are upper limits), with the dashed line showing their best fit. The green line is a refitting of the Howard et al. (2013) data using their upper limits. The blue line is the fit including our sources. Shaded areas are 1σ confidence intervals.</figcaption></figure>

<table><caption>**Table 7**
[O I] vs. 70 µm Continuum Fitting Results</caption><thead><tr><th rowspan="2">Description</th><th colspan="2">Linear Regression</th></tr><tr><th>Intercept</th><th>Slope</th></tr></thead><tbody><tr><td>TT Disks</td><td>0.89 ± 0.18</td><td>-16.74 ± 0.06</td></tr><tr><td>TT + VLMO</td><td>1.14 ± 0.22</td><td>-16.78 ± 0.09</td></tr></tbody></table>

suggests that dust disk outer radii range from 1.3 to 78 au. This is significantly less than the 200 au radius typically assumed for TT disks. Although small disk radii may not be a unique solution, they agree well with the few previously resolved VLMO disks. As such, it becomes necessary to consider the effect of smaller disk sizes when calculating their masses and interpreting the [O I] 63 µm line emission.

### 5.1. Mass Estimates for Small Disks

Disk mass estimates using a $\langle T_d \rangle$ taken from a $L_{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_d \rangle$ can vary by a factor of 3 to 5 times for luminosities between 0.01 and 100 $L_{⊙}$. If we consider VLMOs to be represented by $L_{bol} \in [0.005, 0.2]L_{⊙}$ (the range of our Herschel sample), we find $\langle T_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_{bol} \in [1.5, 8]L_{⊙}$ (comparable to the TT disks considered by Howard et al. 2013) results in an overestimate of $\langle T_d \rangle$ by a factor of ~1.5. In addition, the fixed outer radius relationship gives too low of a $\langle T_d \rangle$ in the VLMO
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Circles are outputs from the RT modeling of our *Herschel* sources. The left panel shows the best fit of these models where the outer disk radius is a free parameter. The right panel shows the best fits where disk radius is fixed at 200 au. Both panels show the best-fit temperature–luminosity relation of our fixed radius parameter study models (lower solid line) and the Andrews et al. (2013) fit (dotted line).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** $R_{\text{out}} \propto M_{\text{star}}$ shows the best fit of our grid of RT models where the disk’s outer radius is scaled to the stellar mass. $R_{\text{out}} = 200$ au shows the best fit of our grid of RT models where the disk’s outer radius is fixed to 200 au (it is the same line seen in the right and left panels of Figure 5). Red filled circles are our best-fit models (Figure 5, left panel), and the blue filled circles are the fixed $R_{\text{out}}$ grid (200 au) models (Figure 5, right panel).</figcaption></figure>

regime. Luminosities below $0.053 L_{\odot}$ (which would include 4 of our 11 sources) result in $\langle T_d \rangle$ under 5 K.

Ultimately these disagreements in $\langle T_d \rangle$ are significant because they result in an even greater discrepancy in mass predictions. For example, using the fixed $\langle T_d \rangle - L_{\text{bol}}$ relationships and applying them to the source J04381486 ($0.005 L_{\odot}$; $F_{890 \mu 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_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 ~887 $\mu$m detections.

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.

## 5.2. [O I] Emission in Small Disks

For our 11 observed sources we had only two [O I] $63 \mu 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 thermo-chemical 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 \mu 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.

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,<fnref n="10" /> 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 \mu 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 \mu m$ line emission falls beyond the taper radius, while the continuum emission remains well within. Thus, one possible contribution to

<footnote n="10">¹⁰ Taper radius is the location of an exponential tapering off of the disk density power law, $R_{\text{taper}} = R_{\text{out}}/8$.</footnote>
<table><caption>**Table 8**
Effect of Model Assumptions on Total Disk Mass (Gas + Dust) Estimates</caption><thead><tr><th rowspan="3">Object</th><th colspan="6">Fixed Radius Estimates</th><th colspan="3">Scaled Radius Estimates</th></tr><tr><th colspan="3">This Work</th><th colspan="3">Andrews et al. (2013)</th><th colspan="3">This Work</th></tr><tr><th colspan="2">Disk Mass</th><th rowspan="2">⟨T<sub>d</sub>⟩<br>(K)</th><th colspan="2">Disk Mass</th><th rowspan="2">⟨T<sub>d</sub>⟩<br>(K)</th><th colspan="2">Disk Mass</th><th rowspan="2">⟨T<sub>d</sub>⟩<br>(K)</th></tr><tr><th>(M<sub>☉</sub>)</th><th>(M<sub>*</sub>)</th><th>(M<sub>☉</sub>)</th><th>(M<sub>*</sub>)</th><th>(M<sub>☉</sub>)</th><th>(M<sub>*</sub>)</th></tr></thead><tbody><tr><td>CIDA-1</td><td>0.017</td><td>0.09</td><td>7.40</td><td>0.0036</td><td>0.019</td><td>16.71</td><td>0.0040</td><td>0.021</td><td>15.67</td></tr><tr><td>J04381486</td><td>0.274</td><td>5.05</td><td>2.45</td><td>0.0039</td><td>0.071</td><td>6.65</td><td>0.0005</td><td>0.008</td><td>20.29</td></tr><tr><td>GM Tau</td><td>0.008</td><td>0.11</td><td>3.25</td><td>0.0003</td><td>0.004</td><td>8.42</td><td>0.0001</td><td>0.001</td><td>18.99</td></tr><tr><td>CFHT-4</td><td>0.004</td><td>0.03</td><td>6.01</td><td>0.0006</td><td>0.005</td><td>14.06</td><td>0.0004</td><td>0.004</td><td>16.45</td></tr><tr><td>ESO Hα 559</td><td>0.069</td><td>0.48</td><td>4.96</td><td>0.0078</td><td>0.055</td><td>11.97</td><td>0.0043</td><td>0.030</td><td>17.21</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Line and continuum emission from thermochemical models by Greenwood et al. (2017) (color filled circles) superimposed on data and fits described in Figure 4. Isarithms of constant mass (1, 4, and 8 × 10⁻⁴ M☉) and flaring index show how smaller disks lead to a rapid decrease in [O I] line flux.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Line and continuum emission from thermochemical models by Greenwood et al. (2017) showing that 70% of [O I] 63 µm emission comes from larger relative distances as the disk size decreases. 70% of the given emission ([O I] or continuum) originates from within the R₇₀ radius.</figcaption></figure>

the observed underluminosity is that the line-emitting area of [O I] 63 µm is disproportionately small in VLMO disks.

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 β = 1.2 levels, respectively, while the continuum emission decreases to 87%, 40%, and 19% of the β = 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.

## 6. Summary and Conclusions

Starting from a sample of 11 luminosity-selected VLMOs and brown dwarfs with disks, we use the PACS spectrometer aboard the Herschel Space Observatory to measure the [O I] 63 µm line fluxes and model the spectral energy distributions. The analysis of 10 out of our 11 sources is consistent with the previous analysis by Riviere-Marichalar et al. (2016); however, using a newer data reduction pipeline, we report a new [O I] 63 µm detection toward FU Tau A.

Our main findings are as follows:

1. We detect only two sources in [O I] 63 µm despite a 3 times better sensitivity than for TT disks, suggesting that disks around VLMO and brown dwarf disks are underluminous in [O I] 63 µm.

2. Assuming that the disk geometry and opacity are stellar mass independent, we find that disk models with outer radii in the range of 1.3–78 au have the highest maximum likelihood values when compared to SEDs. These radii are significantly smaller than those of TT star disks (22–440 au).

3. If VLMO and brown dwarf disks are indeed smaller and hotter than TT disks, they do not follow previously derived temperature–luminosity relationships. A disk outer radius that scales linearly with stellar mass results in almost a flat ⟨T<sub>d</sub>⟩−L<sub>bol</sub> relationship. This results in higher average temperatures for VLMO disks and smaller disk mass estimates.

4. Using thermochemical models, we find that smaller disks result in lower [O I] 63 µm fluxes, which may explain the nondetections in the sample.

A smaller size of protoplanetary disks around substellar objects has a large impact on the derived disk masses from large surveys of star-forming regions. High-resolution ALMA observations are necessary to quantify how the disk radius scales with stellar mass to correctly gauge the planet-forming potential of VLMOs and brown dwarfs.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Bayesian probability distributions of $M_{\text{dust}}$, $R_{\text{in}}$, and $R_{\text{out}}$ for J04381486 and Cha Hα 1. For these examples, we show all three of our free parameters. The vertical dashed blue lines denote our 68.27% confidence interval.</figcaption></figure>

The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of 100 subcontractors. HIPE is a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS, and SPIRE consortia. An allocation of computer time from the UA Research Computing High Performance Computing (HPC) and High Throughput Computing (HTC) at the University of Arizona is gratefully acknowledged. This material is based on work supported by the National Science Foundation under grant no. 1228509. Support for this work, part of the NASA Herschel Science Center Theoretical Research/Laboratory Astrophysics Program, was provided by NASA through a contract issued by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA (grant no. 1483963).

This research has made use of the SIMBAD database and the VizieR catalog access tool, CDS, Strasbourg, France.

Facility: Herschel (PACS).

# Appendix A
## Confidence Interval Calculations

Because some of our probability distributions are asymmetric, or even multimodal, we use a modified version of the shortest interval method for defining our confidence intervals

(e.g., Barlow 1989; Andrae 2010). We restrict our confidence intervals to being continuous and move the intervals away from the probability peak such that higher probabilities are preferred until at least 68.27% of the pdf area is found.

Examples of our confidence intervals are shown for the two sources J04381486 and Cha Hα 1 in Figure 9. These two sources represent two extrema in the SED coverage and illustrate how our confidence intervals prefer areas of higher probability (e.g., the $R_{\text{in}}$ panels illustrate how the confidence intervals do not include the low probabilities to the right of the peak probabilities) and the importance of millimeter photometry to constraining the disk outer radius. Cha Hα 1 lacks photometry data (detections or upper limits) beyond 100 µm, and consequently the outer radius (and dust mass) is unconstrained by our modeling.

# Appendix B
## Effect of Outer Radius on Spectral Energy Distribution

It has been suggested that the outer radius disk parameter has a minimal effect on the SEDs of brown dwarf disks (e.g., Harvey et al. 2012; Liu et al. 2015) and TT disks (e.g., Woitke et al. 2016). Because these previous works used parameterized scale heights and may not have extensively explored disk sizes down to very small radii (e.g., Harvey et al. 2012) explores 50–200 au), we explore the effect here. Figure 10 shows the dependencies of the SED on the three disk parameters investigated in this work ($R_{\text{out}}$, $R_{\text{in}}$, and $M_{\text{dust}}$). At wavelengths shorter than 10 µm the outer radius has some impact on the flux density, and beyond 10 µm it has significant influence on the flux density and the location of the Rayleigh limit. This
<figure><img src="image_13.png" /><figcaption>**Figure 10.** In each panel we show the effect of varying the parameter ($R_{out}$, $R_{in}$, and $M_{dust}$ respectively) while holding all other parameters to the value of the fiducial model (GM Tau; shown in red). Highlighted in the $R_{out}$ plot are models with outer radii of 10 au (orange), 58 au (green), and 105 au (blue).</figcaption></figure>

example demonstrates that the outer disk radius is an important parameter to consider in SED fitting of VLMO disks.

## References

Ahmic, M., Jayawardhana, R., Brandeker, A., et al. 2007, ApJ, 671, 2074

Andrae, R. 2010, arXiv:1009.2755

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241

Apai, D., Pascucci, I., Bouwman, J., et al. 2005, Sci, 310, 834

Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Barlow, R. J. 1989, Statistics: A Guide to the Use of Statistical Methods in the Physical Sciences, Vol. 29 (New York: Wiley)

Bate, M. R. 2012, MNRAS, 419, 3115

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Billot, N., Morales-Calderón, M., Stauffer, J. R., Megeath, S. T., & Whitney, B. 2012, ApJL, 753, L35

Bulger, J., Patience, J., Ward-Duong, K., et al. 2014, A&A, 570, A29

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, 2246

Daemgen, S., Natta, A., Scholz, A., et al. 2016, A&A, 594, A83

Davies, C. L., Gregory, S. G., & Greaves, J. S. 2014, MNRAS, 444, 1157

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

Dent, W. R. F., Thi, W. F., Kamp, I., et al. 2013, PASP, 125, 477

Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075

Feiden, G. A. 2016, A&A, 593, A99

Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Natur, 533, 221

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Natur, 542, 456

Greenwood, A. J., Kamp, I., Waters, L. B. F. M., et al. 2017, arXiv:1702.04744

Guieu, S., Dougados, C., Monin, J.-L., Magnier, E., & Martín, E. L. 2006, A&A, 446, 485

Guieu, S., Pinte, C., Monin, J.-L., et al. 2007, A&A, 465, 855

Guilloteau, S., Dutrey, A., Piétu, V., & Boehler, Y. 2011, A&A, 529, A105

Harvey, P. M., Henning, T., Ménard, F., et al. 2012, ApJL, 744, L1

Harvey, P. M., Jaffe, D. T., Allers, K., & Liu, M. 2010, ApJ, 720, 1374

Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594

Howard, C. D., Sandell, G., Vacca, W. D., et al. 2013, ApJ, 776, 21

Isella, A., Carpenter, J. M., & Sargent, A. I. 2009, ApJ, 701, 260

Kamp, I., Woitke, P., Pinte, C., et al. 2011, A&A, 532, A85

Keane, J. T., Pascucci, I., Espaillat, C., et al. 2014, ApJ, 787, 153

Kelly, B. C. 2007, ApJ, 665, 1489

Klein, R., Apai, D., Pascucci, I., Henning, T., & Waters, L. B. F. M. 2003, ApJL, 593, L57

Liu, Y., Joergens, V., Bayo, A., Nielbock, M., & Wang, H. 2015, A&A, 582, A22

Luhman, K. L. 2000, ApJ, 544, 1044

Luhman, K. L. 2008, in Handbook of Star Forming Regions, Volume II: The Southern Sky ASP Monograph Publications, Vol. 5, ed. B. Reipurth (San Francisco, CA: ASP), 169

Luhman, K. L., Adame, L., D'Alessio, P., et al. 2007, ApJ, 666, 1219

Luhman, K. L., Allen, L. E., Allen, P. R., et al. 2008, ApJ, 675, 1375

Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111

Luhman, K. L., Briceño, C., Stauffer, J. R., et al. 2003, ApJ, 590, 348

Luhman, K. L., Mamajek, E. E., Allen, P. R., & Cruz, K. L. 2009, ApJ, 703, 399

Luhman, K. L., Mamajek, E. E., Shukla, S. J., & Loutrel, N. P. 2017, AJ, 153, 46

Martín, E. L., Dougados, C., Magnier, E., et al. 2001, ApJL, 561, L195

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212

Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155

Mohanty, S., Greaves, J., Mortlock, D., et al. 2013, ApJ, 773, 168

Monin, J.-L., Whelan, E. T., Lefloch, B., Dougados, C., & Alves de Oliveira, C. 2013, A&A, 551, L1

Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012, A&A, 541, A97

Mulders, G. D., & Dominik, C. 2012, A&A, 539, A9

Ott, S. 2010, in ASP Conf. Ser. 434, Astronomical Data Analysis Software and Systems XIX, ed. Y. Mizumoto, K.-I. Morita, & M. Ohishi (San Francisco, CA: ASP), 139

Pascucci, I., Apai, D., Luhman, K., et al. 2009, ApJ, 696, 143

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Pinilla, P., Birnstiel, T., Benisty, M., et al. 2013, A&A, 554, A95

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606

Rebull, L. M., Padgett, D. L., McCabe, C.-E., et al. 2010, ApJS, 186, 259

Ricci, L., Isella, A., Carpenter, J. M., & Testi, L. 2013, ApJL, 764, L27

Ricci, L., Testi, L., Natta, A., et al. 2014, ApJ, 791, 20

Riviere-Marichalar, P., Merín, B., Kamp, I., Eiroa, C., & Montesinos, B. 2016, A&A, 594, A59

Sawicki, M. 2012, PASP, 124, 1208

Schaefer, G. H., Dutrey, A., Guilloteau, S., Simon, M., & White, R. J. 2009, ApJ, 701, 698

Stassun, K. G., Feiden, G. A., & Torres, G. 2014, NewAR, 60, 1

Szűcs, L., Apai, D., Pascucci, I., & Dullemond, C. P. 2010, ApJ, 720, 1668

Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111

Torres, R. M., Loinard, L., Mioduszewski, A. J., & Rodríguez, L. F. 2007, ApJ, 671, 1813

Torres, R. M., Loinard, L., Mioduszewski, A. J., & Rodríguez, L. F. 2009, ApJ, 698, 242

van der Plas, G., Ménard, F., Ward-Duong, K., et al. 2016, ApJ, 819, 102

Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Woitke, P., Pinte, C., Tilling, I., et al. 2010, MNRAS, 405, L26

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868