A millimeter Continuum Size–Luminosity Relationship for Protoplanetary Disks

A sample of 50 nearby ($d\leqslant 200$ pc) disk targets was collated from the archived catalog of ∼340 GHz (880 μm) continuum measurements made with the Submillimeter Array (SMA; Ho et al. 2004), since the start of science operations in 2004. We considered four primary selection criteria to include a given target in this sample: (1) a quality calibration (i.e., suitable phase stability); (2) a 340 GHz continuum flux density ≥20 mJy; (3) baseline lengths ≥200 m (to ensure sufficient resolution for estimating an emission size); and (4) no known companion within a ∼2'' separation.³ The first three of these are practical restrictions to ensure reliable measurements of the metrics of interest. The fourth criterion has a physical motivation: it was designed to avoid targets where dynamical interactions in multiple star systems are known to reduce the continuum sizes and luminosities (e.g., Jensen et al. 1996; Harris et al. 2012).

Of the 50 disks in our survey, 10 were recently observed by us expressly for the purposes of the present study. To our knowledge, the SMA observations of 18 targets have not yet been published elsewhere. The data for the remaining 32 targets have appeared previously in the literature (Section 2.3). Table 1 is a brief SMA observation log, with references for where the data originally appeared. Targets in this sample are primarily located in the Taurus and Ophiuchus star-forming regions, although nine are in other regions or are found in isolation.

We observed 10 disks in the Taurus–Auriga complex (DN Tau, FY Tau, GO Tau, HP Tau, IP Tau, IQ Tau, IRAS 04385+2550, 2MASS J04333905+2227207, 2MASS J04334465+2615005, and MHO 6) using the eight 6 m SMA antennas arranged in their compact, extended, and very extended configurations (9–508 m baseline lengths). The dual-sideband 345 GHz receivers were tuned to a local oscillator (LO) frequency of 341.6 GHz (878 μm). The SMA correlator processed two intermediate frequency (IF) bands, at ±4–6 and ±6–8 GHz from the LO. Each IF band contains 24 spectral chunks of 104 MHz width. One chunk in the lower IF was split into 256 channels; all others were divided into 32 coarser channels. Each observation cycled between multiple targets and the nearest bright quasars, 3C 111 and J0510+180. Additional measurements of 3C 454.3, Uranus, and Callisto were made for bandpass and flux calibration, respectively.

Another eight disks (CI Tau, CY Tau, DO Tau, Haro 6–13, SU Aur, V410 X-ray 2, V836 Tau, and IM Lup) had previously unpublished data that are suitable for this program available in the SMA archive. These observations have slight variations on the correlator setup and calibrator sources noted above, but in the end produce effectively similar continuum data products.

The raw visibility data were calibrated using standard procedures with the facility MIR package. After calibrating the spectral response of the system and setting the amplitude scale, gain variations were corrected based on repeat observations of the nearest quasar. For observations taken over long time baselines, the visibility phases were shifted to align the data: usually this is based on the known proper motion of the stellar host, but in a few cases where this was unavailable we relied on Gaussian fits to individual measurements. All visibility data for a given target were combined, after confirming their consistency on overlapping baselines. Each visibility set was Fourier-inverted, deconvolved with the clean algorithm, and then convolved with a restoring beam to synthesize an image. The imaging process was conducted with the MIRIAD software package. Table 2 lists basic imaging parameters for the composite data sets. Figure 1 shows a gallery of images for all targets in the sample.

Zoom In Zoom Out Reset image size

The bulk of the sample, 32 targets, have SMA continuum data that were already published in the literature. These disks are primarily in the Ophiuchus star-forming region, but also include members of the Taurus and Lupus complexes, as well as assorted isolated sources. The details of these observations and their calibration are provided in the original references (see Table 1), but they are generally similar to those presented in Section 2.2.

High-quality observations of disk continuum emission have previously been modeled with "broken power-law" (e.g., Andrews et al. 2012; de Gregorio-Monsalvo et al. 2013; Hogerheijde et al. 2016) or "similarity solution" surface brightness profiles (e.g., Andrews et al. 2009, 2010; Isella et al. 2009, 2010). For the survey data presented here, experimentation with these brightness profiles demonstrated that some targets were much better described by one of these options, while the other left significant, persistent residuals. Instead of using different models for different targets, we elected to adopt a more flexible prescription to interpret the data for the full sample in a homogeneous context.

The required flexibility in this prescription is a mechanism that treats either smooth or sharp transitions in the brightness profile. Fortunately, analogous studies of brightness profiles of elliptical galaxies provide some useful guidance in this context. Lauer et al. (1995) introduced a prescription (the so-called "Nuker" profile) that is mathematically well suited to this task,

$$\begin{eqnarray}&&{I}_{\nu }(\varrho )\propto {\left(\displaystyle \frac{\varrho }{{\varrho }_{t}}\right)}^{-\gamma }{\left[1+{\left(\displaystyle \frac{\varrho }{{\varrho }_{t}}\right)}^{\alpha }\right]}^{(\gamma -\beta )/\alpha },\end{eqnarray} \tag{ 1 }$$

where ϱ is the radial coordinate projected on the sky. The Nuker profile has five free parameters: (1) a transition radius ϱ_t, (2) an inner disk index γ, (3) an outer disk index β, (4) a transition index α, and (5) a normalization constant, which we recast to be the total flux density F_ν (defined so that ${F}_{\nu }=2\pi \int {I}_{\nu }(\varrho )\varrho \,d\varrho$).

It is instructive to consider some asymptotic behavior. When $\varrho \ll {\varrho }_{t}$ or $\varrho \gg {\varrho }_{t}$, the brightness profile scales like ${\varrho }^{-\gamma }$ or ${\varrho }^{-\beta }$, respectively. The index α controls where these asymptotic behaviors are relevant: a higher α value pushes the profile toward a sharp broken power-law morphology. The same parameterization can reproduce the emission profiles for standard, continuous disk models as well as the ring-like emission noted for "transition" disks (e.g., Andrews et al. 2011). Figure 2 shows Nuker profiles for some representative parameter values.

Zoom In Zoom Out Reset image size

From a given Nuker profile, we compute a set of model visibilities from the Fourier transform of Equation (1), sampled at the same discrete set of spatial frequencies as the data. Those visibilities are modified (stretched, rotated, and shifted) to account for four geometric parameters relevant to the observations: (1) a projected inclination angle i, (2) a rotation of the major axis in the plane of the sky φ (position angle), and (3) + (4) position offsets from the observed phase center {dα, dδ}. Taken together, nine parameters, θ = [ϱ_t, γ, β, α, ${F}_{\nu }$, i, φ, dα, dδ], fully specify a set of model visibilities, ${{ \mathcal V }}_{{\rm{m}}}(\theta )$.

