# Disk Masses and Dust Evolution of Protoplanetary Disks Around Brown Dwarfs

Anneliese M. Rilinger Department of Astronomy and Institute for Astrophysical
Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215
Catherine C. Espaillat Department of Astronomy and Institute for Astrophysical
Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215
Anneliese M. Rilinger<EMAIL_ADDRESS>
(Accepted June 28, 2021 to ApJ)

###### Abstract

We present the largest sample of brown dwarf (BD) protoplanetary disk spectral
energy distributions modeled to date. We compile 49 objects with ALMA
observations from four star-forming regions: $\rho$ Ophiuchus, Taurus, Lupus,
and Upper Scorpius. Studying multiple regions with various ages enables us to
probe disk evolution over time. Specifically, from our models we obtain values
for dust grain sizes, dust settling, and disk mass; we compare how each of
these parameters vary between the regions. We find that disk mass decreases
with age. We also find evidence of disk evolution (i.e., grain growth and
significant dust settling) in all four regions, indicating that planet
formation and disk evolution may begin to occur at earlier stages. We
generally find these disks contain too little mass to form planetary
companions, though we cannot rule out that planet formation may have already
occurred. Finally, we examine the disk mass – host mass relationship and find
that BD disks are largely consistent with previously-determined relationships
for disks around T Tauri stars.

brown dwarfs - planets and satellites: formation - protoplanetary disks -
stars: formation

##  1 Introduction

Young, newly-formed stars are commonly surrounded by disks of gas and dust,
from which planetary companions may form (Williams & Cieza, 2011). Brown
dwarfs (BDs) have also been observed to host these protoplanetary disks
(Comeron et al., 1998; Natta & Testi, 2001; Muench et al., 2001; Natta et al.,
2002). Though usually smaller in size (Hendler et al., 2017), BD
protoplanetary disks are similar to the disks around their young stellar
counterparts (e.g., T Tauri stars, TTS) in many ways. Both BD disks and TTS
disks are typically flat (Scholz et al., 2007) with evidence of dust grain
growth in their inner regions (Sterzik et al., 2004; Apai et al., 2005; Meru
et al., 2013; Pinilla et al., 2013). Some BD disks also show signs of inner
disk clearing (Muzerolle et al., 2006; Rilinger et al., 2019). Furthermore,
the fraction of BDs observed to host disks is similar to the TTS disk
fraction, possibly indicating a shared formation mechanism (Luhman et al.,
2005).

Far-infrared surveys of BDs with the Herschel space observatory (e.g., Harvey
et al., 2012; Liu et al., 2015; Daemgen et al., 2016) greatly enhanced our
understanding of BD disks. Including Herschel far-IR and sub-millimeter
photometry in BD SEDs enabled better constraints on disk parameters via
radiative transfer modeling. Such models show similar disk structure and
geometry to T Tauri disks, including disk scale height (Liu et al., 2015), and
dust settling, though BD disks may be even more settled than their TTS
counterparts (Daemgen et al., 2016). BD disk masses calculated from far-IR and
sub-mm fluxes are consistently lower than typical TTS disk masses, which is
consistent with the positive correlation between disk mass and host mass
reported by, e.g., Andrews et al. (2013) and Pascucci et al. (2016). However,
as noted by Daemgen et al. (2016), the far-IR flux mostly traces the shape and
temperature of the disk, while (sub)-mm fluxes are more sensitive to the disk
mass.

Fortunately, comprehensive millimeter-wavelength surveys of protoplanetary
disks have been made possible with ALMA. Such surveys have been performed in
various star-forming regions for both TTS hosts (e.g., Andrews et al., 2013;
Ansdell et al., 2016; Barenfeld et al., 2016; Eisner et al., 2016; Pascucci et
al., 2016; Ansdell et al., 2017; Eisner et al., 2018; Ruíz-Rodríguez et al.,
2018; Cazzoletti et al., 2019; Cieza et al., 2019; van Terwisga et al., 2019,
2020; Grant et al., 2021) as well as BD hosts in Ophiuchus, (Testi et al.,
2016), Upper Scorpius (van der Plas et al., 2016), Taurus (Ward-Duong et al.,
2018), and Lupus (Sanchis et al., 2020). Though comparisons have been made
between the TTS protoplanetary disks in the different regions, no study has
yet been performed to assess how the BD disks vary across regions.

Disk mass, the amount of gas and dust in the disk, is one of the most
important disk parameters to measure. Protoplanetary disks serve as the mass
reservoir for planet formation and are thought to be the sites where planet
formation occurs (Apai et al., 2005); thus, constraining the amount of
material in the disk directly informs the potential for planet formation. BDs
have been observed to host planetary-mass candidates (Han et al., 2013;
Udalski et al., 2015; Shvartzvald et al., 2017). Since planets are known to be
common around M dwarf stars (Dressing & Charbonneau, 2013; Bonfils et al.,
2013), it follows that planets may also be common around BDs. The formation
mechanism of the BDs themselves, however, may affect their ability to form
planetary companions. BD formation mechanisms can be categorized into two
groups: non-ejection and ejection. In the non-ejection scenario, no
significant disruptions affect the condensation of the BD from its mass
reservoir (Hennebelle & Chabrier, 2009; André et al., 2012; Riaz et al.,
2018). On the other hand, an impulsive reaction with one or more nearby
objects may abruptly eject the BD out of its mass reservoir (Reipurth &
Clarke, 2001; Bate et al., 2003; Basu & Vorobyov, 2012; Stamatellos &
Whitworth, 2009). In these latter ejection scenarios, the disk around the BD
is expected to be truncated or even stripped completely, reducing the amount
of material available for planet formation. Simulations indicate that a BD
disk requires approximately one Jupiter mass of material in order to form a
planet on the order of a few Earth masses; the mass of potential planets and
their likelihood of forming decrease as the mass of the disk decreases (Payne
& Lodato, 2007).

As disks evolve, dust grains are expected to grow and settle: the first steps
towards forming planetesimals and eventually planets. Small ($\sim$0.25
$\mu$m) amorphous silicate dust grains originating in the ISM collide with one
another and can grow to $\sim$micron sizes if the collisions are inelastic. As
grains grow larger, they become more susceptible to the drag force they
experience as they move through the gas in the disk; the dust grains thus
experience the pull of gravity toward the disk midplane and settle out of the
disk atmosphere (Dullemond & Dominik, 2004). Observational studies have shown
that grain growth and dust settling are already underway by $\sim$1 Myr
(Furlan et al., 2009; Ribas et al., 2017; Grant et al., 2018), though older
disks typically show more significant evidence of this evolution.

In addition to the vertical motion of dust grains settling to the midplane,
dust grains also drift radially inwards due to the sub-Keplerian motion of the
gas and the subsequent “headwind” experienced by the dust (Whipple, 1972;
Weidenschilling, 1977). This radial drift is predicted to be more pronounced
for BD disks than TTS since the gas velocity deviates more from the Keplerian
speed as stellar mass decreases (Pinilla et al., 2013, 2017). However, mm-
wavelength observations and measurements of the millimeter spectral index of
BD disks show evidence of large dust grains present in these disks (e.g.,
Ricci et al., 2013, 2014), indicating that some other process prevents radial
drift from depleting the disks of large grains. Possible solutions to the
radial drift problem include pressure bumps within the disk, potentially
created by planetary companions (Pinilla et al., 2012, 2014; Espaillat et al.,
2014; Owen, 2016), lower gas masses than expected within the disks (Pinilla et
al., 2017), or growth of fluffy aggregates, which would remain coupled to the
gas (Kataoka et al., 2013). Evidence of structure (e.g., rings and gaps)
within BD disks could lend support to the pressure bump theory. On the other
hand, fluffy aggregates may have similar opacities to compact dust grains at
mm wavelengths and may be difficult to distinguish from mm-sized dust (Kataoka
et al., 2014).

In this work, we compile the largest sample of BD disks modeled to date, drawn
from the aforementioned studies, in order to search for similarities and
differences between BD protoplanetary disks in various star-forming regions.
We focus on four regions: Ophiuchus (0.5-2 Myr, $\sim$140 pc, Wilking et al.,
2008; Ortiz-León et al., 2017), Taurus (1-2 Myr, $\sim$140 pc, Kenyon &
Hartmann, 1995; Luhman et al., 2010; Andrews et al., 2013), Lupus (1-3 Myr,
150-200 pc, Comerón, 2008; Ansdell et al., 2016), and Upper Scorpius (5-11
Myr, $\sim$145 pc, Preibisch et al., 2002; Pecaut et al., 2012). In
particular, we compare disk masses and disk evolution (characterized by grain
growth and dust settling) in these regions. By considering regions of varying
age, we can probe how BD disks change over their lifetimes.

In Section 2 we present our BD sample selection. Our modeling analysis is
described in Section 3. We discuss the results of our analysis in Section 4,
and summarize our findings in Section 5.

##  2 Sample

We have compiled a sample of 49 BDs with protoplanetary disks located in the
$\rho$ Ophiuchus, Taurus, Lupus, and Upper Scorpius star-forming regions. This
constitutes the largest sample of BD disk spectral energy distributions (SEDs)
modeled to date. These four regions were selected because the BDs in these
regions have been observed in the (sub-)millimeter: Ophiuchus (Testi et al.,
2016, and Testi et al., submitted), Taurus (Andrews et al., 2013; Ward-Duong
et al., 2018), Lupus (Sanchis et al., 2020), and Upper Sco (van der Plas et
al., 2016). We only consider objects with masses less than the hydrogen-
burning limit ($\sim$0.08 M⊙) and/or spectral types of M6 or later since these
objects are most likely to be bona fide BDs (Luhman et al., 2005). We also
required the objects to have millimeter-wavelength photometry measurements or
upper limits, in order to be able to better constrain our SED models.
Additionally, we did not include objects in our sample which are known to be
in binary or multiple systems, since companions can affect disk properties.
Finally, one object, 2MASS J04381486+2611399, has been reported to be a
highly-inclined, edge-on disk (Luhman et al., 2007) so we remove this object
from our sample. Table 1 lists the 49 objects in our final sample, along with
their spectral types, temperatures, extinctions, luminosities, masses, radii,
and distances, compiled from the literature (see references in Table 1). We
note that the relevant properties have been scaled to Gaia DR2 distances for
all objects that had Gaia DR2 parallaxes (Gaia Collaboration et al., 2016,
2018).

Table 1: Parameters of BDs modeled in this work