It is worth highlighting two important points about the choice of this model prescription. First, the specific parameterization adopted for the brightness profile does not matter in the specific context of this study, so long as it accurately describes the data (i.e., leaves no statistically significant residuals). This will be illustrated more directly below, after a description of the full modeling procedure. Second, we have expressly avoided casting the interpretation in a physical context. Converting the results into inferences on optical depths or temperatures implicitly introduces a set of strong degeneracies and a model dependence that is not well motivated from a physical standpoint. Our aim is instead to hew as close to the observations as possible; we will consider some connections to physical properties only in the larger context of the results from the full sample (see Section 5).

To compare a given model to the visibility data, ${{ \mathcal V }}_{{\rm{d}}}$, we employ a standard Gaussian likelihood,

$$\begin{eqnarray}&&\mathrm{log}p({{ \mathcal V }}_{{\rm{d}}}| \theta )=-\displaystyle \frac{1}{2}\displaystyle \sum _{i}{w}_{{\rm{d}},i}| {{ \mathcal V }}_{{\rm{d}},i}-{{ \mathcal V }}_{{\rm{m}},i}(\theta ){| }^{2},\end{eqnarray} \tag{ 2 }$$

the sum of the residual moduli weighted by the standard natural visibility weights. The classic inference problem can be cast with the posterior probability distribution of the model parameters conditioned on the data as

$$\begin{eqnarray}&&\mathrm{log}p(\theta | {{ \mathcal V }}_{{\rm{d}}})=\mathrm{log}p({{ \mathcal V }}_{{\rm{d}}}| \theta )+\mathrm{log}p(\theta )+\mathrm{constant},\end{eqnarray} \tag{ 3 }$$

where $p(\theta )={\prod }_{j}p({\theta }_{j})$, the product of the priors for each parameter (presuming their independence).

We assign uniform priors for most parameters: $p({F}_{\nu })={ \mathcal U }$(0, 10 Jy), $p({\varrho }_{t})={ \mathcal U }$(0, 10''), $p(\varphi )={ \mathcal U }$(0, 180°), and $p(d\alpha )$, $p(d\delta )={ \mathcal U }$(−3, +3''). A simple geometric prior is adopted for the inclination angle, $p(i)=\sin i$.

Priors for the Nuker profile's index parameters are not obvious. They were assigned based on some iterative experimentation, grounded in our findings for the targets that have data with higher sensitivity and resolution. We set $p(\beta )={ \mathcal U }(2,10)$: the low bound forces the intensity profile to decrease at large radii (as is observed), and the high bound is practical—at these resolutions, the data are not able to differentiate between more extreme indices. To appropriately span both smooth and sharp transitions between the two power-law regimes, we set $p({\mathrm{log}}_{10}\alpha )\,={ \mathcal U }(0,2)$ (see Figure 2). We sample the posterior in ${\mathrm{log}}_{10}\alpha$ rather than α, since most of the diversity in Nuker profile shapes happens in the first decade of the prior-space. Finally, we set a softer-edged prior on γ that is approximately uniform over the range (−3, 2), but employs logistic tapers on the boundaries:

$$\begin{eqnarray}&&p(\gamma )\propto \displaystyle \frac{1}{1+{e}^{-5(\gamma +3)}}-\displaystyle \frac{1}{1+{e}^{-15(\gamma -2)}}.\end{eqnarray} \tag{ 4 }$$

The low-γ bound is again practical (like the high-β bound). The steeper bound at high γ values was designed to improve convergence for the poorly resolved cases: when $\gamma \gtrsim 2$, the degeneracy with ϱ_t is severe and both parameters can increase without bound (i.e., very steep gradients can be accommodated for very large transition radii). Since none of the well-resolved targets ended up having $\gamma \gtrsim 1$ (regardless of the γ prior), we consider this adopted upper boundary to be conservative.

A Markov chain Monte Carlo (MCMC) algorithm was used to explore the posterior probability space. We employed the ensemble sampler proposed by Goodman & Weare (2010) and implemented as the open-source code package emcee by Foreman-Mackey et al. (2013). With this algorithm, we used 48 "walkers" to sample in $\mathrm{log}p(\theta | {{ \mathcal V }}_{{\rm{d}}})$ and ran 50,000 steps per walker.

To identify a starting point for parameter estimation, we first crudely estimate the disk geometric parameters {i, φ, $d\alpha$, $d\delta$} with an elliptical Gaussian fit to the visibility data. Using these rough values, we deproject and azimuthally average the visibilities into a one-dimensional profile as a function of baseline length. This is used as a guide to visually probe models with different Nuker profile parameters, {F_ν, ϱ_t, α, β, γ}. When a reasonable match is found (with some experience, this takes about a minute), we assign a conservative (i.e., broad) range around each of the parameters. The ensemble sampler walkers are then initialized with random draws from uniform distributions that span those ranges.

While assessing formal convergence is difficult for this algorithm (since the walkers are co-dependent), we have used a variety of simple rubrics (e.g., trace and autocorrelation examinations, an interwalker test (Gelman & Rubin 1992), and the comparison of sub-chain means and variances advocated by Geweke 1992) that lend confidence that we have reached a stationary target distribution. Autocorrelation lengths for all parameters are of the order of 10² steps, implying that (after excising burn-in steps) we have ≳10⁴ independent samples of the posterior for each target. We find acceptance fractions of 0.2–0.4.

The disk size is not one of the parameters that we infer directly. If we followed the traditional methodology, we would proceed by taking ϱ_t as the size metric. But that approach is problematic; ϱ_t is not really how we think about an emission size. As an illustration, consider the two disk models shown in Figure 3 with identical ϱ_t but different γ values. Which is "larger"? Our instinct suggests the disk with a lower γ is larger, since its emission clearly has a more radially extended morphology.

Zoom In Zoom Out Reset image size

If we consider alternative prescriptions for the surface brightness profile, there is an analogous ambiguity. Any such model will include a "scale" parameter (such as ϱ_t) that will generally have a different value than the one inferred from the data for the Nuker profile. Given the limitations of the data imposed by noise and resolution constraints, we cannot definitively determine which of these profiles, and their corresponding size metrics, is most appropriate. As a more practical concern, a degeneracy with the parameters that describe the gradients of the Nuker profile,⁴ {α, β, γ}, makes it difficult to obtain a precise inference of ϱ_t using data with typical sensitivity and resolution (especially if $\alpha \lesssim 10$, since the resulting smoother profile makes the transition less distinct).

There is a generic definition of size that alleviates these issues. If we construct a cumulative intensity profile,