Object | SpT | T∗ | AV | L∗ | M∗ | R∗ | Distance | Source
---|---|---|---|---|---|---|---|---
|  | (K) |  | log(L⊙) | (M⊙) | (R⊙) | (pc) |
Members of Taurus
J04390396+2544264 | M7.25 | 2837 | 0.5 | -1.336 | 0.035 | 0.552 | 143.99${}^{+4.50}_{-4.24}$ | 1
J04141188+2811535 | M6.25 | 2963 | 2.5 | -1.746 | 0.053 | 0.873 | 131.09${}^{+2.92}_{-2.79}$ | 1
J04292165+2701259 | M6 | 3091 | 2 | -0.115 | 0.058 | 1.884 | 140aaFrom reference listed in Source column. | 1
J04230607+2801194 | M6 | 2990 | 1.5 | -1.332 | 0.058 | 0.942 | 133.88${}^{+2.49}_{-2.40}$ | 1
KPNO Tau 3 | M6 | 2990 | 1.6 | -1.655 | 0.058 | 0.377 | 155.88${}^{+5.72}_{-5.33}$ | 1
GM Tau | M6.5 | 2935 | 0.6 | -1.507 | 0.048 | 0.794 | 138.31${}^{+2.91}_{-2.79}$ | 1
J04390163+2336029 | M6 | 2990 | 0.5 | -1.054 | 0.058 | 1.13 | 127.81${}^{+1.31}_{-1.29}$ | 1
J04400067+2358211 | M6 | 2858 | 0.5 | -1.566 | 0.058 | 0.67 | 120.49${}^{+2.35}_{-2.26}$ | 1
J04414825+2534304 | M7.75 | 2752 | 1.3 | -1.683 | 0.028 | 0.634 | 136.16${}^{+3.84}_{-3.64}$ | 1
J04442713+2512164 | M7.25 | 2838 | 0.0 | -1.553 | 0.05 | 0.69 | 141.0$\pm$2.7 | 2
CFHT Tau 4 | M7 | 2880 | 5.67 | -0.757 | 0.095 | 1.68 | 147.1$\pm$5.1 | 2
CFHT Tau 9 | M5.75 | 3023 | 2.14 | -1.517 | 0.064 | 0.636 | 155.51${}^{+3.21}_{-3.08}$ | 3
KPNO Tau 6 | M8.5 | 2555 | 1.28 | -2.468 | 0.02 | 0.298 | 116.05${}^{+7.65}_{-6.76}$ | 3
KPNO Tau 7 | M8.25 | 2632 | 1.57 | -2.272 | 0.024 | 0.352 | 122.82${}^{+7.35}_{-6.57}$ | 3
J04330945+2246487 | M6 | 2990 | 4.53 | -1.31 | 0.069 | 0.825 | 148.98${}^{+10.06}_{-8.87}$ | 3
J04214631+2659296 | M5.75 | 3023 | 4.72 | -1.529 | 0.063 | 0.627 | 160.46${}^{+7.63}_{-6.97}$ | 3
J04414489+2301513 | M8.5 | 2555 | 1.4 | -2.262 | 0.023 | 0.377 | 120.42${}^{+5.74}_{-5.24}$ | 3
KPNO Tau 12 | M9 | 2400 | 1.88 | -2.973 | 0.011 | 0.189 | 140aaFrom reference listed in Source column. | 3
J04201611+2821325 | M6.5 | 2935 | 2.3 | -1.87 | 0.041 | 0.449 | 128.99${}^{+3.43}_{-3.25}$ | 3
J04263055+2443558 | M8.75 | 2478 | 0.9 | -2.26 | 0.02 | 0.402 | 119.88$\pm 10.05$ | 3
J04290068+2755033 | M8.25 | 2700 | 1.71 | -1.98 | 0.023 | 0.467 | 145.65${}^{+7.85}_{-7.08}$ | 3
J04554969+3019400 | M6 | 2990 | 0.09 | -1.695 | 0.052 | 0.529 | 156.51${}^{+4.76}_{-4.49}$ | 3
J04574903+3015195 | M9.25 | 2350 | 0.2 | -2.96 | 0.015 | 0.2 | 140aaFrom reference listed in Source column. | 4
Members of Ophiuchus
CFHTWIR-Oph 16 | M8 | 2710 | 11.71 | -2.00 | 0.03 | 0.45 | 139.4aaFrom reference listed in Source column. | 6
ISO-Oph 23 | M6.5 | 2930 | 7.91 | -1.40 | 0.07 | 0.69 | 139.4aaFrom reference listed in Source column. | 5
ISO-Oph 32 | M7.25 | 2840 | 2.09 | -1.18 | 0.08 | 1.06 | 151.26${}^{+4.88}_{-4.58}$ | 5
ISO-Oph 33 | M87 | 2880 | 6.82 | -2.23 | 0.06 | 0.31 | 139.4aaFrom reference listed in Source column. | 5
GY92 141 | M8.5 | 2550 | 1.16 | -2.54 | 0.02 | 0.28 | 142.7aaFrom reference listed in Source column. | 6
CFHTWIR-Oph 58 | M8 | 2710 | 6.20 | -2.7 | 0.04 | 0.2 | 139.4aaFrom reference listed in Source column. | 6
CFHTWIR-Oph 66 | M7.75 | 2750 | 10.62 | -2.65 | 0.05 | 0.21 | 139.4aaFrom reference listed in Source column. | 6
CFHTWIR-Oph 77 | M9.75 | 2260 | 7.29 | -3.02 | 0.02 | 0.2 | 139.4aaFrom reference listed in Source column. | 6
GY92 264 | M8 | 2710 | 0.31 | -1.74 | 0.04 | 0.61 | 139.15${}^{+2.88}_{-2.77}$ | 5
CFHTWIR-Oph 90 | L0 | 2250 | 1.94 | -2.65 | 0.01 | 0.31 | 139.4aaFrom reference listed in Source column. | 6
GY92 320 | M7.5 | 2753 | 2.09 | -2.15 | 0.04 | 0.37 | 140.03${}^{+15.11}_{-12.43}$ | 5
CFHTWIR-Oph 98 | M9.75 | 2260 | 3.72 | -2.96 | 0.01 | 0.22 | 139.4aaFrom reference listed in Source column. | 6
CFHTWIR-Oph 107 | M6.25 | 2960 | 1.71 | -2.02 | 0.07 | 0.37 | 159.6aaFrom reference listed in Source column. | 6
Members of Lupus
J15451851-3421246 | M6.5 | 2935 | 0 | -1.40 | 0.09 | 0.774 | 152$\pm$4 | 7
SONYC Lup 3-7 | M8.5 | 2600 | 0 | -2.0 | 0.02 | 0.493 | 150$\pm$6 | 7
Lup 706 | M7.5 | 2795 | 0 | -2.70 | 0.05 | 0.191 | 158.5aaFrom reference listed in Source column. | 7
AKC 2006-18 | M6.5 | 2935 | 0 | -2.0 | 0.07 | 0.387 | 149$\pm$8 | 7
Lup 818s | M6 | 2990 | 0 | -1.70 | 0.09 | 0.527 | 157$\pm$3 | 7
J16101984-3836065 | M6.5 | 2935 | 0 | -1.40 | 0.09 | 0.774 | 159$\pm$3 | 7
J16085529-3848481 | M6.5 | 2935 | 0.5 | -1.30 | 0.09 | 0.865 | 158$\pm$3 | 7
Lup 607 | M6.5 | 2935 | 0 | -1.30 | 0.1 | 0.865 | 175$\pm$6 | 7
Members of Upper Scorpius
J15555600-2045187 | M6.5 | 3400 | 3.38 | -1.56 | 0.07 | 0.479 | 145.65${}^{+4.66}_{-4.38}$ | 8
J15560104-2338081 | M6.5 | 2800 | 1.3 | -1.99 | 0.07 | 0.432 | 139.88${}^{+5.65}_{-5.23}$ | 8
J15591135-2338002 | M7 | 2500 | 1.3 | -2.20 | 0.06 | 0.421 | 139.94${}^{+7.39}_{-6.69}$ | 8
J16100541-1919362 | M7 | 2600 | 0.7 | -2.19 | 0.06 | 0.397 | 149.69${}^{+7.40}_{-6.73}$ | 8
J16060391-2056443 | M7.5 | 2800 | 1.5 | -1.83 | 0.04 | 0.518 | 137.15${}^{+5.22}_{-4.85}$ | 8

The Source column gives the reference for all other properties listed here:
(1) Ward-Duong et al. (2018); (2) Ricci et al. (2014) and references therein;
(3) Andrews et al. (2013); (4) Luhman (2004); (5) Testi et al. (2016) and
Testi et al., submitted; (6) Testi et al., submitted; (7) Sanchis et al.
(2020); (8) mass and spectral type from van der Plas et al. (2016), other
properties from Liu et al. (2015).

Note. — Distances were computed from Gaia parallaxes (Gaia Collaboration et
al., 2016, 2018) unless otherwise noted. These distances were used to scale
the relevant properties obtained from the literature.

##  3 Analysis and Results

###  3.1 SED Modeling

The disk properties for each object were determined by fitting disk structure
models to the SED of each BD. To construct the SEDs, photometric data points
from visible to millimeter wavelengths were taken from the Vizier catalogue
access tool (Ochsenbein et al., 2000). Of the 49 objects in this sample, 31
also have low resolution (R $\sim$ 60–130) spectra from the InfraRed
Spectrograph (IRS) on the Spitzer Space Telescope (Houck et al., 2004). The
SEDs for each object are presented in Appendix A, Figures A.1, A.3, A.4, and
A.5.

We use the D’Alessio Irradiated Accretion Disk (DIAD) radiative transfer
models (D’Alessio et al., 1998, 1999, 2001, 2005, 2006) to model the vertical
and radial structure of the disks. Previous works (e.g., Morrow et al., 2008;
Adame et al., 2011; Rilinger et al., 2019) have used these models to fit BD
SEDs. The disks are assumed to be irradiated by the central BD. DIAD
calculates the vertical and radial structure of each disk self-consistently
while enforcing hydrostatic equilibrium. Two populations of dust grains are
modeled: a population of small dust grains in the atmosphere of the disk and
another consisting of large grains in the disk midplane. In each dust
population, the minimum grain size is fixed at 0.005 $\mu$m; particle sizes
are distributed following a power law with a power of –3.5 (Mathis et al.,
1977) up to maximum sizes that are left as free parameters in the model
(amax,atm and amax,mid). The disk mass is calculated by integrating over the
surface mass density ($\Sigma$), which in turn is determined by the mass
accretion rate ($\dot{M}$) and disk viscosity ($\alpha$;
$\Sigma\propto\dot{M}\alpha^{-1}$). An important feature of the DIAD models is
that $\Sigma$ is calculated self-consistently with the vertical and radial
structure of the disk, as opposed to assuming a fixed power-law surface
density distribution. Of the 49 objects in this sample, 18 have measured
$\dot{M}$ values in the literature (see Table 2). We fix $\dot{M}$ for these
18 objects; for the remaining objects we adopt a value of $1\times 10^{-10}$
${\rm{\it M}_{\odot}\,yr^{-1}}$, which is approximately the average of the 18
measured $\dot{M}$ values. Disk viscosity ($\alpha$) is left as a free
parameter in the DIAD model.

Table 2: BDs with literature accretion rates

Object | $\dot{M}$ | Source
---|---|---
| (M⊙ yr-1) |
J04390396+2544264 | 6.9 x 10-12 | 1
J04141188+2811535 | 9.0 x 10-10 | 1
GM Tau | 1.9 x 10-9 | 1
J04414825+2534304 | 1.9 x 10-11 | 1
J04442713+2512164 | 2.0 x 10-10 | 2
KPNO Tau 6 | 3.6 x 10-11 | 1
KPNO Tau 12 | 1.6 x 10-11 | 1
ISO-Oph 23 | 1.4 x 10-10 | 3
ISO-Oph 32 | 3.2 x 10-11 | 3
ISO-Oph 33 | 2.0 x 10-11 | 3
J15451851-3421246 | 1.77 x 10-11 | 4
SONYC Lup 3-7 | 1.39 x 10-11 | 4
Lup 706 | 1.52 x 10-12 | 4
AKC 2006-18 | 5.48 x 10-12 | 4
Lup 818s | 1.15 x 10-11 | 4
J16101984-3836065 | 2.74 x 10-11 | 4
J16085529-3848481 | 1.91 x 10-11 | 4
Lup 607 | 3.34 x 10-12 | 4

Note. — Sources: (1) Herczeg & Hillenbrand (2008); (2) Bouy et al. (2008); (3)
Manara et al. (2015); (4) Alcalá et al. (2017), scaled to Gaia distances.

The inner edge or “wall” of the disk is set by Twall,in, a temperature that is
a free parameter in the DIAD model. For a full-disk model, Twall,in is set to
the dust destruction temperature, which we assume to be 1400 K. To model a
transitional disk (i.e., a disk with a hole), this temperature can be
decreased, effectively moving the inner edge of the disk further from the
central BD. We can also model a pre-transitional disk (i.e., a disk with a
large gap separating the inner and outer disk) by modifying the disk structure
to include an annular gap within the disk (Espaillat et al., 2011); in this
case, the disk has both an inner wall and an outer wall, set by temperature
Twall,out. Other free parameters of the DIAD model include the inner wall
scale height (Hwall,in), disk outer radius (Rout), inclination (i), and dust
settling ($\epsilon$, see Section 3.1.3). If the disk is a pre-transitional
disk, the scale height of the outer wall, Hwall,out, is also a free parameter
in the model. We assume a fixed dust-to-gas mass ratio of 0.01.

In Figure 1, we explore the effect each of these disk parameters on the DIAD
model. All but one parameter is held constant in each panel to show how that
parameter changes the shape of the resulting SED. In all cases, we fix M∗, T∗,
L∗, and distance to the average values from Table 1. We also fix inclination
to 60∘, $\dot{M}$ to 1x10-10 M⊙ yr-1, and the dust-to-gas mass ratio to 0.01.
The remaining default parameter values are as follows: $\alpha$ = 0.001,
$\epsilon$ = 0.01, Rout = 50 au, amax,atm = 1.0 $\mu$m, amax,mid = 500 mm,
Twall = 1400 K, and Hwall = 1. Figure 1 shows the importance of including mm
photometry in the SEDs: one or two mm-wavelength photometry points can help
break the degeneracy between the $\alpha$ and $\epsilon$ parameters in the
far-IR, resulting in more accurate measurements of the vertical structure
(settling) as well as the disk mass.

Figure 1: Effect of different parameters in the DIAD SED models. In each
panel, all but one parameter is held constant. See text for the values of the
fixed parameters.

Photospheres are constructed for each BD using the Kenyon & Hartmann (1995)
color table and the reported T∗, L∗, and R∗ values in Table 1. Colors are
scaled to the dereddened observed J-band magnitude of the BD which is in turn
calculated from the L∗ value in Table 1. We then interpolate to determine the
photospheric emission at each wavelength. Many objects in the sample have
spectral types later than M6 (see Table 1), the latest spectral type in the
Kenyon & Hartmann (1995) color table. A spectral type of M6 was thus adopted
for these objects. Note that while the Pecaut & Mamajek (2013) color table
includes later spectral types, spectral types later than M5 are missing
infrared wavelengths beyond K-band. Since the two color tables show good
agreement (within 10% uncertainty) between the Pecaut & Mamajek (2013)
spectral type M7 colors and Kenyon & Hartmann (1995) spectral type M6 colors,
we opted to use the more complete Kenyon & Hartmann (1995) color table to
construct the BD photospheres.

DIAD disk models are computationally expensive to run: each model typically
takes 1.5–2 hours to complete. Therefore, it is unworkable to employ common
statistical methods of estimating best-fit model parameters and their
uncertainties (e.g., Markov-Chain Monte Carlo methods, Levenberg-Marquardt
$\chi^{2}$ minimization, etc.). To obtain best-fit SED models, we first ran a
coarse grid of models for each object to explore the parameter space (see
Table 3 for the initial grid parameters). We then iteratively refined the
grids as necessary based on visually inspecting the models with minimum
$\chi^{2}$. Upper limit points were not included in the $\chi^{2}$
calculation, but rather served as further visual guides for assessing the
goodness of the fit (i.e., we do not accept a model that is inconsistent with
the upper limits, even if it fits the other points well). The reported values
for the best-fit parameters thus represent qualitative approximations to the
range of possible values of the parameters rather than accurate quantitative
values.

Table 3: Initial Grid Parameters Parameter | Range | Units
---|---|---
$\epsilon$ | {0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1} | …
$\alpha$ | {0.0001, 0.0005, 0.001, 0.005, 0.01} | …
amax,atm | {0.25, 1.0, 10.0} |  $\mu$m
amax,mid | {0.5, 1.0, 5.0, 10.0} | mm
Hin | {1, 2, 3, 4, 5} | …
Twall,in | 400–1400, in steps of 100 | K

In addition to the parameters in Table 3, which were varied, three additional
parameters were fixed for each object: $\dot{M}$, inclination, and Rout. We
fix $\dot{M}$ as described above. Inclination was fixed at 60o except for
J04442713 and CFHT Tau 4, whose inclinations were determined by Rilinger et
al. (2019). Disk radii were adopted from the literature where possible. We
adopt the Rout values determined via ALMA modeling by Rilinger et al. (2019)
for J04442713 and CFHT Tau 4. Testi et al. (2016) constrained the disk radii
for the 5 BDs in their sample: the Rout values for their BD disks are $\leq$
25 au. We adopt these Rout values for those objects, and assume the same value
for the additional 8 objects observed by Testi et al., submitted. Sanchis et
al. (2020) report disk radii for five of the disks in Lupus; we adopt these
Rout values rounded to the nearest ten. For the remaining Lupus disks, we
adopt the average of the $R_{out}$ values for this region, 50 au. van der Plas
et al. (2016) did not resolve any of the disks in their Upper Sco sample. An
Rout of 50 au is consistent with their angular resolution and the distance to
these objects, and was adopted for consistency with the Lupus disks.
Similarly, Ward-Duong et al. (2018) do not measure disk radii for the objects
in their Taurus sample. Adopting Rout values of 50 au for these objects is
consistent with our assumptions for other regions. An Rout of 50 au was
assumed for the remaining disks in the sample, with two additional exceptions:
the SEDs of J04330945 and J04263055 appear to show truncated disks. Rout
values of 0.1, 0.2, 0.3, 0.4, 0.5, 1, 5, 10, and 15 au were tested for these
two objects. The best-fit SED models for each object are shown in Appendix A,
Figures A.1, A.3, A.4, and A.5. The best-fit parameters for each of these
models are given in Table 4.

Table 4: Best-fit parameters for BD SED models