$$\begin{eqnarray}&&{f}_{\nu }(\varrho )=2\pi {\int }_{0}^{\varrho }{I}_{\nu }({\varrho }^{{\prime} })\,{\varrho }^{{\prime} }\,d{\varrho }^{{\prime} },\end{eqnarray} \tag{ 5 }$$

then we can assign an effective radius, ${\varrho }_{\mathrm{eff}}$, that encircles a fixed fraction, x, of the total flux: ${f}_{\nu }({\varrho }_{\mathrm{eff}})={{xF}}_{\nu }$ for some $x\in [0,1]$ (note that ${F}_{\nu }={f}_{\nu }(\infty )$ by definition). The inference of ${\varrho }_{\mathrm{eff}}$ is considerably less affected by imprecise constraints on the gradient parameters. Moreover, ${\varrho }_{\mathrm{eff}}$ more faithfully captures the intuitive intent of a size metric, as shown in Figure 3. Most importantly, it has the same straightforward meaning regardless of the underlying surface brightness profile; any profile that accurately reproduces the data will have the same ${\varrho }_{\mathrm{eff}}$.

That said, the selection of x in the definition of ${\varrho }_{\mathrm{eff}}$ is technically arbitrary. It makes sense to fix x and therefore homogenize the analysis for a sample. However, there are some practical concerns to take into consideration. If x is too low, then ${\varrho }_{\mathrm{eff}}$ would be uncomfortably reliant on a sub-resolution extrapolation of the Nuker profile. And if x is too high, then ${\varrho }_{\mathrm{eff}}$ would simply reflect the quality of the constraint on β (and in some cases α) based on the part of the brightness profile where we have the poorest sensitivity (large ϱ). Some of the more intuitive choices are x = 0.50 (so ${\varrho }_{\mathrm{eff}}$ is the "half-light" radius) or 0.68 (so ${\varrho }_{\mathrm{eff}}$ is comparable to a "standard deviation" in the admittedly poor approximation of a Gaussian brightness profile). Given the concern about extrapolation expressed above, we prefer to define ${\varrho }_{\mathrm{eff}}$ based on x = 0.68. An alternative choice in the range $x\in [0.5,0.8]$ makes little difference in the analysis or results that follow.

We derive posterior samples for ${\varrho }_{\mathrm{eff}}$ by integrating the cumulative intensity profiles (i.e., solving Equation (5)) for each brightness profile sampled by the θ posteriors. We think of ${\varrho }_{\mathrm{eff}}$ as a "distilled" size metric, since it reduces a complicated (and intrinsically uncertain) brightness profile into a straightforward, intuitive benchmark value.

To examine the relationship between size and luminosity, we need to work with physical (rather than observational) parameters (e.g., radii in au rather than arcseconds, and luminosities rather than flux densities) that account for both systematics in the flux calibration and the different target distances in the sample. Our procedure for the former is to multiply each posterior sample of F_ν by a factor s, drawn from a normal distribution $\sim { \mathcal N }(1.0,0.1)$, that mimics the ∼10% uncertainty in the absolute flux calibration of the data. For the latter, we start with a parallax measurement for each target and, presuming a normal probability distribution, $\sim { \mathcal N }(\varpi ,{\sigma }_{\varpi })$, we adopt the formalism of Astraatmadja & Bailer-Jones (2016) to infer a distance posterior, $p(d| \varpi ,{\sigma }_{\varpi })$, for a uniform distance prior.

Trigonometric parallaxes from the revised Hipparcos (van Leeuwen 2007) and Gaia DR1 (Gaia Collaboration et al. 2016) catalogs are available for AS 209, CQ Tau, UX Tau A, SU Aur, HD 163296, IM Lup, MWC 758, SAO 206462, and TW Hya. We assign the same parallax for HP Tau that was measured for its wide companion, HP Tau/G2, using very long baseline interferometry (VLBI) radio observations (Torres et al. 2009). The remaining 22 targets in Taurus are assigned parallaxes based on the values of their nearest neighbors. To do that, we compiled a list of 38 Taurus members with either Gaia DR1 or VLBI (Loinard et al. 2007; Torres et al. 2009) trigonometric parallaxes (A. L. Kraus 2017, private communication). For each target, we calculated the mean and standard deviation of the parallaxes from those members within a 5° radius. We associate WaOph 6 with its neighbor AS 209, albeit with an inflated ${\sigma }_{\varpi }$. For the remaining 14 targets in Ophiuchus, we adopt the mean ϖ measured from VLBI data for L1688 or L1689 sources by Ortiz-León et al. (2017). The parallax uncertainty for DoAr 33 was increased, since it lies north of L1688. LkHα 330, RX J1604.3−2130, and RX J1615.3−3255 are assigned parallaxes based on their affiliated cluster means (Per OB2, Upper Sco, and Lupus, respectively; see de Zeeuw et al. 1999; Galli et al. 2013).

For each posterior sample of {F_ν, ${\varrho }_{\mathrm{eff}}$}, we draw a distance from $p(d| \varpi ,{\sigma }_{\varpi })$ and use it to calculate a posterior sample of the logarithms of {${L}_{\mathrm{mm}}$, ${R}_{\mathrm{eff}}\}=\{{F}_{\nu }\times {(d/140)}^{2}\times s$, ${\varrho }_{\mathrm{eff}}\times d\}$, where d is in parsecs.⁵ The associated posteriors are summarized in Table 4.

To illustrate more concretely the procedure outlined above, we present a step-by-step analysis of the IQ Tau disk, which has a continuum luminosity that is roughly the median of the sample distribution.

We initialize the parameters (see Section 3.2) using ${F}_{\nu }\sim { \mathcal U }$(0.15, 0.19 Jy), ${\varrho }_{t}\sim { \mathcal U }$(0.3, 06), $\mathrm{log}\alpha \sim { \mathcal U }$(0.3, 1.3), $\beta \sim { \mathcal U }$(2, 6), $\gamma \sim { \mathcal U }$(−1, 1), $i\sim { \mathcal U }$(44, 58°), $\varphi \sim { \mathcal U }$(30, 45°), $d\alpha \sim { \mathcal U }$(0.01, 011), $d\delta \sim { \mathcal U }$(−0.39, −029). We then sample the posterior with the MCMC algorithm, and reach a stationary target distribution after 4000 steps. The walkers have acceptance fractions of 0.26–0.30; the autocorrelation times are 80–90 steps. Figure 4 shows the inferred covariances and marginalized posteriors for the Nuker profile parameters. We omit the geometric parameters for clarity, but they agree with independent estimates (e.g., Guilloteau et al. 2014). Confidence intervals for each parameter are derived from the marginalized posteriors and noted in Table 3. Figure 5 compares the simulated observations constructed with random draws from these posteriors to the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Surface Brightness Profile Parameters