Object | $\epsilon$ | $\alpha$ | amax,atm | amax,mid | zwall,in | Twall,in | Rwall,inaaThis is not a free parameter and was calculated using Twall,in following D’Alessio et al. (2005). | Rout bbConstrained by literature values where possible as described in Section 3.1. | i ccFixed to 60∘ except for two objects with previous ALMA modeling. See Section 3.1. | MdiskddThis is not a free parameter and was
---|---|---|---|---|---|---|---|---|---|---
|  |  | ($\mu$m) | (mm) | (au) | (K) | (au) | (au) | (∘) | (MJup)
Members of Taurus
J04390396+2544264 | 0.01 | 6 x 10-5 | 10.0 | 0.5 | 0.0212 | 500 | 0.1 | 50 | 60 | 0.283
J04141188+2811535 | 0.0005 | 0.01 | 1.0 | 0.5 | 0.0531 | 400 | 0.25 | 50 | 60 | 0.188
J04292165+2701259 | 0.05 | 0.0005 | 0.25 | 0.5 | 0.0079 | 1400 | 0.06 | 50 | 60 | 0.366
J04230607+2801194 | 0.0005 | 0.0001 | 10.0 | 10.0 | 0.0011 | 1400 | 0.02 | 50 | 60 | 2.964
KPNO Tau 3 | 0.001 | 0.0001 | 10.0 | 10.0 | 0.0016 | 1400 | 0.01 | 50 | 60 | 2.964
GM Tau | 0.005 | 0.0005 | 10.0 | 5.0 | 0.0377 | 600 | 0.11 | 50 | 60 | 0.450
J04390163+2336029 | 0.0001 | 0.001 | 0.25 | 0.5 | 0.0016 | 1400 | 0.03 | 50 | 60 | 0.209
J04400067+2358211 | 0.001 | 0.0001 | 10.0 | 10.0 | 0.0005 | 1400 | 0.01 | 50 | 60 | 3.121
J04414825+2534304 | 0.01 | 7 x 10-5 | 0.25 | 0.5 | 0.0007 | 1400 | 0.01 | 50 | 60 | 0.628
J04442713+2512164eeBest fit as a pre-transitional disk. See Table 5 for additional parameters.ffParameters from previous modeling by Rilinger et al. (2019). | 0.05 | 0.0005 | 1.0 | 5.0 | 0.0083 | 1400 | 0.07 | 100 | 40 | 2.050
CFHT Tau 4ffParameters from previous modeling by Rilinger et al. (2019). | 0.007 | 0.0007 | 0.25 | 3.0 | 0.0020 | 1400 | 0.07 | 80 | 70 | 0.420
CFHT Tau 9eeBest fit as a pre-transitional disk. See Table 5 for additional parameters. | 0.005 | 0.005 | 10.0 | 5.0 | 0.0006 | 1400 | 0.02 | 50 | 60 | 0.052
KPNO Tau 6 | 0.0005 | 0.0005 | 10.0 | 1.0 | 0.0004 | 1400 | 0.01 | 50 | 60 | 0.262
KPNO Tau 7 | 0.0005 | 0.0005 | 10.0 | 0.5 | 0.0008 | 1400 | 0.01 | 50 | 60 | 0.524
J04330945+2246487 | 0.0001 | 0.001 | 1.0 | 0.5 | 0.0010 | 1400 | 0.02 | 0.2 | 60 | 0.002
J04214631+2659296 | 0.01 | 0.005 | 10.0 | 10.0 | 0.0011 | 1400 | 0.02 | 50 | 60 | 0.052
J04414489+2301513 | 0.005 | 0.001 | 10.0 | 0.5 | 0.0003 | 1400 | 0.01 | 50 | 60 | 0.220
KPNO Tau 12 | 0.001 | 0.0005 | 1.0 | 0.5 | 0.0005 | 1400 | 0.003 | 50 | 60 | 0.095
J04201611+2821325 | 0.001 | 0.005 | 1.0 | 0.5 | 0.0005 | 1400 | 0.01 | 50 | 60 | 0.042
J04263055+2443558 | 0.1 | 0.005 | 1.0 | 1.0 | 0.0003 | 1400 | 0.01 | 0.1 | 60 | 1.37 x 10-4
J04290068+2755033 | 0.001 | 0.001 | 1.0 | 1.0 | 0.0005 | 1400 | 0.01 | 50 | 60 | 0.157
J04554969+3019400 | 0.001 | 0.005 | 0.25 | 0.5 | 0.0011 | 1400 | 0.02 | 50 | 60 | 0.042
J04574903+3015195 | 0.1 | 0.0001 | 1.0 | 1.0 | 0.0001 | 1400 | 0.003 | 50 | 60 | 3.048
Members of Ophiuchus
CFHTWIR-Oph 16 | 0.05 | 0.001 | 1.0 | 5.0 | 0.0037 | 1400 | 0.01 | 25 | 60 | 0.157
ISO-Oph 23 | 0.001 | 0.001 | 10.0 | 10.0 | 0.0074 | 800 | 0.05 | 25 | 60 | 0.241
ISO-Oph 32 | 0.1 | 0.0005 | 10.0 | 1.0 | 0.0082 | 800 | 0.08 | 25 | 60 | 0.115
ISO-Oph 33 | 0.1 | 0.0002 | 0.25 | 10.0 | 0.0401 | 500 | 0.09 | 25 | 60 | 0.272
GY92 141 | 0.0005 | 0.001 | 10.0 | 5.0 | 0.0002 | 1400 | 0.004 | 25 | 60 | 0.131
CFHTWIR-Oph 58 | 0.01 | 0.005 | 10.0 | 5.0 | 0.0001 | 1400 | 0.003 | 25 | 60 | 0.036
CFHTWIR-Oph 66 | 0.1 | 0.0001 | 0.25 | 5.0 | 0.0011 | 1400 | 0.01 | 25 | 60 | 3.184
CFHTWIR-Oph 77 | 0.001 | 0.001 | 1.0 | 10.0 | 0.0001 | 1400 | 0.003 | 25 | 60 | 0.114
GY92 264 | 0.005 | 0.0001 | 1.0 | 10.0 | 0.0029 | 1400 | 0.01 | 25 | 60 | 1.56
CFHTWIR-Oph 90 | 0.05 | 0.005 | 10.0 | 1.0 | 0.0002 | 1400 | 0.004 | 25 | 60 | 0.018
GY92 320 | 0.005 | 0.005 | 10.0 | 10.0 | 0.0012 | 1400 | 0.01 | 25 | 60 | 0.031
CFHTWIR-Oph 98 | 0.1 | 0.001 | 10.0 | 10.0 | 0.0006 | 1400 | 0.003 | 25 | 60 | 0.135
CFHTWIR-Oph 107 | 0.0005 | 0.01 | 1.0 | 1.0 | 0.0003 | 1400 | 0.01 | 25 | 60 | 0.021
Members of Lupus
J15451851-3421246 | 0.01 | 5 x 10-5 | 0.25 | 0.5 | 0.0164 | 1000 | 0.05 | 15 | 60 | 0.443
SONYC Lup 3-7 | 0.01 | 0.0001 | 10.0 | 10.0 | 0.0013 | 1400 | 0.01 | 30 | 60 | 0.199
Lup 706 | 0.05 | 5 x 10-5 | 10.0 | 0.5 | 0.0045 | 700 | 0.02 | 95 | 60 | 0.268
AKC 2006-18 | 0.0005 | 0.0005 | 1.0 | 10.0 | 0.0003 | 1400 | 0.01 | 50 | 60 | 0.052
Lup 818s | 0.1 | 1 x 10-5 | 1.0 | 5.0 | 0.0004 | 1400 | 0.01 | 20 | 60 | 2.138
J16101984-3836065 | 0.0005 | 0.0005 | 0.25 | 5.0 | 0.0007 | 1400 | 0.02 | 50 | 60 | 0.157
J16085529-3848481 | 0.0005 | 0.0001 | 10.0 | 5.0 | 0.0007 | 1400 | 0.02 | 40 | 60 | 0.578
Lup 607 | 0.0005 | 0.0001 | 1.0 | 5.0 | 0.0008 | 1400 | 0.02 | 50 | 60 | 0.147
Members of Upper Scorpius
J15555600-2045187 | 0.0005 | 0.001 | 10.0 | 10.0 | 0.0022 | 1400 | 0.02 | 50 | 60 | 0.325
J15560104-2338081 | 0.0005 | 0.0005 | 10.0 | 1.0 | 0.0007 | 1400 | 0.01 | 50 | 60 | 0.733
J15591135-2338002 | 0.001 | 0.0005 | 10.0 | 0.5 | 0.0009 | 1200 | 0.01 | 50 | 60 | 0.721
J16100541-1919362 | 0.01 | 0.005 | 10.0 | 0.5 | 0.0012 | 1200 | 0.01 | 50 | 60 | 0.063
J16060391-2056443 | 0.0005 | 0.005 | 1.0 | 10.0 | 0.0005 | 1400 | 0.01 | 50 | 60 | 0.042

calculated by DIAD. See Section 3.1.1

Note. — $\epsilon$, $\alpha$, amax,atm, amax,mid, Hin, and Twall,in, in were
varied to obtain the best fit starting with the initial grid in Table 3. The
zwall values reported here are the Hin values converted to units of au.

####  3.1.1 Disk Mass

One of the most important disk properties to measure is the disk mass, since
the mass of the disk constrains the potential for planetary companion
formation. As described above, the total (gas and dust) disk mass is
determined by integrating over the surface mass density ($\Sigma$) calculated
by DIAD. For each BD in the sample, the disk mass was computed in this way
using the $\Sigma$, $\alpha$, and $\dot{M}$ of the best-fit SED model. A
histogram of the disk masses found for this sample is shown in Figure 2. The
vast majority of disks in this sample have total disk masses below 1 MJ. Seven
disks (J04230607, KPNO Tau 3, J04400067, J04442713, CFHTWIR-Oph 66, GY92 264,
and Lup 818s) are clear outliers, with more massive disks.

####  3.1.2 Dust Growth

Dust grains in the interstellar medium (ISM) are typically assumed to have a
maximum size of 0.25 $\mu$m (Mathis et al., 1977; Weingartner & Draine, 2001).
The presence of larger ($>$1 $\mu$m) grains thus indicates that dust growth
has taken place within the disk and therefore that the first steps towards
planet formation have occurred.

The 10 $\mu$m silicate feature commonly observed in TTS and BD disks contains
significant information about grain growth in protoplanetary disks. In these
disks, this feature is often observed to be in emission (Cohen & Witteborn,
1985; Gürtler et al., 1999; Natta et al., 2000). The presence of an emission
line indicates that a (super)heated optically thin dust layer exists above the
optically thick disk (Calvet et al., 1991). In other words, the 10 $\mu$m
emission feature is produced by the disk atmosphere. The shape of the feature
corresponds to the size of grains in this atmosphere. Submicron-sized grains
show a more narrow, peaked shape; as grains grow to be larger than $\sim$ 1
$\mu$m, the 10 $\mu$m feature broadens and flattens.

Figure 2: Histogram of the calculated disk masses for the BDs in this sample.
Only a few outliers have masses greater than 1 MJ, indicating that most of the
BD disks are not currently massive enough to form planetary companions
according to simulations by Payne & Lodato (2007).

The maximum size of dust grains in the disk atmosphere, amax,atm, is one of
the parameters varied in our SED models. Furthermore, the wavelength coverage
of IRS spectra encompasses the 10 $\mu$m spectral feature. Therefore, for the
34 objects in our sample with IRS spectra, we can use the model fit to the 10
$\mu$m feature to probe grain growth. Note that we do not attempt to fit the
spectral feature in detail (e.g., Sargent et al., 2009); rather, we consider
its width and height to generally constrain the grain size. We assume the dust
particles to be a standard mixture of olivine silicates and graphite. Relative
to the gas, silicates have a fractional mass abundance of 0.004 and graphite
has a fractional mass abundance of 0.0025. Much work has been done to
understand how different dust compositions and opacities affect the SED (e.g.,
Espaillat et al., 2010, 2012); however varying the species and/or relative
abundances is outside the scope of this work. We assume the standard opacities
calculated from Mie theory; optical constants are taken from Dorschner et al.
(1995) and Draine & Lee (1984) for the silicates and graphite, respectively.

Table 5: Outer wall parameters Object | zwall,out | Twall,out | Rwall,out | Gap width
---|---|---|---|---
| au | K | au | au
J04442713 | 0.0827 | 400 | 0.27 | 0.20
CFHT Tau 9 | 0.0275 | 300 | 0.36 | 0.34

We ran models with three different amax,atm for atmospheric dust grains: 0.25
$\mu$m (i.e., ISM grains), 1.0 $\mu$m (moderate grain growth), and 10.0 $\mu$m
(significant grain growth). Figure 3 shows a histogram of the amax,atm values
from the best-fit models to each BD SED. The majority of disks are best fit
with grains larger than standard ISM grains, indicating that significant grain
growth has taken place in these disks.

Figure 3: Histogram of the best-fit amax,atm values for the BDs in this
sample. The majority of disks in this sample have amax,atm greater than 0.25
$\mu$m, indicating grain growth in these disks.

####  3.1.3 Dust Settling

Another expected process in disk evolution is the settling of dust grains
towards the disk midplane as the grains grow larger. To parameterize the
amount of dust settling, we assume the grain size distributions and the dust-
to-gas mass ratio ($\zeta$) to be constant in radius and vary only with height
above and below the disk midplane. In order to keep the total $\zeta$ constant
while still allowing for settling of large grains, we decrease the dust-to-gas
ratio of small grains ($\zeta_{small}$) in the atmosphere and increase the
dust-to-gas ratio of larger grains ($\zeta_{big}$) in the midplane. We thus
define the settling parameter $\epsilon$ such that

$\epsilon=\zeta_{small}/\zeta_{std}$ (1)

where $\zeta_{std}$ is the standard assumed dust-to-gas mass ratio of 0.01.
Therefore, lower $\epsilon$ values correspond to more settled disks.

In our grid of models, we vary $\epsilon$ between 0.0005 and 0.1. The
histogram of best-fit $\epsilon$ values for the disks in this sample is shown
in Figure 4. The distribution peaks towards lower $\epsilon$ values,
indicating a high degree of settling in BD disks. In conjunction with the
grain growth present in many BD disks modeled here, this high degree of
settling indicates that BD disks show significant evidence of dust evolution.

Figure 4: Distribution of the dust settling parameter $\epsilon$ for all BDs
in this sample.

####  3.1.4 Disk Clearing

Disks that show evidence of clearing (i.e., pre-transitional and transitional
disks; PTDs and TDs, respectively) are not uncommon around around T Tauri
stars (e.g., Espaillat et al., 2007; Muzerolle et al., 2010; Kim et al.,
2013). BD disks may also exhibit similar substructures.

In young stellar objects (YSOs), dust sublimation temperatures typically range
between 1000–2000 K (Monnier & Millan-Gabet, 2002). Therefore, inner wall
temperatures below 1000 K indicate additional disk clearing beyond the dust
sublimation radius expected in the inner disk. In other words, inner wall
temperatures below 1000 K correspond to a larger inner disk radius and thus a
dust-depleted hole in the disk. In our sample of BDs, 7 disks were best fit
with inner wall temperatures below 1000 K and are therefore classified as TDs
(see Table 4, column 7).