Name	Fν (Jy)	ϱt (arcsec)	$\mathrm{log}\alpha$	β	γ	i (deg)	φ (deg)	${\varrho }_{\mathrm{eff}}$ (aecsec)
LkHα 330	0.208 ${}_{-0.002}^{+0.004}$	0.38 ${}_{-0.26}^{+0.45}$	0.83 ${}_{-0.10}^{+0.90}$	5.2 ${}_{-0.3}^{+2.7}$	−1.7 ${}_{-1.0}^{+0.3}$	39 ${}_{-5}^{+3}$	76 ${}_{-5}^{+7}$	0.44 ${}_{-0.01}^{+0.03}$
FN Tau	0.044 ${}_{-0.002}^{+0.002}$	$\lt 0.32$	$p(\mathrm{log}\alpha );\downarrow$	$\gt 3.5$	$p(\gamma );\uparrow$	p(i)	$p(\varphi )$	$\lt 0.17$
IRAS 04125+2902	0.042 ${}_{-0.001}^{+0.001}$	0.27 ${}_{-0.02}^{+0.03}$	$\gt 0.74$	6.8 ${}_{-1.0}^{+2.4}$	−2.5 ${}_{-0.3}^{+1.1}$	26 ${}_{-13}^{+7}$	155 ${}_{-26}^{+41}$	0.30 ${}_{-0.01}^{+0.02}$
CY Tau	0.207 ${}_{-0.005}^{+0.012}$	0.26 ${}_{-0.01}^{+0.16}$	$p(\mathrm{log}\alpha );\uparrow$	4.3 ${}_{-0.2}^{+3.7}$	0.3 ${}_{-2.6}^{+0.1}$	27 ${}_{-4}^{+3}$	136 ${}_{-8}^{+8}$	0.36 ${}_{-0.02}^{+0.02}$
V410 X-ray 2	0.026 ${}_{-0.001}^{+0.002}$	0.14 ${}_{-0.04}^{+0.38}$	$p(\mathrm{log}\alpha );\downarrow$	$p(\beta );\downarrow$	$p(\gamma );\uparrow$	62 ${}_{-27}^{+10}$	106 ${}_{-20}^{+19}$	0.16 ${}_{-0.02}^{+0.05}$
IP Tau	0.030 ${}_{-0.001}^{+0.002}$	0.25 ${}_{-0.04}^{+0.09}$	$p(\mathrm{log}\alpha );\uparrow$	$\gt 3.0$	−1.2 ${}_{-1.3}^{+1.4}$	41 ${}_{-20}^{+12}$	40 ${}_{-17}^{+93}$	0.27 ${}_{-0.04}^{+0.05}$
IQ Tau	0.163 ${}_{-0.002}^{+0.004}$	0.52 ${}_{-0.04}^{+0.39}$	0.40 ${}_{-0.13}^{+0.57}$	4.4 ${}_{-0.1}^{+4.5}$	0.6 ${}_{-0.5}^{+0.1}$	60 ${}_{-1}^{+1}$	43 ${}_{-2}^{+2}$	0.55 ${}_{-0.02}^{+0.01}$
UX Tau A	0.147 ${}_{-0.002}^{+0.001}$	0.27 ${}_{-0.01}^{+0.01}$	$\gt 1.25$	$\gt 8.0$	−2.9 ${}_{-0.4}^{+0.3}$	39 ${}_{-2}^{+2}$	167 ${}_{-2}^{+3}$	0.28 ${}_{-0.01}^{+0.01}$
DK Tau A	0.071 ${}_{-0.006}^{+0.004}$	0.08 ${}_{-0.01}^{+0.11}$	$p(\mathrm{log}\alpha );\downarrow$	2.6 ${}_{-0.1}^{+2.0}$	0.9 ${}_{-2.8}^{+0.1}$	26 ${}_{-12}^{+7}$	106 ${}_{-32}^{+24}$	0.32 ${}_{-0.02}^{+0.15}$
Haro 6-13	0.441 ${}_{-0.022}^{+0.069}$	0.06 ${}_{-0.01}^{+0.08}$	$p(\mathrm{log}\alpha );\uparrow$	2.3 ${}_{-0.1}^{+0.3}$	$p(\gamma );\uparrow$	42 ${}_{-3}^{+3}$	145 ${}_{-6}^{+6}$	0.48 ${}_{-0.05}^{+0.43}$
MHO 6	0.047 ${}_{-0.002}^{+0.002}$	0.23 ${}_{-0.05}^{+0.17}$	$p(\mathrm{log}\alpha );\downarrow$	$p(\beta );\uparrow$	$p(\gamma );\uparrow$	78 ${}_{-14}^{+7}$	110 ${}_{-8}^{+11}$	0.26 ${}_{-0.04}^{+0.06}$
FY Tau	0.023 ${}_{-0.001}^{+0.002}$	<0.50	$p(\mathrm{log}\alpha );\downarrow$	$p(\beta );\uparrow$	$p(\gamma );\uparrow$	$\gt 25$	134 ${}_{-62}^{+20}$	0.07 ${}_{-0.04}^{+0.04}$
UZ Tau E	0.401 ${}_{-0.007}^{+0.008}$	0.55 ${}_{-0.02}^{+0.01}$	0.83 ${}_{-0.08}^{+0.13}$	3.3 ${}_{-0.1}^{+0.2}$	0.6 ${}_{-0.1}^{+0.1}$	55 ${}_{-1}^{+1}$	90 ${}_{-1}^{+1}$	0.80 ${}_{-0.02}^{+0.03}$
J04333905+2227207	0.076 ${}_{-0.001}^{+0.002}$	0.58 ${}_{-0.04}^{+0.14}$	0.66 ${}_{-0.09}^{+0.83}$	$p(\beta );\uparrow$	0.3 ${}_{-0.3}^{+0.1}$	79 ${}_{-1}^{+1}$	115 ${}_{-2}^{+1}$	0.56 ${}_{-0.02}^{+0.02}$
J04334465+2615005	0.043 ${}_{-0.002}^{+0.002}$	0.46 ${}_{-0.15}^{+0.20}$	$p(\mathrm{log}\alpha );\downarrow$	$\gt 4.1$	1.2 ${}_{-1.2}^{+0.1}$	72 ${}_{-17}^{+11}$	158 ${}_{-7}^{+5}$	0.34 ${}_{-0.05}^{+0.04}$
DM Tau	0.254 ${}_{-0.008}^{+0.014}$	0.15 ${}_{-0.04}^{+0.07}$	$\lt 0.31$	3.0 ${}_{-0.3}^{+0.3}$	−1.4 ${}_{-1.2}^{+0.6}$	34${}_{-2}^{+2}$	158 ${}_{-5}^{+5}$	0.92 ${}_{-0.06}^{+0.15}$
CI Tau	0.440 ${}_{-0.014}^{+0.050}$	1.45 ${}_{-0.30}^{+0.69}$	0.38 ${}_{-0.20}^{+0.12}$	$\gt 4.1$	0.8 ${}_{-0.2}^{+0.1}$	44 ${}_{-3}^{+1}$	11 ${}_{-2}^{+3}$	0.89 ${}_{-0.04}^{+0.13}$
DN Tau	0.179 ${}_{-0.003}^{+0.005}$	0.35 ${}_{-0.02}^{+0.29}$	$p(\mathrm{log}\alpha );\downarrow$	$p(\beta );\downarrow$	0.8 ${}_{-0.8}^{+0.1}$	28 ${}_{-8}^{+2}$	80 ${}_{-9}^{+11}$	0.38 ${}_{-0.01}^{+0.02}$
HP Tau	0.106 ${}_{-0.002}^{+0.001}$	0.44 ${}_{-0.04}^{+0.05}$	$p(\mathrm{log}\alpha );\uparrow$	$\gt 3.8$	1.6 ${}_{-0.1}^{+0.1}$	64 ${}_{-3}^{+2}$	70 ${}_{-3}^{+3}$	0.22 ${}_{-0.01}^{+0.02}$
DO Tau	0.250 ${}_{-0.003}^{+0.012}$	0.12 ${}_{-0.01}^{+0.28}$	$p(\mathrm{log}\alpha );\downarrow$	${3.7}_{-0.7}^{+4.8}$	$p(\gamma );\uparrow$	37 ${}_{-7}^{+4}$	159 ${}_{-10}^{+7}$	0.19 ${}_{-0.01}^{+0.01}$
LkCa 15	0.396 ${}_{-0.004}^{+0.005}$	0.62 ${}_{-0.02}^{+0.05}$	0.62 ${}_{-0.10}^{+0.08}$	5.3 ${}_{-0.3}^{+1.0}$	−1.6 ${}_{-0.4}^{+0.2}$	51 ${}_{-1}^{+1}$	62 ${}_{-1}^{+1}$	0.77 ${}_{-0.01}^{+0.01}$
IRAS 04385+2550	0.061 ${}_{-0.003}^{+0.004}$	$p({\varrho }_{t});\downarrow$	$p(\mathrm{log}\alpha );\downarrow$	$p(\beta )$	>1.4	p(i)	148 ${}_{-19}^{+6}$	0.10 ${}_{-0.01}^{+0.08}$
GO Tau	0.179 ${}_{-0.007}^{+0.019}$	0.89 ${}_{-0.05}^{+2.73}$	0.20 ${}_{-0.08}^{+1.05}$	$p(\beta );\downarrow$	1.1 ${}_{-0.1}^{+0.1}$	53 ${}_{-3}^{+2}$	26 ${}_{-3}^{+4}$	1.16 ${}_{-0.07}^{+0.23}$
GM Aur	0.632 ${}_{-0.008}^{+0.011}$	0.92 ${}_{-0.27}^{+0.52}$	$\lt 0.26$	6.7 ${}_{-1.1}^{+2.3}$	−1.1 ${}_{-0.7}^{+0.4}$	55 ${}_{-1}^{+1}$	64 ${}_{-1}^{+1}$	0.87 ${}_{-0.02}^{+0.03}$
SU Aur	0.052 ${}_{-0.001}^{+0.002}$	$\lt 0.27$	$p(\mathrm{log}\alpha );\downarrow$	$\gt 3.5$	$p(\gamma );\uparrow$	p(i)	79 ${}_{-42}^{+57}$	$\lt 0.16$
V836 Tau	0.061 ${}_{-0.002}^{+0.002}$	0.13 ${}_{-0.01}^{+0.05}$	$p(\mathrm{log}\alpha );\uparrow$	$\gt 4.7$	$p(\gamma )$	61 ${}_{-8}^{+11}$	137 ${}_{-6}^{+9}$	0.13 ${}_{-0.01}^{+0.02}$
MWC 758	0.177 ${}_{-0.001}^{+0.002}$	0.50 ${}_{-0.01}^{+0.01}$	$\gt 1.3$	$\gt 7.9$	−3.0 ${}_{-0.4}^{+0.4}$	40 ${}_{-1}^{+1}$	168 ${}_{-1}^{+1}$	0.53 ${}_{-0.01}^{+0.01}$
CQ Tau	0.444 ${}_{-0.009}^{+0.013}$	0.34 ${}_{-0.01}^{+0.03}$	$\gt 0.87$	$\gt 6.6$	−2.3 ${}_{-0.5}^{+1.1}$	36 ${}_{-3}^{+3}$	53 ${}_{-6}^{+6}$	0.36 ${}_{-0.02}^{+0.01}$
TW Hya	1.311 ${}_{-0.002}^{+0.003}$	0.98 ${}_{-0.02}^{+0.01}$	0.95 ${}_{-0.03}^{+0.05}$	$\gt 8.6$	0.6 ${}_{-0.1}^{+0.1}$	7 ${}_{-1}^{+1}$	155 ${}_{-1}^{+1}$	0.77 ${}_{-0.01}^{+0.01}$
SAO 206462	0.596 ${}_{-0.022}^{+0.024}$	0.57 ${}_{-0.04}^{+0.02}$	0.88 ${}_{-0.08}^{+0.59}$	$\gt 5.9$	−1.0 ${}_{-0.7}^{+0.1}$	27 ${}_{-4}^{+3}$	33 ${}_{-9}^{+7}$	0.56 ${}_{-0.01}^{+0.02}$
IM Lup	0.587 ${}_{-0.006}^{+0.007}$	2.25 ${}_{-0.42}^{+0.48}$	0.32 ${}_{-0.07}^{+0.12}$	$\gt 5.5$	0.9 ${}_{-0.1}^{+0.1}$	49 ${}_{-2}^{+1}$	142 ${}_{-3}^{+2}$	1.11 ${}_{-0.02}^{+0.04}$
RX J1604.3−2130	0.194 ${}_{-0.006}^{+0.009}$	0.63 ${}_{-0.01}^{+0.01}$	$\gt 1.2$	5.8 ${}_{-0.4}^{+0.5}$	−2.5 ${}_{-0.3}^{+0.3}$	6 ${}_{-1}^{+1}$	77 ${}_{-1}^{+1}$	0.72 ${}_{-0.02}^{+0.01}$
RX J1615.3−3255	0.434 ${}_{-0.003}^{+0.005}$	1.39${}_{-0.51}^{+0.24}$	$\lt 0.31$	$\gt 5.1$	−0.1 ${}_{-0.3}^{+0.3}$	44 ${}_{-2}^{+2}$	149 ${}_{-3}^{+3}$	0.73 ${}_{-0.02}^{+0.02}$
SR 4	0.148 ${}_{-0.003}^{+0.003}$	0.19 ${}_{-0.01}^{+0.05}$	$\gt 0.45$	$\gt 5.2$	−0.9 ${}_{-1.5}^{+0.9}$	43 ${}_{-6}^{+3}$	27 ${}_{-6}^{+8}$	0.21 ${}_{-0.01}^{+0.01}$
Elias 20	0.264 ${}_{-0.004}^{+0.006}$	0.41 ${}_{-0.05}^{+0.05}$	$\gt 0.43$	7.0 ${}_{-1.5}^{+2.1}$	1.0 ${}_{-0.2}^{+0.1}$	54 ${}_{-2}^{+2}$	162 ${}_{-2}^{+2}$	0.32 ${}_{-0.01}^{+0.02}$
DoAr 25	0.565 ${}_{-0.008}^{+0.008}$	1.18 ${}_{-0.18}^{+0.09}$	0.58 ${}_{-0.07}^{+0.14}$	$\gt 5.5$	0.5 ${}_{-0.1}^{+0.1}$	63 ${}_{-1}^{+1}$	109 ${}_{-1}^{+1}$	0.85 ${}_{-0.02}^{+0.02}$
Elias 24	0.912 ${}_{-0.011}^{+0.016}$	0.82 ${}_{-0.03}^{+0.09}$	$\gt 0.76$	5.9 ${}_{-0.5}^{+2.2}$	1.1 ${}_{-0.1}^{+0.1}$	21 ${}_{-5}^{+3}$	62 ${}_{-11}^{+13}$	0.68 ${}_{-0.02}^{+0.01}$
GSS 39	0.654 ${}_{-0.010}^{+0.012}$	1.61 ${}_{-0.05}^{+0.05}$	$\gt 1.04$	$\gt 7.2$	1.1 ${}_{-0.1}^{+0.1}$	55 ${}_{-1}^{+1}$	120 ${}_{-1}^{+1}$	1.19 ${}_{-0.09}^{+0.02}$
WL 18	0.051 ${}_{-0.002}^{+0.004}$	0.13 ${}_{-0.03}^{+0.08}$	$p(\mathrm{log}\alpha );\uparrow$	$\gt 3.8$	$p(\gamma )$	60 ${}_{-29}^{+17}$	68 ${}_{-23}^{+54}$	0.13 ${}_{-0.02}^{+0.04}$
SR 24 S	0.509 ${}_{-0.006}^{+0.006}$	0.43 ${}_{-0.02}^{+0.05}$	0.65 ${}_{-0.08}^{+0.87}$	6.2 ${}_{-0.4}^{+2.9}$	−0.4 ${}_{-1.1}^{+0.1}$	46 ${}_{-2}^{+2}$	23 ${}_{-2}^{+2}$	0.45 ${}_{-0.08}^{+0.08}$
SR 21	0.405 ${}_{-0.003}^{+0.003}$	0.44 ${}_{-0.02}^{+0.02}$	$\gt 0.95$	$\gt 6.9$	−0.9 ${}_{-0.6}^{+0.2}$	18 ${}_{-9}^{+5}$	10 ${}_{-49}^{+45}$	0.43 ${}_{-0.01}^{+0.02}$
DoAr 33	0.070 ${}_{-0.004}^{+0.008}$	0.19 ${}_{-0.04}^{+0.29}$	$p(\mathrm{log}\alpha );\downarrow$	$p(\beta );\uparrow$	$p(\gamma );\uparrow$	68 ${}_{-27}^{+11}$	76 ${}_{-19}^{+21}$	0.21 ${}_{-0.03}^{+0.13}$
WSB 52	0.154 ${}_{-0.005}^{+0.016}$	0.93 ${}_{-0.24}^{+1.26}$	$p(\mathrm{log}\alpha );\downarrow$	$\gt 3.1$	1.6 ${}_{-0.3}^{+0.1}$	47 ${}_{-24}^{+9}$	165 ${}_{-26}^{+35}$	0.42 ${}_{-0.06}^{+0.12}$
WSB 60	0.258 ${}_{-0.006}^{+0.015}$	0.26 ${}_{-0.01}^{+0.01}$	0.37 ${}_{-0.02}^{+1.25}$	3.5 ${}_{-0.2}^{+3.2}$	−0.3 ${}_{-2.1}^{+0.1}$	32 ${}_{-12}^{+4}$	121 ${}_{-17}^{+17}$	0.34 ${}_{-0.02}^{+0.03}$
SR 13	0.153 ${}_{-0.003}^{+0.004}$	0.31 ${}_{-0.03}^{+0.03}$	$\gt 0.67$	$\gt 5.1$	−2.4 ${}_{-0.4}^{1.3}$	54 ${}_{-6}^{+4}$	72 ${}_{-7}^{+6}$	0.32 ${}_{-0.01}^{+0.03}$
DoAr 44	0.211 ${}_{-0.005}^{+0.006}$	0.38 ${}_{-0.03}^{+0.05}$	0.69 ${}_{-0.06}^{+0.87}$	5.6 ${}_{-0.4}^{+3.3}$	−1.8 ${}_{-0.8}^{+0.7}$	16 ${}_{-6}^{+10}$	55 ${}_{-27}^{+92}$	0.42 ${}_{-0.01}^{+0.02}$
RX J1633.9−2422	0.225 ${}_{-0.005}^{+0.011}$	0.34 ${}_{-0.03}^{+0.01}$	0.97 ${}_{-0.07}^{+0.80}$	6.7 ${}_{-0.7}^{+2.3}$	−1.5 ${}_{-0.9}^{+0.2}$	49 ${}_{-1}^{+1}$	83 ${}_{-1}^{+2}$	0.35 ${}_{-0.01}^{+0.01}$
WaOph 6	0.445 ${}_{-0.022}^{+0.024}$	0.47 ${}_{-0.03}^{+0.09}$	$\gt 0.33$	2.9 ${}_{-0.2}^{+0.5}$	1.1 ${}_{-0.1}^{+0.1}$	41 ${}_{-3}^{+4}$	173 ${}_{-3}^{+10}$	0.68${}_{-0.05}^{+0.15}$
AS 209	0.604 ${}_{-0.010}^{+0.007}$	1.65 ${}_{-0.50}^{+0.13}$	$\lt 0.16$	$\gt 6.4$	−0.2 ${}_{-0.1}^{+0.3}$	31 ${}_{-5}^{+3}$	73 ${}_{-9}^{+6}$	0.70 ${}_{-0.02}^{+0.02}$
HD 163296	1.822 ${}_{-0.006}^{+0.004}$	1.60 ${}_{-0.07}^{+0.04}$	0.54 ${}_{-0.04}^{+0.04}$	$\gt 8.6$	0.8 ${}_{-0.1}^{+0.1}$	47 ${}_{-1}^{+1}$	133 ${}_{-1}^{+1}$	0.96 ${}_{-0.01}^{+0.01}$