Additionally, two objects in this sample were best fit as PTDs, having an
inner disk and outer disk separated by a dust-depleted gap. These objects and
their outer disk parameters are listed in Table 5. One of these PTDs,
J04442713, was previously reported by Rilinger et al. (2019). We find that 9
objects, 18 percent of the disks studied here, show some evidence of disk
clearing. This fraction of disks with evidence of clearing is consistent with
the fraction of TDs that has been reported for TTS (between 10 and 20 percent
of objects, see Muzerolle et al., 2010).

It is important to note that these holes and gaps were inferred from the SEDs
of these objects and have not been directly detected through mm-wavelength
imaging (e.g., with ALMA). Most of the TDs presented here do not have clear
dips in their SEDs as are seen in TTS TDs. Ercolano et al. (2009) note that
the unambiguous detection of TDs around late-M stars is complicated by their
lower stellar luminosities: the ratio between the dust and stellar
contributions to the SED is smaller (i.e., the contrast is lower) than for
earlier-type stars. This challenge in interpreting the SED means that high-
resolution imaging of these objects is necessary to confirm the presence of
these substructures. However, given the small inner hole radii and gap widths
(see Tables 4 and 5), these features may be below current resolution
limitations.

##  4 Discussion

Having obtained models for all 49 objects in our sample, we can compare how
the four star-forming regions differ from one another. Since the four regions
studied here vary in age, differences between regions correspond to how disks
evolve over time. We also consider how the BD disks differ from previously-
studied TTS disks to determine what affect host mass has on certain disk
properties. Disk mass is of particular interest; we assess the implications
the BD disk masses have for planet formation and measure the relationship
between the host masses and disk masses.

###  4.1 Comparison Between Star-Forming Regions

####  4.1.1 Disk Mass

Figure 5: Calculated disk masses for all BDs in this sample, separated by
region. The four regions are ordered by age, increasing left to right. The
blue lines show the total range of disk masses in each region; the width of
the shaded blue regions represent the distributions of the disk mass values.

Average disk masses for TTS have been reported to decrease with age of the
star-forming region (Ansdell et al., 2016; Barenfeld et al., 2016; Pascucci et
al., 2016; Cieza et al., 2019; van Terwisga et al., 2019). Figure 5 shows how
disk mass varies with age for BD disks to determine whether this trend holds
for substellar objects as well. Two main deviations from the trend of
decreasing disk mass with age can be seen: (1) Upper Sco, the oldest region in
the sample, has only slightly lower disk masses than the other regions and (2)
Ophiuchus, the youngest region in the sample, does not have significantly
higher disk masses than the older regions.

The unexpectedly high mass of the Upper Sco disks is likely due to selection
effects and the substellar masses of the host BDs. In order to maximize the
likelihood of detecting disk emission in their ALMA survey of BDs and low mass
stars in Upper Sco, van der Plas et al. (2016) only observed objects with the
largest mid-IR excesses and required the objects to have detections in both
the Herschel PACS 70 $\mu$m and 160 $\mu$m channels. Therefore, the Upper Sco
disks considered here do not represent a complete sample of BD disks in the
region; moreover, the disks that did not meet the mid-IR excess or PACS
detection requirements are likely less massive than the disks included in the
sample.

However, Carpenter et al. (2006), Bayo et al. (2012), and Luhman & Mamajek
(2012) find evidence that disk lifetimes are longer for lower mass hosts.
Though Upper Sco is older than the other regions studied here, we would still
expect to see massive disks around BDs in that region if their disks dissipate
more slowly than their TTS counterparts.

As the youngest region in the sample, disks in Ophiuchus should be the most
massive. However, we find that the disk masses in Ophiuchus are relatively
comparable to the older Taurus and Lupus regions. All three regions have at
least one disk greater than 1 MJ. In fact, Taurus has a slightly greater
fraction of massive disks than Ophiuchus. Williams et al. (2019) also found
TTS disks in Ophiuchus to be less massive than expected. In their study, they
find the Ophiuchus disks to be slightly less massive than TTS disks in the
somewhat older Lupus region. These authors suggest the local cloud environment
could have a substantial effect on the disk masses in a given region. For
example, Kuffmeier et al. (2020) find that ionization can decrease disk size:
disks may be “born” with smaller masses in star-forming regions with higher
ionization rates. Figure 6 plots the locations of the Ophiuchus BDs in this
sample along with the locations of the TTS disks from Williams et al. (2019).
The two populations overlap, so it is reasonable to assume that any local
cloud environment effects that influence the disk masses of the TTS would also
influence the disk masses of the BDs.

Figure 6: Location of the TTS studied by Williams et al. (2019) (black
circles) and the BD disks observed by Testi et al. (2016) and Testi et al.,
submitted and studied in this work (red crosses). Any local cloud environment
effects that influence the disk size of the TTS would affect the BDs as well
since their positions in space overlap.

Of the four star-forming regions studied here, Taurus contains the most
massive disks, including four disks with masses greater than 1 MJup (KPNO Tau
3, J04400067, J04442713, and J04230607), though Lupus and Ophiuchus both
contain at least one equally massive disk (Lup 818s in Lupus, and CFHTWIR-Oph
66 and GY92 264 in Ophiuchus). The high mass of J04442713 was previously
reported by Rilinger et al. (2019); the finding of three additional objects in
the same region with similar masses may indicate that there is a population of
high-mass BD disks in Taurus. As the youngest region at $\sim$0.5 Myr, disks
in Ophiuchus are expected to have not yet dissipated. At $\sim$1 Myr, Taurus
is also a relatively young star-forming region; Lupus is slightly older, at
$\sim$1–3 Myr, but still young enough that significant disk dissipation may
not yet be common.

Figure 7: Comparison of disk dust masses calculated in this work to those
reported by (a) Ward-Duong et al. (2018), (b) Testi et al. (2016), (c) Sanchis
et al. (2020), and (d) van der Plas et al. (2016). In each panel the red line
indicates a 1:1 relationship. Triangles indicate upper limits.

Figure 7 compares the dust disk masses calculated in this work to previously
reported masses in the literature. In particular, we compare our results to
Ward-Duong et al. (2018), Testi et al. (2016), Sanchis et al. (2020), and van
der Plas et al. (2016) for objects in Taurus, Ophiuchus, Lupus, and Upper Sco,
respectively. Ward-Duong et al. (2018) and van der Plas et al. (2016) also use
SED fitting to determine disk masses for their sample of BDs. Both works use
the MCFOST radiative transfer code (Pinte et al., 2006, 2009). As shown in the
panels (a) and (d) of Figure 7, we determine disk masses that are generally in
good agreement, though we find three disks (J04230607, KPNO3, and J04400067)
to be $\sim$2-3 times more massive (note that Ward-Duong et al. (2018) did not
study the fourth massive disk in Taurus, J04442713). None of these three
massive disks have measured $\dot{M}$ values. As discussed in Section 3.1.1,
$\dot{M}$ directly impacts our disk mass calculation; thus, if the $\dot{M}$
values for these objects are significantly lower than the assumed value of
10-10 M⊙ yr-1, these disk masses would be overestimated. Overestimating the
radii of these disks could also result in an overestimate of the disk mass,
since we integrate the surface mass density over the disk radius. High-
resolution ALMA images of these disks would help constrain their radii and
thus the disk mass.

We note that the fourth massive disk reported here, J04442713, was determined
by Rilinger et al. (2019) to have a large disk radius of 100 au based on ALMA
observations, so its high mass is consistent with the disk being large. This
object also has a measured accretion rate of 2 x 10-10 M⊙ yr-1, so its disk
mass is more certain than the other three massive disks.

Testi et al. (2016) and Sanchis et al. (2020) do not use SED fitting to
determine disk mass; rather, they use measurements of the mm flux (Fν) of each
BD to calculate the disk mass according to the following relationship:

$M_{dust}=\frac{F_{\nu}d^{2}}{\kappa_{\nu}B_{\nu}(T)}$ (2)

where d is the distance to the BD, $\kappa_{\nu}$ is the opacity at frequency
$\nu$, and Bν(T) is the Planck function at temperature (T). This relationship
assumes the emission is optically thin and isothermal at T. For most objects
in the Testi et al. (2016) sample, we find disks that are less massive than
their reported values (see Figure 7, panel (b)). This discrepancy could
indicate that either or both of the assumptions of optically thin and/or
isothermal emission are incorrect. The DIAD model does not assume an
isothermal disk, but rather calculates the temperature profile of the disk.
Using radiative transfer models of protoplanetary disks in Taurus, Ballering &
Eisner (2019) find that the optically thin assumption frequently does not hold
and that the disks tend to be somewhat optically thick. This conclusion is
also supported by neural network-based modeling with DIAD performed by Ribas
et al. (2020). Thus, the assumption of optically thin emission in Equation 2
is likely the source of the difference with the Testi et al. (2016) disk mass
estimates. However, the results from Testi et al. (2016) do agree that GY92
264 is a substantially more massive disk than the other disks in Ophiuchus. We
note that Testi et al. (2016) did not study CFHTWIR-Oph 66, the other massive
Ophiuchus disk in our sample.

The masses we determine for the Lupus objects are in somewhat better agreement
with the Sanchis et al. (2020) masses, though our mass for Lup 818s is about 2
times more massive than they report. Unlike J04230607, KPNO3, and J04400067
discussed above, Lup 818s has a measured accretion rate, so uncertainty in
$\dot{M}$ is likely not the cause for the large mass. Furthermore, Sanchis et
al. (2020) determined an upper limit to the disk radius of $\sim$20 au for
this object from their ALMA observations. This compact yet massive disk may be
an interesting object for follow-up ALMA observations.

####  4.1.2 Dust Growth

Figure 8: Distribution of the maximum size of dust grains in the disk
atmosphere, amax,atm, for all BDs in this sample, separated by region. The
four regions are ordered by age, increasing left to right.

As discussed in Section 3.1.2, dust grains are expected to increase in size as
disks evolve. Figure 8 shows the distribution of the best-fit atmospheric dust
grain maximum sizes from our models for the BD disks in each of the four star-
forming regions studied here. The histograms show that all regions, even the
youngest, have undergone some amount of grain growth. All disks in Upper Sco,
the oldest region, have grains larger than typical ISM grains; in fact all but
one of the disks in this region have grains about two orders of magnitude
larger than ISM-sized grains.

The result that the oldest star-forming region in our sample has the largest
dust grains is consistent with the idea that the oldest disks should be the
most evolved. It is important to note, however, that all four regions contain
at least some disks with evidence of grain growth. Grain growth occuring even
in young star-forming regions is consistent with previous studies of TTS disks
in Taurus, Ophiuchus, Lupus, and Chamaeleon, which report grain growth to
millimeter sizes in the disk midplanes based on spectral indices (e.g., Lommen
et al., 2007; Ricci et al., 2010b, a; Ribas et al., 2017). Grain growth at
young ages supports the hypothesis that planets may form quickly (e.g., Manara
et al., 2018, see Section 4.2 of this work.).