Note. Quoted values are the peaks of the posterior distributions, with uncertainties that represent the 68% confidence interval. Limits represent the 95% confidence boundary. A note of "p(X)" means the posterior is identical to the prior at <68% confidence; "$p(X);\downarrow$" or "$p(X);\uparrow$" means the posterior is consistent with the prior at 95% confidence, but has a marginal preference toward the lower or upper bound of the prior on X, respectively. Note that the F_ν posterior does not include systematic uncertainties in flux calibration (see Section 3.4) and that the ${\varrho }_{\mathrm{eff}}$ posterior is not directly inferred, but rather constructed from the joint posterior on {${F}_{\nu }$, ϱ_t, $\mathrm{log}\alpha$, β, γ}.

Download table as:  ASCIITypeset image

To construct posterior samples of ${\varrho }_{\mathrm{eff}}$, we calculate ${f}_{\nu }(\varrho )/{F}_{\nu }$ (see Equation (5)) for each posterior sample of θ. As an illustration, Figure 6 shows ${I}_{\nu }(\varrho )$ and ${f}_{\nu }(\varrho )/{F}_{\nu }$ for 200 random draws from the posteriors.

Zoom In Zoom Out Reset image size

If we circle back to our arguments advocating for a size metric like ${\varrho }_{\mathrm{eff}}$ (see Section 3.3), rather than ϱ_t (or its equivalent in other model prescriptions), IQ Tau provides a clear demonstration of how the former is more precise than the latter in the face of poorly constrained gradient parameters. Because {α, β, γ} are determined with relatively poor precision in this case, the 68% confidence interval on ϱ_t is (0.46, 090). In contrast, the inference on ${\varrho }_{\mathrm{eff}}$ is considerably narrower, (0.53, 056). The ambiguity in the indices does not matter much for ${\varrho }_{\mathrm{eff}}$: even for a diversity of ${I}_{\nu }(\varrho )$, the models that satisfactorily reproduce the data essentially have similar ${f}_{\nu }(\varrho )$.