####  4.1.3 Dust Settling

Figure 9: Distribution of the dust settling parameter $\epsilon$ for all BDs
in this sample, separated by region. The four regions are ordered by age,
increasing left to right.

Figure 9 shows the distribution of the best-fit settling parameter $\epsilon$
from our models for the BD disks in each of the four star-forming regions
studied here. Lower values of $\epsilon$ correspond to higher degrees of
settling. We expect the degree of settling to increase with age and amount of
grain growth. Indeed, Ophiuchus has the largest fraction of high $\epsilon$
values (i.e., the least settled disks) and Upper Sco shows the lowest average
$\epsilon$ value in its disks, which is consistent with its larger average
grain size and older age. However, the three younger regions contain some
settled disks as well, though they also have less-settled disks than Upper
Sco. Evidence of settling in young disks has previously been reported by
Furlan et al. (2009), Lewis & Lada (2016), and Grant et al. (2018), among
others. Along with evidence of grain growth (see previous section) and disk
substructures (e.g., HL Tau ALMA Partnership et al., 2015) in young ($\sim$1
Myr old or younger) objects, these low $\epsilon$ values are consistent with
disk evolution occurring sooner than previously thought.

Figure 10 compares the settling measured for Taurus BDs in this sample to the
settling of TTS disks in Taurus as modeled by Grant et al. (2018). Grant et
al. (2018) used a forward modeling process and DIAD to predict the
distribution of $\epsilon$ for TTS disks in the Taurus Molecular Cloud. This
distribution is not a histogram of the $\epsilon$ values from best-fit SEDs
for individual objects. However, it is still useful to compare this
distribution to our histogram of best-fit $\epsilon$ values for BDs in Taurus
in order to search for similarities or differences in the amount of settling
for these two types of objects.

For both BDs and TTS in Taurus, the distribution of $\epsilon$ values peaks at
lower epsilon values. Thus, both populations of objects appear to display
similarly high degrees of settling. BDs appear to evolve similarly to TTS over
time, with perhaps longer lifetimes as reported by Luhman & Mamajek (2012).

Figure 10: Distribution of the dust settling parameter $\epsilon$ for BDs
(left) and TTS (right) in Taurus. Right panel adapted from Grant et al.
(2018). Error bars are standard deviations of the forward-modeling
realizations performed by Grant et al. (2018).

###  4.2 Disk Mass and Planet Formation

The amount of mass in a disk is directly linked to the potential for planet
formation. In their simulations of planet formation in the disks around BDs,
Payne & Lodato (2007) found that Earth-mass planets are able to form easily in
disks of $\sim$1 Jupiter mass or more, but lower mass disks were very unlikely
to form planets of any mass. These authors ascribe the lack of planets in
lower mass disks to the lower disk surface density. Taking the 1 MJ lower
limit as a “cutoff” for planet formation, very few disks in our sample are
likely to form Earth-mass planets and the possibility of planets forming in
the bulk of BD disks is therefore slim. Only seven disks, discussed in Section
4.1, have masses above this cutoff.

As shown in Figure 7, previous works also report small disk masses for BD
disks. Furthermore, studies of protoplanetary disks around TTS, e.g., Najita &
Kenyon (2014), find masses too small to account for observed planetary
systems. However, though the current masses of these disks are too small for
planet formation, our models do not rule out the possibility that planets may
have already formed in the disks. Manara et al. (2018) suggested rapid planet
formation within $<$0.1–1 Myr could explain the discrepancy between disk mass
and observed planetary systems. Observations of disk substructures in 0.1–1
Myr protostellar disks (Sheehan et al., 2020) as well as grain growth and
settling in $\sim$1 Myr-old disks (Grant et al., 2018) are evidence of early-
onset disk evolution and possibly planet formation.

Alternatively, if grain growth in BD disks is efficient (i.e., grains stick
together upon collision without splintering, and photoevaporation is
inefficient), a significant amount of the disk material could be composed of
bodies larger than $\sim$1 m but smaller than planets. These bodies would not
contribute to the emissivity at millimeter wavelengths, and thus would not
have an effect on the SEDs of the disks. Yet another possibility is that disks
are replenished over their lifetimes by material from the surrounding
interstellar medium (Manara et al., 2018), so that the total amount of planet-
forming material available is greater than the amount of material in the disk
at any given time. However, the exact mechanism by which planets form in these
disks is still unknown, and the origin of observed planetary-mass companions
to BDs remains mysterious.

###  4.3 Disk Mass – Host Mass Relationship

Sub-millimeter surveys of disks in various star-forming regions have
established a positive correlation between the masses of the disks and the
masses of the stars that host them (e.g., Andrews et al., 2013; Ansdell et
al., 2016; Barenfeld et al., 2016; Pascucci et al., 2016). For the most part,
these surveys have focused on stellar-mass hosts and have not included many
substellar-mass (i.e., BD) hosts. In this section, we explore how this mass
relationship changes when BD disks are included in the sample.

Following Pascucci et al. (2016), we use the Bayesian linear fitting procedure
developed by Kelly (2007). This procedure accounts for the measurement
uncertainties, upper limit values, and intrinsic scatter in both the disk mass
and host mass values to obtain a linear relationship between disk mass and
host mass in the log-log plane. We adopt the disk mass and host mass values
for objects in Taurus, Lupus, and Upper Scorpius reported by Pascucci et al.
(2016), who compiled data from Andrews et al. (2013), Ansdell et al. (2016)
and Barenfeld et al. (2016) for these three regions, respectively. Pascucci et
al. (2016) did not study the disk mass–host mass relationship for Ophiuchus,
so we did not explore how BDs affect the relationship in that region. For
objects in both the Pascucci et al. (2016) sample and our sample of BDs, we
adopt the disk mass obtained by our modeling. The linear fits to this sample,
with and without our sample of BDs included, are shown in Figure 11 and
summarized in Table 6.

As Figure 11 shows, including BD disks does not substantially change the disk
mass – host mass relationship. The two relationships (with and without our BD
sample) are almost exactly the same in Lupus. The BDs appear to slightly
flatten the relationship in Taurus and Upper Sco, but the linear fits are
still within the uncertainties of the relationships without the BD sample, so
the effect is not significant. Furthermore, as discussed in Section 4.1.1, the
disks in Upper Sco modeled here may constitute a biased sample; a more
complete sample would likely include objects with lower disk masses that would
bring the disk mass – host mass relationship into closer agreement with the
Pascucci et al. (2016) result. A continuous Mdisk – M∗ relationship between
stellar-mass hosts and BD-mass hosts supports the theory of a single shared
formation mechanism for stars and BDs (e.g., Hennebelle & Chabrier, 2009;
Ricci et al., 2013, 2014; Rilinger et al., 2019).

Table 6: Mdust – M∗ relationships  Region | $\alpha$ | $\beta$
---|---|---
Without our BD sample
Taurus | 1.31 $\pm$ 0.24 | 1.13 $\pm$ 0.11
Lupus | 1.10 $\pm$ 0.30 | 1.34 $\pm$ 0.18
Upper Sco | 1.73 $\pm$ 0.39 | 0.77 $\pm$ 0.23
With our BD sample
Taurus | 0.99 $\pm$ 0.14 | 0.98 $\pm$ 0.10
Lupus | 1.13 $\pm$ 0.24 | 1.35 $\pm$ 0.16
Upper Sco | 1.10 $\pm$ 0.31 | 0.53 $\pm$ 0.20
Figure 11: Mdisk versus M∗ relationships for disks in (a) Taurus, (b) Lupus,
and (c) Upper Sco. Open circles represent the sample reported in Pascucci et
al. (2016), filled circles represent the BDs studied in this work. The black
line and grey shaded region represent the best linear fit and uncertainties
for the Pascucci et al. (2016) sample; the red line and red shaded region
represent the best linear fit and uncertainties for the Pascucci et al. (2016)
sample plus the BDs studied in this work. Note that some BDs have very similar
host and disk masses, so their points overlap completely.

##  5 Summary and Conclusions

We present here the largest sample of SED models of BD protoplanetary disks to
date. We consider 49 objects located in four different star-forming regions in
the Milky Way. The DIAD radiative transfer code is used to fit models to the
SEDs of these objects. From the best-fit SED models, we obtain dust grain
sizes, the degree of dust settling, and the disk mass, among other parameters,
for each disk. We compare the four star-forming regions to each other to
assess how disk mass, grain growth, and dust settling vary with age. Our
results are summarized as follows:

  1. 1.

Though Ophiuchus is the youngest region, and thus would be expected to have
the most massive disks, its disk masses are comparable to the older Taurus and
Lupus regions. This result is consistent with the TTS disk masses reported by
Williams et al. (2019).

  2. 2.

Upper Scorpius is the oldest region, but still contains some relatively
massive disks. This may be a selection effect, since these disks were chosen
due to their high mid-IR fluxes. However, Luhman & Mamajek (2012) note that
disk lifetimes may be longer for substellar hosts, and the disk may not
dissipate as quickly.

  3. 3.

All four regions show some evidence of grain growth and dust settling,
indicating that dust evolution may begin to occur within 1 Myr of disk
formation.

  4. 4.

Very few of the disks studied here are massive enough to form Earth-mass
planetary companions according to the Payne & Lodato (2007) simulations. Rapid
planet formation and/or disk replenishment from the interstellar medium may be
necessary to explain observed planetary companions to BDs. We do not rule out
the possibility that planets may form earlier and may already be present in
the disks.

  5. 5.

The disk masses obtained in our modeling are consistent with previously
published disk mass – host mass relationships for stellar-mass hosts, which
lends support to the theory of a shared formation mechanism between stars and
BDs.

Unbiased millimeter-wavelength studies of BD disks in these and additional
star-forming regions are needed to obtain a more complete picture of how BD
disks evolve over time. The launch of the James Webb Space Telescope will also
enable detailed studies of the dust and gas in these disks, which will improve
upon the disk properties presented here.

We thank the referee for their careful reading of the manuscript and their
helpful feedback. This paper utilizes the D’Alessio Irradiated Accretion Disk
(DIAD) code. We wish to recognize the work of Paola D’Alessio, who passed away
in 2013. Her legacy and pioneering work live on through her substantial
contributions to the field. This work was funded by NASA ADAP 80NSSC20K0451.

##  Appendix A SED Models

Figure A.1: Spectral energy distribution models for the 23 objects in Taurus.
Filled blue triangles represent upper limits (continued on next page). (The
observed photometry shown in this figure are available electronically.) Figure
A.1: Continued: Spectral energy distribution models for the 23 objects in
Taurus. Filled blue triangles represent upper limits. (The observed photometry
shown in this figure are available electronically.) Figure A.3: Spectral
energy distribution models for the 13 objects in Ophiuchus. Filled blue
triangles represent upper limits. (The observed photometry shown in this
figure are available electronically.) Figure A.4: Spectral energy distribution
models for the 8 objects in Lupus. Filled blue triangles represent upper
limits. (The observed photometry shown in this figure are available
electronically.) Figure A.5: Spectral energy distribution models for the 5
objects in Upper Scorpius. Filled blue triangles represent upper limits. (The
observed photometry shown in this figure are available electronically.)