To emphasize the point that the functional form of the brightness profile does not matter, we compared these Nuker profile inferences to the more standard models involving a broken power-law and a similarity solution. For the broken power-law prescription, the radius marking the transition between the two gradients lies in the range (0.40, 052); for the similarity solution, the radius marking the transition between the inner power law and the outer exponential taper is constrained to (0.22, 038). Despite these differences, all three models find statistically indistinguishable constraints on the effective size metric (${\varrho }_{\mathrm{eff}}$), with a 68% confidence interval (0.53, 056).

While the existence and qualities of a millimeter continuum size–luminosity relationship are now emerging, it is not particularly obvious why we might expect such behavior. The preliminary speculations on its potential origins can be differentiated into two broad categories: (1) the scaling is a natural consequence of the initial conditions, potentially coupled with some dispersion introduced by evolutionary effects; or (2) the scaling reflects some universal structural configuration. The concepts behind these categories are not always distinct, and certainly not mutually exclusive, but in light of the new data analysis presented here we will explore them separately in the following sections.

5.1.1. Initial Conditions, Evolutionary Dispersion

An appeal to a specific distribution of initial conditions may seem like a fine-tuning solution, but it cannot be easily dismissed given the strong links expected between young disks and the star formation process. Andrews et al. (2010) noted that the size–luminosity relationship is nearly perpendicular to the direction expected for viscously evolving disks that conserve angular momentum, and because of that they speculated that it may point to the underlying spread in the specific angular momenta in molecular cloud cores (see also Isella et al. 2009). A few other possibilities along these lines were raised by Piétu et al. (2014), including tidal truncation (unlikely for this sample of primarily single stars in low-density clusters; but see Testi et al. 2014) and significantly broad age and/or viscosity distributions. But neither of these studies considered how those initial conditions in the gas disks would be manifested a few million years later in observational tracers of the solid particles.