## References

  * Adame et al. (2011) Adame, L., Calvet, N., Luhman, K. L., et al. 2011, ApJ, 726, L3
  * Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  * André et al. (2012) André, P., Ward-Thompson, D., & Greaves, J. 2012, Science, 337, 69
  * Andrews et al. (2013) Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129
  * Ansdell et al. (2017) Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46
  * Apai et al. (2005) Apai, D., Pascucci, I., Bouwman, J., et al. 2005, Science, 310, 834
  * Ballering & Eisner (2019) Ballering, N. P., & Eisner, J. A. 2019, AJ, 157, 144
  * Barenfeld et al. (2016) Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142
  * Basu & Vorobyov (2012) Basu, S., & Vorobyov, E. I. 2012, ApJ, 750, 30
  * Bate et al. (2003) Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577
  * Bayo et al. (2012) Bayo, A., Barrado, D., Huélamo, N., et al. 2012, A&A, 547, A80
  * Bonfils et al. (2013) Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109
  * Bouy et al. (2008) Bouy, H., Huélamo, N., Pinte, C., et al. 2008, A&A, 486, 877
  * Calvet et al. (1991) Calvet, N., Patino, A., Magris, G. C., & D’Alessio, P. 1991, ApJ, 380, 617
  * Carpenter et al. (2006) Carpenter, J. M., Mamajek, E. E., Hillenbrand, L. A., & Meyer, M. R. 2006, ApJ, 651, L49
  * Cazzoletti et al. (2019) Cazzoletti, P., Manara, C. F., Baobab Liu, H., et al. 2019, A&A, 626, A11
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698
  * Cohen & Witteborn (1985) Cohen, M., & Witteborn, F. C. 1985, ApJ, 294, 345
  * Comerón (2008) Comerón, F. 2008, in Handbook of Star Forming Regions, Volume II, ed. B. Reipurth, Vol. 5 (ASP), 295
  * Comeron et al. (1998) Comeron, F., Rieke, G. H., Claes, P., Torra, J., & Laureijs, R. J. 1998, A&A, 335, 522
  * Daemgen et al. (2016) Daemgen, S., Natta, A., Scholz, A., et al. 2016, A&A, 594, A83
  * D’Alessio et al. (2001) D’Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321
  * D’Alessio et al. (2006) D’Alessio, P., Calvet, N., Hartmann, L., Franco-Hernández, R., & Servín, H. 2006, ApJ, 638, 314
  * D’Alessio et al. (1999) D’Alessio, P., Calvet, N., Hartmann, L., Lizano, S., & Cantó, J. 1999, ApJ, 527, 893
  * D’Alessio et al. (1998) D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411
  * D’Alessio et al. (2005) D’Alessio, P., Hartmann, L., Calvet, N., et al. 2005, ApJ, 621, 461
  * Dorschner et al. (1995) Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503
  * Draine & Lee (1984) Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89
  * Dressing & Charbonneau (2013) Dressing, C. D., & Charbonneau, D. 2013, ApJ, 767, 95
  * Dullemond & Dominik (2004) Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075
  * Eisner et al. (2016) Eisner, J. A., Bally, J. M., Ginsburg, A., & Sheehan, P. D. 2016, ApJ, 826, 16
  * Eisner et al. (2018) Eisner, J. A., Arce, H. G., Ballering, N. P., et al. 2018, ApJ, 860, 77
  * Ercolano et al. (2009) Ercolano, B., Clarke, C. J., & Robitaille, T. P. 2009, MNRAS, 394, L141
  * Espaillat et al. (2007) Espaillat, C., Calvet, N., D’Alessio, P., et al. 2007, ApJ, 670, L135
  * Espaillat et al. (2011) Espaillat, C., Furlan, E., D’Alessio, P., et al. 2011, ApJ, 728, 49
  * Espaillat et al. (2010) Espaillat, C., D’Alessio, P., Hernández, J., et al. 2010, ApJ, 717, 441
  * Espaillat et al. (2012) Espaillat, C., Ingleby, L., Hernández, J., et al. 2012, ApJ, 747, 103
  * Espaillat et al. (2014) Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 497
  * Furlan et al. (2009) Furlan, E., Watson, D. M., McClure, M. K., et al. 2009, ApJ, 703, 1964
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Grant et al. (2018) Grant, S. L., Espaillat, C. C., Megeath, S. T., et al. 2018, ApJ, 863, 13
  * Grant et al. (2021) Grant, S. L., Espaillat, C. C., Wendeborn, J., et al. 2021, arXiv e-prints, arXiv:2103.16453
  * Gürtler et al. (1999) Gürtler, J., Schreyer, K., Henning, T., Lemke, D., & Pfau, W. 1999, A&A, 346, 205
  * Han et al. (2013) Han, C., Jung, Y. K., Udalski, A., et al. 2013, ApJ, 778, 38
  * Harvey et al. (2012) Harvey, P. M., Henning, T., Liu, Y., et al. 2012, ApJ, 755, 67
  * Hendler et al. (2017) Hendler, N. P., Mulders, G. D., Pascucci, I., et al. 2017, ApJ, 841, 116
  * Hennebelle & Chabrier (2009) Hennebelle, P., & Chabrier, G. 2009, ApJ, 702, 1428
  * Herczeg & Hillenbrand (2008) Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594
  * Houck et al. (2004) Houck, J. R., Roellig, T. L., Van Cleve, J., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5487, Proc. SPIE, ed. J. C. Mather, 62–76
  * Kataoka et al. (2014) Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42
  * Kataoka et al. (2013) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4
  * Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
  * Kenyon & Hartmann (1995) Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117
  * Kim et al. (2013) Kim, K. H., Watson, D. M., Manoj, P., et al. 2013, ApJ, 769, 149
  * Kuffmeier et al. (2020) Kuffmeier, M., Zhao, B., & Caselli, P. 2020, A&A, 639, A86
  * Lewis & Lada (2016) Lewis, J. A., & Lada, C. J. 2016, ApJ, 825, 91
  * Liu et al. (2015) Liu, Y., Joergens, V., Bayo, A., Nielbock, M., & Wang, H. 2015, A&A, 582, A22
  * Lommen et al. (2007) Lommen, D., Wright, C. M., Maddison, S. T., et al. 2007, A&A, 462, 211
  * Luhman (2004) Luhman, K. L. 2004, ApJ, 617, 1216
  * Luhman et al. (2010) Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111
  * Luhman & Mamajek (2012) Luhman, K. L., & Mamajek, E. E. 2012, ApJ, 758, 31
  * Luhman et al. (2005) Luhman, K. L., Lada, C. J., Hartmann, L., et al. 2005, ApJ, 631, L69
  * Luhman et al. (2007) Luhman, K. L., Adame, L., D’Alessio, P., et al. 2007, ApJ, 666, 1219
  * Manara et al. (2018) Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3
  * Manara et al. (2015) Manara, C. F., Testi, L., Natta, A., & Alcalá, J. M. 2015, A&A, 579, A66
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * Meru et al. (2013) Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 774, L4
  * Monnier & Millan-Gabet (2002) Monnier, J. D., & Millan-Gabet, R. 2002, ApJ, 579, 694
  * Morrow et al. (2008) Morrow, A. L., Luhman, K. L., Espaillat, C., et al. 2008, ApJ, 676, L143
  * Muench et al. (2001) Muench, A. A., Alves, J., Lada, C. J., & Lada, E. A. 2001, ApJ, 558, L51
  * Muzerolle et al. (2010) Muzerolle, J., Allen, L. E., Megeath, S. T., Hernández, J., & Gutermuth, R. A. 2010, ApJ, 708, 1107
  * Muzerolle et al. (2006) Muzerolle, J., Adame, L., D’Alessio, P., et al. 2006, ApJ, 643, 1003
  * Najita & Kenyon (2014) Najita, J. R., & Kenyon, S. J. 2014, MNRAS, 445, 3315
  * Natta et al. (2000) Natta, A., Meyer, M. R., & Beckwith, S. V. W. 2000, ApJ, 534, 838
  * Natta & Testi (2001) Natta, A., & Testi, L. 2001, A&A, 376, L22
  * Natta et al. (2002) Natta, A., Testi, L., Comerón, F., et al. 2002, A&A, 393, 597
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23
  * Ortiz-León et al. (2017) Ortiz-León, G. N., Loinard, L., Kounkel, M. A., et al. 2017, ApJ, 834, 141
  * Owen (2016) Owen, J. E. 2016, PASA, 33, e005
  * Pascucci et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125
  * Payne & Lodato (2007) Payne, M. J., & Lodato, G. 2007, MNRAS, 381, 1597
  * Pecaut & Mamajek (2013) Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9
  * Pecaut et al. (2012) Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154
  * Pinilla et al. (2013) Pinilla, P., Birnstiel, T., Benisty, M., et al. 2013, A&A, 554, A95
  * Pinilla et al. (2012) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114
  * Pinilla et al. (2017) Pinilla, P., Pohl, A., Stammler, S. M., & Birnstiel, T. 2017, ApJ, 845, 68
  * Pinilla et al. (2014) Pinilla, P., Benisty, M., Birnstiel, T., et al. 2014, A&A, 564, A51
  * Pinte et al. (2009) Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967
  * Pinte et al. (2006) Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797
  * Preibisch et al. (2002) Preibisch, T., Brown, A. G. A., Bridges, T., Guenther, E., & Zinnecker, H. 2002, AJ, 124, 404
  * Reipurth & Clarke (2001) Reipurth, B., & Clarke, C. 2001, AJ, 122, 432
  * Riaz et al. (2018) Riaz, R., Vanaverbeke, S., & Schleicher, D. R. G. 2018, MNRAS, 478, 5460
  * Ribas et al. (2020) Ribas, Á., Espaillat, C. C., Macías, E., & Sarro, L. M. 2020, A&A, 642, A171
  * Ribas et al. (2017) Ribas, Á., Espaillat, C. C., Macías, E., et al. 2017, ApJ, 849, 63
  * Ricci et al. (2013) Ricci, L., Isella, A., Carpenter, J. M., & Testi, L. 2013, ApJ, 764, L27
  * Ricci et al. (2010a) Ricci, L., Testi, L., Natta, A., & Brooks, K. J. 2010a, A&A, 521, A66
  * Ricci et al. (2010b) Ricci, L., Testi, L., Natta, A., et al. 2010b, A&A, 512, A15
  * Ricci et al. (2014) —. 2014, ApJ, 791, 20
  * Rilinger et al. (2019) Rilinger, A. M., Espaillat, C. C., & Macías, E. 2019, ApJ, 878, 103
  * Ruíz-Rodríguez et al. (2018) Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2018, MNRAS, 478, 3674
  * Sanchis et al. (2020) Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114
  * Sargent et al. (2009) Sargent, B. A., Forrest, W. J., Tayrien, C., et al. 2009, ApJS, 182, 477
  * Scholz et al. (2007) Scholz, A., Jayawardhana, R., Wood, K., et al. 2007, ApJ, 660, 1517
  * Sheehan et al. (2020) Sheehan, P. D., Tobin, J. J., Federman, S., Megeath, S. T., & Looney, L. W. 2020, ApJ, 902, 141
  * Shvartzvald et al. (2017) Shvartzvald, Y., Yee, J. C., Calchi Novati, S., et al. 2017, ApJ, 840, L3
  * Stamatellos & Whitworth (2009) Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 392, 413
  * Sterzik et al. (2004) Sterzik, M. F., Pascucci, I., Apai, D., van der Bliek, N., & Dullemond, C. P. 2004, A&A, 427, 245
  * Testi et al. (2016) Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111
  * Udalski et al. (2015) Udalski, A., Jung, Y. K., Han, C., et al. 2015, ApJ, 812, 47
  * van der Plas et al. (2016) van der Plas, G., Ménard, F., Ward-Duong, K., et al. 2016, ApJ, 819, 102
  * van Terwisga et al. (2019) van Terwisga, S. E., Hacar, A., & van Dishoeck, E. F. 2019, A&A, 628, A85
  * van Terwisga et al. (2020) van Terwisga, S. E., van Dishoeck, E. F., Mann, R. K., et al. 2020, A&A, 640, A27
  * Ward-Duong et al. (2018) Ward-Duong, K., Patience, J., Bulger, J., et al. 2018, AJ, 155, 54
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Wilking et al. (2008) Wilking, B. A., Gagné, M., & Allen, L. E. 2008, in Handbook of Star Forming Regions, Volume II, ed. B. Reipurth, Vol. 5 (ASP), 351
  * Williams et al. (2019) Williams, J. P., Cieza, L., Hales, A., et al. 2019, ApJ, 875, L9
  * Williams & Cieza (2011) Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67