To assess this category of explanation for the size–luminosity relation, we considered the expected behavior of disks in the {${L}_{\mathrm{mm}}$, ${R}_{\mathrm{eff}}$} plane for simple models of the evolution of their solids when embedded in a standard, viscously evolving gas disk. We adopted the methodology of Birnstiel et al. (2015) to roughly estimate how the radial surface density profiles for solids of different sizes evolve with time for a range of initial conditions, parameterized through the total disk masses (M_d), (initial) characteristic radii (R_c), and (fixed) viscosity coefficients (${\alpha }_{t}$). In these calculations, we assumed a representative 0.5 M_⊙ host star and a fragmentation velocity of 10 m s⁻¹. For each time and combination of parameters from a coarse grid, we constructed a radial profile of 340 GHz optical depths as the summed (over particle size) product of the computed surface densities and representative opacities (the latter using the prescription of Ricci et al. 2010). We then converted that to an intensity profile with a crude approximation, ${I}_{\nu }(r)={B}_{\nu }[T(r)](1-{e}^{-{\tau }_{\nu }(r)})$, using a (fixed) representative temperature profile appropriate for the presumed stellar host (see Equation (9)). The luminosities and effective radii corresponding to those intensity profiles were calculated as in Sections 3.3 and 3.4.

Figure 11 shows a highly condensed summary of the models compared with the size and luminosity metrics inferred from the data. These models actually do a reasonably good job at reproducing the observed mean trend with a ${M}_{d}/{M}_{* }$ distribution from ∼0.5%–10%, a modest range of initial R_c of ∼60–150 au, and relatively low turbulence levels, ${\alpha }_{t}\lesssim 0.001$, but preferentially at early points in the evolutionary sequence (<1 Myr).

Zoom In Zoom Out Reset image size

While the ages of individual young stars are notoriously uncertain (e.g., Soderblom et al. 2014), the nominal cluster ages from which this sample is drawn are thought to be more like 2–3 Myr. By those times, the model predictions have diverged up and to the left in the {${L}_{\mathrm{mm}}$, ${R}_{\mathrm{eff}}$} plane (away from the mean relation, with larger ${R}_{\mathrm{eff}}$ than >90% of the sample targets for a given ${L}_{\mathrm{mm}}$) due to both viscous spreading and a relative depletion of particles that emit efficiently in the millimeter continuum. That discrepancy in timescale is not surprising; it is another manifestation of a generic feature of such models that assume a smooth gas disk structure (e.g., Takeuchi & Lin 2002, 2005; Brauer et al. 2007).

Putting aside the timescale problem, it is interesting that these models tend to predict a shallower slope than what is inferred from the data. This might suggest a correlation between the initial M_d and R_c, perhaps as a signature of the disk formation process. Alternatively, it might indicate an evolutionary timescale dependent on disk mass (perhaps a relationship between M_d and ${\alpha }_{t}$). Unfortunately, it would be difficult to disentangle such effects from the biases due to stellar age and mass likely present in such an inhomogeneously selected sample.

5.1.2. Small-scale, Optically Thick Substructures

The discrepancy in timescale highlighted above may be a telling failure of the underlying assumptions used to set up such models. The remedy often proposed for this inconsistency is the presence of fine-scale, localized maxima in the radial pressure profile of the gas disk (e.g., Whipple 1972; Pinilla et al. 2012). Those pressure peaks attract drifting particles, can slow or stop their inward motion, and thereby produce high solid concentrations in small areas that may promote rapid growth to larger bodies. Those concentrations would likely be manifested in the millimeter continuum as small regions of optically thick emission. Indeed, recent studies have found that narrow rings of bright emission accompanied by darker gaps, among other features, are prevalent in the few disks that have already been observed at very high spatial resolution (ALMA Partnership et al. 2015; Andrews et al. 2016; Cieza et al. 2016; Isella et al. 2016; Pérez et al. 2016; Zhang et al. 2016; Loomis et al. 2017). Such fine-scale optically thick features were considered by Piétu et al. (2014) as a potential contributor to a size–luminosity relation. Having firmly established the character of that relation here, it makes sense to revisit that possibility.

In Figure 12 we illustrate how optically thick emission could mimic the inferred distribution of disks in the {${L}_{\mathrm{mm}}$, ${R}_{\mathrm{eff}}$} plane. A fiducial optically thick surface brightness profile was constructed by setting ${I}_{\nu }(r)={B}_{\nu }[{T}_{d}(r)]$, where T_d(r) is specified from the approximation in Equation (9). We then consider a suite of such models where the disk extends out to some sharp cutoff radius, and for each of those profiles we calculate a luminosity and effective size as described in Section 3. For the median stellar luminosity of the sample (∼1 L_⊙), this results in the red curve in Figure 12. In this context, it makes sense that we do not see disks much to the right of that boundary (since they would all be completely optically thick and therefore could not produce more emission than allowed by the Planck function). If a given disk has an optically thin contribution or small-scale substructures that reduce the filling factor of optically thick emission below unity, it will lie to the left of this fiducial curve (as illustrated by the blue curve for the case of a 10% filling factor).

Zoom In Zoom Out Reset image size

In this scenario, the inferred distribution of disks in the size–luminosity plane could be explained if disks have optically thick filling fractions (inside a given ${R}_{\mathrm{eff}}$) of a few tens of percent. As one might expect, that is naturally in line with the rough estimate of the mean optical depth (≈1/3) inside ${R}_{\mathrm{eff}}$ implied by the inferred size–luminosity relation earlier in this section. Adopting again the IQ Tau disk as an example, we find that its location in the size–luminosity plane could be explained with an optically thick filling factor of ∼10%–30% inside 80 au.

Of course, the shift away from the pure optically thick curve for any disk can be achieved in myriad ways, depending on the detailed spatial distribution of the emission. That said, it is instructive to note that the mean size–luminosity relation can, in principle, be reproduced with a toy model of rings and gaps, as shown by the green curve in Figure 12. There, we have intentionally tuned the model to match the inferred relation by simply removing regions of emission from the purely thick model that span segments of a few au at regular spacings. The intent is to crudely mimic recent high-resolution images of disks. In this case, brighter disks represent systems where the local pressure maxima continue to larger distances from the host star (i.e., brighter disks are just systems that have more large rings). It is worth highlighting that one of the key hints that, for us, makes this explanation more than just plausible is that both "normal" and "transition" disks populate the same size–luminosity trend. The latter systems could be considered an especially simple form of substructure.

A robust assessment of the viability of this explanation will require measurements with much higher angular resolution than are yet available. But if such data end up demonstrating that optically thick substructures are important contributors to the continuum emission, the size–luminosity relation could end up being a useful (albeit indirect) demographic metric for understanding their origins in the broader disk population (especially because it is far less observationally expensive to reconstruct the {${L}_{\mathrm{mm}}$, ${R}_{\mathrm{eff}}$} plane than to measure fine-scale emission distributions for large samples).

If this explanation is relevant, it creates a series of other peripheral issues. For example, if optical depths are large, it means previous assessments of the disk masses might instead just reflect the underlying spatial distribution of substructures. It could imply that the inferences of particle growth from continuum spectral indices are contaminated; indeed, Ricci et al. (2012) already warned of essentially this same scenario. And it might indicate that the interpretation of optically thin spectral line emission from key tracer molecules near the inner disk midplane could be complicated (it could be obscured by the optically thick continuum, depending on its vertical distribution; e.g., Öberg et al. 2015; Cleeves et al. 2016).

One concern that limits our conclusions stems from the inhomogeneous sample of targets available in the SMA archive. This sample is notable for its size, but is composed of disks in the bright half of the millimeter luminosity distribution, which belong to various clusters (and thereby potentially different ages/environments) and consider a limited range of host masses. It would be beneficial to repeat this analysis for targets with lower continuum luminosities and host masses, and those formed in the same environment at the same time (or nearly so). An investigation for an older cluster (say ∼5 Myr) would be particularly compelling. The evolutionary models described above predict that the size–luminosity relationship should flatten out: older disks have more radially concentrated solids, and so smaller particles that tend to emit less overall in the millimeter continuum, but over more extended distributions, end up dominating the behavior in the {${L}_{\mathrm{mm}}$, ${R}_{\mathrm{eff}}$} plane. In the substructure scenario, we might expect little change in the shape of the relation, although effects such as growth in the trapped regions could make it difficult to characterize any trend.

Another important point to keep in mind is that the relationship we have quantified here is expressly valid only at ∼340 GHz. While Piétu et al. (2014) appear to have tentatively found a similar correlation at ∼230 GHz, we know empirically that the spatial extents of continuum emission from disks tend to increase with the observing frequency (e.g., Pérez et al. 2012, 2015; Menu et al. 2014; Tazzari et al. 2016). The frequency-dependent variation in the size–luminosity relation that this behavior might imply could prove to be crucial in better understanding its origins. In addition to improving the sample quality, complementing the analysis with analogous measurements at a lower (≲100 GHz) frequency (albeit at higher angular resolution, given the observed tendency for more compact emission at lower frequencies) should also be viewed as highly desirable.