Implications of PSR J1614-2230 for NJL hybrid star

César Henrique Lenzi
Universidade Tecnológica Federal do Paraná - Campus Medianeira
85884-000 - Vila São Cristovão
Medianeira, PR
Brazil
Email<EMAIL_ADDRESS>

Germán Lugones
Universidade Federal do ABC - Campus Santo André
09210-580- Bangu
Santo André, SP
Brazil
Email<EMAIL_ADDRESS>

##  1 Introduction

The recent determination of the mass of the pulsar PSR J1614-2230 with
$1.97\pm 0.04M_{\odot}$ by [1], renewed the discussions about the possibility
of exotic matter being present at the core of neutron stars. Since the
description of matter at densities beyond nuclear saturation is model
dependent, several works have explored different aspects of the fact that the
maximum neutron star mass implied by any equation of state (EoS) must exceed
the mass of PSR J1614-2230 [2, 3]

In this article we present an extensive study of hybrid star masses using
several parametrizations of a relativistic mean-field hadronic EoS together
with a typical three-flavor NJL model with scalar, vector and ’t Hooft
interactions as realized in [2].

##  2 The model

To describe the quark matter phase we use the SU(3) NJL model with scalar-
pseudoscalar, isoscalar-vector and ’t Hooft six fermion interaction. The
Lagrangian density of the model is:

$\displaystyle{\cal L}_{Q}$ $\displaystyle=$
$\displaystyle\overline{\psi}(i\gamma_{\mu}\partial^{\mu}-\hat{m})\psi$ (1)
$\displaystyle+$ $\displaystyle
g_{s}\sum_{a=0}^{8}[(\overline{\psi}\lambda^{a}\psi)^{2}+(\overline{\psi}i\gamma_{5}\lambda^{a}\psi)^{2}]$
$\displaystyle-$ $\displaystyle
g_{v}\sum_{a=0}^{8}[(\overline{\psi}\gamma_{\mu}\lambda^{a}\psi)^{2}+(\overline{\psi}\gamma_{5}\gamma_{\mu}\lambda^{a}\
\psi)^{2}]$ $\displaystyle+$ $\displaystyle
g_{t}\\{\det[\overline{\psi}(1+\gamma_{5})\psi]+\det[\overline{\psi}(1-\gamma_{5})\psi]\\},$

where $\psi=(u,d,s)$ denotes the quark fields, $\lambda^{a}(0\leq a\leq 8)$
are the U(3) flavour matrices, $\hat{m}=\mbox{diag}({m}_{u},{m}_{d},{m}_{s})$
is the quark current mass, and $g_{s}$, $g_{v}$ and $g_{t}$ are coupling
constants. The mean-field thermodynamic potential density $\Omega$ for a given
baryon chemical potential $\mu$ at $T=0$, is given by

$\displaystyle\Omega=$ $\displaystyle-$ $\displaystyle\eta
N_{c}\sum_{i}\int_{k_{Fi}}^{\Lambda}{\frac{p^{2}\,dp}{2\pi^{2}}}\sqrt{p^{2}+M_{i}^{2}}+2g_{s}\sum_{i}\langle\overline{\psi}\psi\rangle_{i}^{2}$
(2) $\displaystyle-$ $\displaystyle
2g_{v}\sum_{i}\langle\psi^{\dagger}\psi\rangle_{i}^{2}+4g_{t}\langle\overline{u}u\rangle\langle\overline{d}d\rangle\langle\overline{s}s\rangle$
$\displaystyle-$ $\displaystyle\eta
N_{c}\sum_{i}\mu_{i}{\int_{0}^{k_{Fi}}\frac{p^{2}\,dp}{2\pi^{2}}}-\Omega_{0},$

where the sum is over the quark flavor $(i=u,d,s)$, the constants $\eta=2$ and
$N_{c}=3$ are the spin and color degeneracies, and $\Lambda$ is a
regularization ultraviolet cutoff to avoid divergences in the medium
integrals. The Fermi moment of the particle $i$ is given by
$k_{Fi}=\theta(\mu^{\ast}_{i}-M_{i})\sqrt{(}\mu^{\ast 2}_{i}-M_{i}^{2})$,
where $\mu^{\ast}_{i}$ is the quark chemical potential modified by the
vectorial interaction, i.e.
$\mu^{\ast}_{u,d,s}=\mu_{u,d,s}-4g_{v}\langle\psi^{\dagger}\psi\rangle_{u,d,s}$.

The conventional procedure for fixing the $\Omega_{0}$ term in Eq. (2) is to
assume that the grand thermodynamic potential $\Omega$ must vanish at zero
$\mu$ and $T$. Nevertheless, this prescription is no more than an arbitrary
way to uniquely determine the EoS of the NJL model without any further
assumptions [4]. In view of this, [5] adopt a different strategy. They fix a
bag constant for the hadron-quark deconfinement to occur at the same chemical
potential as the chiral phase transition. This method leads to a significant
change in the EoS with respect to the conventional procedure. Differently, in
[2] we explore the above possibility of having chiral restoration and
deconfinement occurring at different densities. To this end, we shall
substitute $\Omega_{0}$ in Eq. (2) by the new value
$\Omega_{0}+\delta\Omega_{0}$, where $\delta\Omega_{0}$ is a free parameter:

$\displaystyle\Omega_{0}\longrightarrow\Omega_{0}+\delta\Omega_{0}\ \ \ \ \
\mathrm{in}\ \mathrm{Eq.}\ (\ref{eq3}).$ (3)

With this change, the thermodynamic potential $\Omega$ can be non-vanishing at
zero $\mu$ and $T$, and the $\mu$ of the deconfinement transition can be
tuned. For more details, see [2].

##  3 Results

We have solved the Tolman-Oppenheimer-Volkoff equations for spherically
symmetric and static stars in order to investigate the influence of $g_{v}$
and $\delta\Omega_{0}$ on the maximum mass of hybrid stars. In Figure 1 we
have represented the maximum mass of hybrid stars for different
parametrizations of the NJL model (more deals see the label of Figure). An
interesting feature of Figure 1, is that large masses are situated on the
right-upper corner but stable configurations are located on the left-lower
corner of the figure (or left side of the figure in the case of NL3). This
clearly illustrates the difficulty of obtaining stable hybrid stars with
arbitrarily large masses. Concerning the effect of the hadronic model we see
that stable hybrid stars have higher values of the maximum mass for the
stiffer hadronic EoS.

Figure 1: Background colors represent the maximum mass of hybrid stars for
different parametrizations of the NJL model (i.e. different values of $g_{v}$
and $\delta\Omega_{0}$). In each panel we use a different hadronic EoS
(without hyperons): (a) GM1, (b) TM1 and (c) NL3. Notice that the color scale
is different for each panel. The solid contour lines indicate specific values
of the maximum mass. The black solid line represents the boundary between
parametrizations that allow for stable hybrid stars and parametrizations that
do not. The red dashed line indicates the value $1.97{M}_{\odot}$
corresponding to the observed mass of PSR J1614-2230 [1]. The region between
the red dashed line and the solid black line allows to explain the mass of PSR
J1614-2230. Font: Reference [2].

The observed mass of PSR J1614-2230 can be explained by parameters within the
large region located between the red dashed line and the solid black line in
each panel of Figure 1. However, a hypothetical future observation of a
neutron star with a mass a $\sim 10\%$ larger than the mass of PSR J1614-2230
will be hard to explain within hybrid star models using the GM1 and TM1 EoS
and will require a very stiff hadronic model such as NL3.

##  4 Conclusions

In this work we performed a systematic study of hybrid star configurations
using a relativistic mean-field hadronic EoS and the NJL model for three-
flavor quark matter. For the hadronic phase we used the stiff GM1 and TM1
parametrizations, as well as the very stiff NL3 model. In the NJL Lagrangian
we included scalar, vector and ’t Hooft interactions. The vector coupling
constant $g_{v}$ was treated as a free parameter. We also considered that
there is an arbitrary split between the deconfinement and the chiral phase
transitions. This split can be adjusted by a redefinition of the constant
parameter $\Omega_{0}$ in the NJL thermodynamic potential. The effect of this
redefinition is clear in Figure 1, where we show the maximum mass of static
spherically symmetric stars in the parameter space of $g_{v}$ and
$\delta\Omega_{0}$. The larger masses are situated on the right-upper corner
of the diagram, where both $g_{v}$ and $\delta\Omega_{0}$ are larger. On the
other hand, stable configurations are placed on the opposite part of the
diagram; i.e. on the left-lower corner for the GM1 and TM1 models and on the
left side for the NL3 model. As a consequence, stable configurations with a
maximum mass compatible with PSR J1614-2230 are located halfway, specifically,
between the red dashed line and the solid black line of Figure 1. The effect
of the hadronic model (with nucleons only) is also clear from Figure 1, where
we see that stable hybrid stars have higher values of the maximum mass for the
stiffer hadronic EoS. More details about this work can be seen in [2].

## References

  * [1] P.B. Demorest, T. Pennucci, S.M. Ransom, M.S.E. Roberts, J.W.T. Hessels, Nature, 467, 1081 (2010).
  * [2] C.H. Lenzi, G. Lugones, ApJ, 759, 57 (2012).
  * [3] L. Bonanno, A. Sedrakian, A&A, 539, A16 (2012).
  * [4] K. Schertler, S. Leupold, J. Schaffner-Bielich, Phys. Rev. C, 60, 025801 (1999).
  * [5] G. Pagliara, J. Schaffner-Bielich, Phys. Rev. D, 77, 063004 (2008).