An ALMA Survey of Protoplanetary Disks in the $\sigma$ Orionis Cluster

TL;DR

This ALMA survey of 92 protoplanetary disks in the $\sigma$ Orionis cluster reveals rapid dust dissipation and the rarity of giant planet formation within a few million years, highlighting the influence of external photoevaporation.

Contribution

First high-sensitivity ALMA survey of dust and gas in $\sigma$ Orionis disks, providing new insights into disk evolution and planet formation timescales.

Findings

Only 11 disks have dust mass >10 Earth masses.

Dust masses decline with proximity to the O9 star.

External photoevaporation influences disk dissipation.

Abstract

The $\sigma$ Orionis cluster is important for studying protoplanetary disk evolution, as its intermediate age ($\sim$3-5 Myr) is comparable to the median disk lifetime. We use ALMA to conduct a high-sensitivity survey of dust and gas in 92 protoplanetary disks around $\sigma$ Orionis members with $M_{\ast}\gtrsim0.1 M_{\odot}$. Our observations cover the 1.33 mm continuum and several CO $J=2-1$ lines: out of 92 sources, we detect 37 in the mm continuum and six in $^{12}$CO, three in $^{13}$CO, and none in C$^{18}$O. Using the continuum emission to estimate dust mass, we find only 11 disks with $M_{\rm dust}\gtrsim10 M_{\oplus}$, indicating that after only a few Myr of evolution most disks lack sufficient dust to form giant planet cores. Stacking the individually undetected continuum sources limits their average dust mass to 5$\times$ lower than that of the faintest detected disk, supporting theoretical models that indicate rapid dissipation once disk clearing begins. Comparing the protoplanetary disk population in $\sigma$ Orionis to those of other star-forming regions supports the steady decline in average dust mass and the steepening of the $M_{\rm dust}$-$M_{\ast}$ relation with age; studying these evolutionary trends can inform the relative importance of different disk processes during key eras of planet formation. External photoevaporation from the central O9 star is influencing disk evolution throughout the region: dust masses clearly decline with decreasing separation from the photoionizing source, and the handful of CO detections exist at projected separations $>1.5$ pc. Collectively, our findings indicate that giant planet formation is inherently rare and/or well underway by a few Myr of age.

Figures

c in Eistrup et al., [2016](#bib.bib35)); these reactions have typical timescales of a few Myr (depending on the ionization rate) and thus could be more significant in older systems. Finally, volatile carbon may be locked up in large icy bodies in the disk midplane (Bergin et al., [2010](#bib.bib18); Ros & Johansen, [2013](#bib.bib98); Guidi et al., [2016](#bib.bib46)). These large pebbles cannot diffuse upward and thus would no longer participate in gas-phase chemistry (Du et al., [2015](#bib.bib34); Kama et al., [2016](#bib.bib64)). Such a process would “dry out” the CO from the warm molecul

Tables

Table 1. Table 1 Continuum Properties

Source	RAJ2000	DecJ2000	SpT	Ref	$M_{\ast }$	$F_{1.33\mathrm{mm}}$	rms	$M_{\mathrm{dust}}$
					($M_{\odot }$)	(mJy)	(mJy beam-1)	($M_{\oplus }$)
1036	05:39:25.206	-02:38:22.09	K7.5 $\pm$ 0.5	H14	0.67 $\pm$ 0.06	5.94 $\pm$ 0.15	0.16	26.45 $\pm$ 0.66
1050	05:39:26.330	-02:28:37.70	M3.9 $\pm$ 2.0	VJ	0.19 $\pm$ 0.12	-0.15 $\pm$ 0.15	0.15	-0.67 $\pm$ 0.66
1075	05:39:29.350	-02:27:21.02	M0.0 $\pm$ 1.5	H14	0.62 $\pm$ 0.14	1.48 $\pm$ 0.15	0.16	6.57 $\pm$ 0.65
1152	05:39:39.377	-02:17:04.50	M0.0 $\pm$ 1.0	H14	0.62 $\pm$ 0.16	8.57 $\pm$ 0.17	0.18	38.16 $\pm$ 0.77
1153	05:39:39.828	-02:31:21.89	K5.5 $\pm$ 1.0	H14	0.91 $\pm$ 0.12	13.62 $\pm$ 0.16	0.18	60.66 $\pm$ 0.72
1154	05:39:39.833	-02:33:16.08	M3.8 $\pm$ 2.0	VJ	0.25 $\pm$ 0.14	1.44 $\pm$ 0.15	0.16	6.43 $\pm$ 0.65
1155	05:39:39.900	-02:43:09.00	K1.0 $\pm$ 2.5	H14	1.71 $\pm$ 0.25	-0.12 $\pm$ 0.15	0.15	-0.54 $\pm$ 0.65
1156	05:39:40.171	-02:20:48.04	K5.0 $\pm$ 1.0	H14	0.96 $\pm$ 0.16	5.66 $\pm$ 0.15	0.18	25.21 $\pm$ 0.68
1182	05:39:43.190	-02:32:43.30	M4.0 $\pm$ 2.0	VJ	0.20 $\pm$ 0.13	0.18 $\pm$ 0.15	0.15	0.82 $\pm$ 0.66
1193	05:39:44.510	-02:24:43.20	M5.0 $\pm$ 2.0	VJ	0.10 $\pm$ 0.10	-0.02 $\pm$ 0.15	0.15	-0.09 $\pm$ 0.66

Table 2. Table 2 Gas Properties

Source	$F_{12CO}$	$F_{13CO}$	$F_{\mathrm{C18O}}$	$M_{\mathrm{gas}}$	$M_{\mathrm{gas},\min }$	$M_{\mathrm{gas},\max }$
	(mJy km s-1)	(mJy km s-1)	(mJy km s-1)	($M_{\mathrm{Jup}}$)	($M_{\mathrm{Jup}}$)	($M_{\mathrm{Jup}}$)
Gas Detections
540	1204 $\pm$ 85	276 $\pm$ 54		2.4	1.0	10.5
1274	861 $\pm$ 88	326 $\pm$ 68		5.5	1.0	31.4
1152	633 $\pm$ 82	314 $\pm$ 65		7.1	1.0	31.4
1153	557 $\pm$ 57			…	…	1.0
818	514 $\pm$ 58			…	…	1.0
1075	165 $\pm$ 33			…	…	10.5
Gas Non-detections
1036				…	…	3.1
1050				…	…	3.1
1154				…	…	3.1
1155				…	…	3.1
1156				…	…	3.1
1182				…	…	3.1
1193				…	…	3.1
1230				…	…	3.1
1248				…	…	3.1
1260				…	…	3.1

Table 3. Table 3 Gas Properties derived from Miotello et al. ( 2016a ) models

Source	$M_{\mathrm{gas}}$ ($M_{\mathrm{Jup}}$)	Gas-to-dust ratio	Line
540			C18O
1274			C18O
1152			C18O
1153			13CO
818			13CO
1075			13CO

Table 4. Table 4 M dust subscript 𝑀 dust M_{\rm dust} – M ∗ subscript 𝑀 ∗ M_{\ast} Bayesian Fit Parameters

Region	Age (Myr)	${\alpha }^{‡}$	${\beta }^{‡}$	$\delta$
Taurus	1–2	1.2$\pm$0.1	1.7$\pm$0.2	0.7$\pm$0.1
Lupus†	1–3	1.2$\pm$0.2	1.8$\pm$0.4	0.9$\pm$0.1
Cha I	2–3	1.0$\pm$0.1	1.8$\pm$0.3	0.8$\pm$0.1
$\sigma$ Orionis	3–5	1.0$\pm$0.2	2.0$\pm$0.4	0.6$\pm$0.1
Upper Sco	5–11	0.8$\pm$0.2	2.4$\pm$0.4	0.7$\pm$0.1

Equations

$$M_{\mathrm{dust}}=\frac{F_{\nu }{d}^2}{{\kappa }_{\nu }{B}_{\nu }(T_{\mathrm{dust}})}.$$

Full text

††slugcomment: Accepted to AJ: 23 March 2017

I. INTRODUCTION

Planets are thought to form in so-called “protoplanetary” disks around young stars within $\sim$5–10 Myr (Williams & Cieza, 2011). The resulting exoplanet population is diverse, as revealed by the Kepler transit survey (Borucki et al., 2010) and long-term radial velocity (RV) surveys (Howard et al., 2010; Mayor et al., 2011). However, certain trends are emerging; for example, intermediate-mass planets (i.e., “super-Earths” with masses between that of Earth and Neptune) appear to be an order of magnitude more abundant than gas giants (i.e., planets with masses several times that of Jupiter), at least for short orbital periods (Howard et al., 2012; Fressin et al., 2013; Petigura et al., 2013; Gaidos et al., 2016). To identify the origins of these trends, and thus better understand planet formation, we must survey the preceding protoplanetary disks. Indeed, exoplanet population synthesis models indicate that planetary properties and architectures are dictated by the initial dust and gas content of protoplanetary disks and their subsequent evolution (e.g., Thommes et al., 2008; Mordasini et al., 2012; Bitsch et al., 2015; Mordasini et al., 2016).

Sub-mm and mm wavelength surveys are particularly useful for probing the bulk dust and gas content of protoplanetary disks, as disk emission at these longer wavelengths can be optically thin. The first (sub-)mm surveys of star-forming regions made the initial steps in identifying trends in protoplanetary disk evolution that could potentially explain correlations seen in the exoplanet population (Andrews et al., 2009; Lee et al., 2011; Andrews et al., 2013; Williams et al., 2013; Ansdell et al., 2015). Most notably, early surveys of Taurus disks (Natta et al., 2000; Andrews et al., 2013) revealed a positive dependence between disk dust mass ($M_{\rm dust}$) and host star mass ($M_{\ast}$), which could fundamentally explain the positive correlation between giant planet frequency and stellar mass (Endl et al., 2006; Johnson et al., 2007; Bowler et al., 2010; Bonfils et al., 2013). However, these initial disk surveys were often incomplete and limited by dust mass sensitivities of a few Earth masses. These constraints meant that it remained unclear whether (sub-)mm continuum emission systematically declines with age, reflecting steady disk dispersal and/or grain growth in protoplanetary disks (e.g., Williams, 2012). Moreover, none of these initial surveys probed bulk gas mass, as contemporary facilities lacked the sensitivity to detect faint line emission.

Measuring both dust and gas content independently is essential for studying planet formation, as growing dust grains decouple from the gas and evolve differently, yet both components determine what types of planets may form in a disk. However, due to the challenges associated with estimating disk gas masses, the canonical interstellar medium (ISM) gas-to-dust ratio of $\sim$100 (Bohlin et al., 1978) is often used to infer gas mass from dust mass, requiring an extrapolation of two orders of magnitude. Recent observations suggest that the inherited ISM ratio may actually decrease by an order of magnitude after just a few Myr of evolution (Williams & Best, 2014; Ansdell et al., 2016), although these calculations may be underestimated due to carbon depletion (e.g., Miotello et al., 2016b). If gas is being depleted in disks (e.g., due to winds) then this may help to explain the lack of gas giants and prevalence of super-Earths seen in the exoplanet population (e.g., Fressin et al., 2013; Gaidos et al., 2016). In this scenario, super-Earths would result when giant planet cores form in gas-depleted disks, prohibiting the cores from rapidly accreting gaseous envelopes (e.g., Lee & Chiang, 2016) as predicted by core accretion theory (e.g., Pollack et al., 1996; Ida & Lin, 2004).

The enhanced sensitivity of the Atacama Large Millimeter/sub-millimeter Array (ALMA) now allows for efficient surveys of both dust and gas for large samples of protoplanetary disks across star-forming regions spanning the expected disk lifetime ($\sim$1–10 Myr). The first large-scale protoplanetary disk surveys conducted with ALMA include: Ansdell et al. (2016), who carried out a near-complete survey of 89 protoplanetary disks in the young Lupus star-forming region ($\sim$1–3 Myr at $\sim$150 pc; Comerón, 2008; Alcalá et al., 2014) with continuum and line sensitivities corresponding to M_{\rm dust}$$\sim$$0.3 $M_{\oplus}$ and M_{\rm gas}$$\sim$$1.0 $M_{\rm Jup}$, respectively; Barenfeld et al. (2016), who observed 106 disks in the more evolved Upper Sco region ($\sim$5–10 Myr at 145 pc; Preibisch et al., 2002; Pecaut et al., 2012) with sensitivities of M_{\rm dust}$$\sim$$0.1 $M_{\oplus}$; and Pascucci et al. (2016), who surveyed 93 protoplanetary disks in the young Chamaeleon I region ($\sim$2–3 Myr at 160 pc; Luhman, 2008) with sensitivities of M_{\rm dust}$$\sim0.2–0.8 $M_{\oplus}$.111These dust mass sensitivities correspond to $3\times$ the rms and use the assumptions described in Section V.1 to convert (sub-)mm flux to dust mass.

These ALMA disk surveys are beginning to reveal trends in protoplanetary disk evolution that can help to constrain planet formation theory and explain correlations seen in the exoplanet population. One of their clearest findings is that average disk dust mass does indeed decrease with age. Ansdell et al. (2016) showed that Lupus disks have a mean dust mass $\sim$3$\times$ higher than that of the older Upper Sco region, but are statistically indistinguishable from disks in the similarly aged Taurus region. Equivalently, Barenfeld et al. (2016) found that the average $M_{\rm dust}/M_{\ast}$ ratio in Upper Sco is $\sim$4.5$\times$ lower than in Taurus. Interestingly, even in the younger regions, most disks lack sufficient dust to form the solid cores needed to build gas giants (e.g., only 26% of protoplanetary disks in Lupus have $M_{\rm dust}\gtrsim 10~{}M_{\oplus}$; Ansdell et al., 2016). These findings point to significant global disk evolution during the first few Myr and imply that giant planet formation occurs rapidly and/or is rare.

Moreover, these ALMA disk surveys have confirmed the aforementioned $M_{\rm dust}$–$M_{\ast}$ relation initially seen in the pre-ALMA surveys of Taurus, while also revealing a steepening of the relation with age (Ansdell et al., 2016; Pascucci et al., 2016). This steepening would indicate that dust evolution occurs more rapidly around lower-mass stars, and can be compared to theoretical models to constrain the relative importance of different disk evolution processes (Pascucci et al., 2016). Finally, Manara et al. (2016) have combined estimates of disk mass ($M_{\rm disk}$) from ALMA with spectroscopic measurements of stellar mass accretion rates ($\dot{M}_{\rm acc}$) from VLT/X-Shooter in the Lupus star-forming region to provide the first observational confirmation of the $M_{\rm disk}$–$\dot{M}_{\rm acc}$ relation predicted by viscous disk evolution theory.

The $\sigma$ Orionis region (Walter et al., 2008) is a particularly important target for studying disk evolution due to its intermediate age of $\sim$3–5 Myr (Oliveira et al., 2002, 2004), which is comparable to the median disk lifetime (Williams & Cieza, 2011). Only one-third of cluster members (92 sources) exhibit strong infrared (IR) excess indicative of a protoplanetary disk (Hernández et al., 2007). Williams et al. (2013) surveyed the region at 850 $\mu$m with JCMT/SCUBA-2, detecting just 8 disks at \sim$$4~{}M_{\oplus} sensitivity; they also noted the remarkable diversity in the spectral energy distributions (SEDs) of their detections (see their Figure 3), indicating substantial and ongoing disk evolution. By stacking their individually non-detected disks, Williams et al. (2013) found a mean signal of 1.3 mJy at 4$\sigma$ significance, motivating more sensitive follow-up observations of this region.

We therefore use ALMA to conduct a high-sensitivity mm wavelength survey of all known protoplanetary disks in $\sigma$ Orionis in both dust and gas. We describe the sample in Section II and our ALMA observations in Section III. The continuum and line measurements are presented in Section IV, then converted to dust and gas masses in Section V. We interpret our findings within the context of disk evolution in Section VI by identifying correlations with stellar and cluster properties as well as comparing our results to those found in other star-forming regions. This work is summarized in Section VII.

II. SAMPLE

The $\sigma$ Orionis cluster consists of several hundred young stellar objects (YSOs) ranging from brown dwarfs to OB-type stars (see review in Walter et al., 2008). The cluster is named after its brightest member, $\sigma$ Ori, a trapezium-like system whose most massive component is an O9 star. Cluster membership is defined by the Mayrit catalog (Caballero, 2008), which identifies 241 stars and brown dwarfs that are located within 30 arcmin of the $\sigma$ Ori system and exhibit known features of youth (X-ray emission, Li 6708Å absorption, etc.). We adopt a cluster distance of 385 pc based on recent orbital parallax measurements of the $\sigma$ Ori triple system (Schaefer et al., 2016). The low reddening towards this cluster, estimated at $E(B-V)\sim 0.05$ (e.g. Sherry et al., 2008), makes it a valuable site for analyzing the evolution of young stars.

Our sample consists of the 92 YSOs in $\sigma$ Orionis with IR excesses consistent with the presence of a protoplanetary disk. These sources are identified by cross-matching the Class II and transition disk (TD) candidates from the Spitzer survey of Hernández et al. (2007) with the aforementioned Mayrit catalog (Caballero, 2008). Both catalogs are expected to be complete down to the brown dwarf limit. Disk classifications are based on the Spitzer/IRAC SED slope, as described in Hernández et al. (2007). We also include in our sample a Class I disk (1153), as it is located near the Spitzer/IRAC color cutoff for Class II disks.

The sources in our sample are presented in Table 1 with their stellar spectral types (SpT) and stellar masses ($M_{\ast}$). Spectral types were primarily taken from the homogenous sample of low-resolution optical spectra analyzed in Hernández et al. (2014), but supplemented with those from medium-resolution VLT/X-Shooter spectra when available from Rigliaco et al. (2012). For the 23 sources that lack spectroscopic information, we estimate their spectral types using an empirical relation between $V-J$ color and stellar spectral type; the relation was derived by measuring synthetic photometry from flux-calibrated VLT/X-Shooter spectra of YSOs with spectral types from G5 to M9.5, then performing a non-parametric fit of the $V-J$ color versus spectral type relation (Manara et al. 2017, submitted). For these sources with photometrically derived spectral types, we cautiously assume uncertainties of $\pm$ two spectral subtypes. Figure 1 shows the stellar spectral type distribution of our sample.

We estimate $M_{\ast}$ values for our sample by comparing their positions on the Hertzsprung-Russel (HR) diagram to the evolutionary models of Siess et al. (2000). In order to place our targets on the HR diagram, we convert their spectral types to stellar effective temperatures ($T_{\rm eff}$) and derive their stellar luminosities ($L_{\ast}$) from $J$-band magnitudes using the relations in Herczeg & Hillenbrand (2015). The uncertainties on $L_{\ast}$ are obtained by propagating the uncertainties on spectral type and bolometric correction, and thus on distance and optical extinction ($A_{V}$). We then calculate the uncertainties on $M_{\ast}$ using a Monte Carlo (MC) method, where we take the standard deviation of 1000 estimates of $M_{\ast}$, each calculated after randomly perturbing the derived values of $T_{\rm eff}$ and $L_{\ast}$ by their uncertainties.

III. ALMA OBSERVATIONS

Our Band 6 ALMA observations were obtained on 2016 July 30 and 31 during Cycle 3 (Project ID: 2015.1.00089.S; PI: Williams). The array configuration used 36 and 37 12 m antennas on July 30 and 31, respectively, with baselines of 15–1124 m on both runs. The correlator setup included two broadband continuum windows centered on 234.293 and 216.484 GHz with bandwidths of 2.000 and 1.875 GHz and channel widths of 15.625 and 0.976 MHz, respectively. The bandwidth-weighted mean continuum frequency was 225.676 GHz (1.33 mm). The spectral windows covered the 12CO (230.538 GHz), 13CO (220.399 GHz), and C18O (219.560 GHz) $J=2$–1 transitions at velocity resolutions of 0.16-0.17 km s*-1*. These spectral windows were centered on 230.531, 220.392, and 219.554 GHz with bandwidths of 11.719 MHz and channel widths of 0.122 MHz.

On-source integration times were 1.2 min per object for an average continuum rms of 0.15 mJy beam*-1* (Table 1). This sensitivity was based on the JCMT/SCUBA-2 survey of $\sigma$ Orionis disks by Williams et al. (2013), who found that stacking their individual non-detections revealed a mean 850 $\mu$m continuum signal of 1.3 mJy at 4$\sigma$ significance. The sensitivity of our ALMA survey was therefore chosen to provide $\sim$3–4$\sigma$ detections of such disks at 1.3 mm, based on an extrapolation of the 850 $\mu$m mean signal using a spectral slope of $\alpha=2$–3.

The raw data were pipeline calibrated by NRAO staff using the CASA package (version 4.5.3). The pipeline calibration included: absolute flux calibration with observations of J0522-3627 or J0423-0120; bandpass calibration with observations of J0510+1800 or J0522-3627; and gain calibration with observations of J0532-0307. We estimate an absolute flux calibration error of $\sim$10% based on the amplitude variations of gain calibrators over time.

We extract the continuum images from the calibrated visibilities by averaging over the continuum channels and cleaning with a Briggs robust weighting parameter of +0.5 for an average beam size of 0$\farcs$31$\times$0$\farcs$25 ($\sim$120$\times$95 AU at 385 pc). We extract 12CO, 13CO, and C18O $J=2$–1 line channel maps from the calibrated visibilities by subtracting the continuum from the spectral windows containing line emission using the uvcontsub routine in CASA. Sources showing clear line emission were cleaned with a Briggs robust weighting parameter of +0.5. We find average rms values of 13, 14, and 11 mJy beam*-1* within 1 km s*-1* velocity channels for the 12CO, 13CO, and C18O lines, respectively.

IV. ALMA RESULTS

IV.1. 1.33 mm Continuum Emission

Nearly all sources in our sample are unresolved, thus we measure continuum flux densities by fitting point-source models to the visibility data using the uvmodelfit routine in CASA. The point-source model has three free parameters: integrated flux density ($F_{\lambda}$), right ascension offset from the phase center ($\Delta\alpha$), and declination offset from the phase center ($\Delta\delta$). For the five resolved sources (1036, 1152, 1153, 1274, 540), we fit an elliptical Gaussian model instead, which has three additional free parameters: FWHM along the major axis ($a$), aspect ratio of the axes ($r$), and position angle (PA). We scale the uncertainties on the fitted parameters by the square root of the reduced $\chi^{2}$ value of the fit.

Table 1 presents the 1.33 mm continuum flux densities and associated uncertainties ($F_{\rm 1.33mm}$), where the uncertainties are statistical errors and do not include the 10% flux calibration error (Section III). We detect only 37 out of the 92 observed sources at $>3\sigma$ significance (Figure 2). For detections, the source locations in Table 1 are the fitted source centers output by uvmodelfit, while for non-detections they are simply the phase centers of the ALMA observations, which were chosen based on 2MASS positions. The average offsets from the phase centers for the detections are $\Delta\alpha=0\farcs 057$ and $\Delta\delta=-0\farcs 096$ (1.9 and $-3.2$ pixels), both much smaller than the average beam size (Section III). We also note that only 5 out of the 37 continuum detections have photometrically derived spectral types, which are less precise than the spectroscopically determined spectral types (Section II).

IV.2. CO Line Emission

To search for objects exhibiting significant line emission, we first extract the 12CO spectrum for each source. When creating the spectrum, we use 1 km s*-1* velocity sampling and measure fluxes in each channel using a circular aperture 0$\farcs$30 in radius and centered on the continuum emission (for detections) or the expected stellar position (for non-detections). We measure the image rms using a 4″–9″radius annulus centered on the fitted or expected source position. Candidate detections are identified as those with emission exceeding 3$\times$ the rms in multiple nearby channels within 0–25 km s*-1* (LSRK frame), which covers the range of RVs found for $\sigma$ Orionis members (e.g., Jeffries et al., 2006).

For each candidate detection, we create zero-moment maps by integrating across the velocity range where the emission exceeds the noise. The integrated flux ($F_{\rm 12CO}$) is then measured using circular aperture photometry, where the aperture radius for each source is determined by a curve-of-growth method in which successively larger apertures are applied until the flux stabilizes to within errors. Uncertainties ($E_{\rm 12CO}$) are estimated by taking the standard deviation of the fluxes measured within the same sized aperture placed randomly within the field of view but away from the source. We consider sources as detections when $F_{\rm 12CO}>4\times E_{\rm 12CO}$. We adopt this high detection threshold because, as pointed out in Barenfeld et al. (2016), this procedure selects both the velocity range and aperture size that maximize the signal, thus can produce false detections at lower significance levels. We detect only six sources in 12CO using this procedure. For these sources, we also search for 13CO and C18O emission using the same velocity range and aperture photometry method as for 12CO; we detect three of these sources in 13CO and none in C18O.

For sources with no significant line emission found using the above procedure, we create zero-moment maps by integrating across the channels $\pm$1 km s*-1* from their known RVs, when available in the literature (Maxted et al., 2008; Sacco et al., 2008). For sources with unknown RVs, we integrate around the average value for $\sigma$ Orionis members with known RVs. Jeffries et al. (2006) showed that $\sigma$ Orionis members are divided into two kinematically distinct subgroups differentiated by their RVs (Group 1 and 2, by their convention). However, the region is dominated by Group 2 sources at $\delta<-02$:18:00, where all our gas non-detections with unknown RVs are located. Thus we adopt the average RV of Group 2 (13 km s*-1* in the LSRK frame) when creating zero-moment maps for our gas non-detections with unknown RVs. We measure 12CO, 13CO, and C18O integrated fluxes from these zero-moment maps using the aforementioned aperture photometry method, but with an aperture size fixed to the beam size. We found no additional detections, thus took upper limits as 3$\times$ the image rms.

Table 2 gives our integrated line fluxes or upper limits. Of the 92 targets, only six are detected in 12CO, three are detected in 13CO, and none are detected in C18O with $>4\sigma$ significance. All sources detected in 12CO are detected in the continuum, and all sources detected in 13CO are detected in 12CO. The zero- and first-moment maps of the gas detections are shown in Figure 3.

V. PROPERTIES OF $\sigma$ ORIONIS DISKS

V.1. Dust Masses

Because dust emission at (sub-)mm wavelengths can be optically thin, the bulk dust mass of a disk ($M_{\rm dust}$) can be estimated from its (sub-)mm continuum emission at a given wavelength ($F_{\nu}$), as shown in Hildebrand (1983):

[TABLE]

where $B_{\nu}(T_{\rm dust})$ is the Planck function for a characteristic dust temperature of $T_{\rm dust}=20$ K (the median for Taurus disks; Andrews & Williams 2005). We take the dust grain opacity, $\kappa_{\nu}$, as 10 cm2 g*-1* at 1000 GHz and use an opacity power-law index of $\beta=1$ (Beckwith et al., 1990). The source distance, $d$, is taken as 385 pc based on the updated parallax of the $\sigma$ Ori triple system (Schaefer et al., 2016). Equation 1 can therefore be approximated as $M_{\rm dust}\approx 1.34\times 10^{-5}~{}F_{1.33{\rm mm}}$, where $F_{1.33{\rm mm}}$ is in mJy and $M_{\rm dust}$ is in $M_{\odot}$.

We use this simplified approach, which assumes a single grain opacity and an isothermal disk temperature, to ease comparisons with other disk surveys (Sections VI.2 & VI.3). Moreover, although Andrews et al. (2013) derived a $T_{\rm dust}=25{\rm K}\times(L_{\ast}/L_{\odot})^{0.25}$ relation from two-dimensional continuum radiative transfer models, we use an isothermal disk temperature as more detailed modeling of resolved disks suggests that $T_{\rm dust}$ is independent of stellar parameters. In particular, Tazzari et al. (2017, submitted) fit continuum observations of 36 resolved disks in Lupus directly in the uv-plane to two-layer disk models by solving energy balance equations at each disk radius, finding no dependence of $T_{\rm dust}$ as a function of stellar parameters.

Table 1 presents our $M_{\rm dust}$ estimates, derived using Equation 1 with our $F_{\rm 1.33mm}$ measurements (Section IV.1). Figure 4 shows the continuum-detected disks in order of increasing $M_{\rm dust}$ as well as the typical 3$\sigma$ upper limit of $\sim$2.0 $M_{\oplus}$ ($F_{\rm 1.33mm}\sim 0.45$ mJy). Only 4 disks have $M_{\rm dust}>30~{}M_{\oplus}$ ($F_{\rm 1.33mm}>6.7$ mJy), thus nearly all protoplanetary disks in $\sigma$ Orionis have dust masses well below the minimum mass of solids needed to form the planets in our Solar System (Weidenschilling, 1977). Moreover, only 11 disks have $M_{\rm dust}>10~{}M_{\oplus}$ ($F_{\rm 1.33mm}>2.2$ mJy), thus by $\sim$3–5 Myr most protoplanetary disks appear to lack sufficient dust to form giant planet cores. Note that significant amounts of solids may still exist in objects larger than a few cm in size, which do not produce detectable (sub-)mm emission.

V.2. Gas Masses

CO line emission can be used to roughly estimate bulk gas masses ($M_{\rm gas}$) independently from the dust, assuming simple CO chemistry and adopting an ISM-like CO/H2 abundance. Williams & Best (2014) (hereafter WB14) used parametrized gas disk models to demonstrate that the majority of CO may exist in the warm molecular layer, where it is sufficiently warm to survive freeze-out onto the disk midplane as well as adequately shielded from UV radiation to avoid photodissociation (except for particularly cold or low-mass disks). The CO isotopologue lines, specifically 13CO and C18O, are especially useful for constraining $M_{\rm gas}$ as their moderate-to-low optical depths mean that they trace the bulk gas content rather than the temperature profile of the disk. Ansdell et al. (2016) and Miotello et al. (2016b) have used these CO isotopologues lines to estimate $M_{\rm gas}$ for protoplanetary disks in the Lupus clouds that have been surveyed by ALMA.

Unfortunately, because we find no C18O detections in $\sigma$ Orionis, we cannot use the same combination of CO isotopologue lines to estimate $M_{\rm gas}$ in this region. However, we can still place rough constraints on $M_{\rm gas}$ by comparing our measured 12CO and 13CO line luminosities or upper limits to the WB14 model grid. The uncertainties on $M_{\rm gas}$ are larger for this line combination because 12CO is optically thick and therefore more sensitive to other disk parameters, such as the temperature profile. Nevertheless, $M_{\rm gas}$ can still be estimated using this method because the combination of integrated line fluxes still primarily depends on bulk gas mass rather than these other disk parameters (see the parameter exploration described in WB14 as well as the separation of gas masses in Figure 5).

Figure 5 compares the WB14 model grid to our measured 12CO and 13CO line luminosities for sources detected in at least one of these lines. Table 2 provides our $M_{\rm gas}$ constraints derived from the WB14 model grid. To estimate $M_{\rm gas}$ for the three sources detected in both 12CO and 13CO (540, 1274, 1152), we calculate the mean (in log space) of the WB14 model grid points consistent with our measured line luminosities and their associated errors; these values span 2–7 $M_{\rm Jup}$. We also set upper ($M_{\rm gas,max}$) and lower ($M_{\rm gas,min}$) limits based on the maximum and minimum WB14 model grid points consistent with the data, respectively. For the three sources with 12CO detections and 13CO upper limits (1153, 818, 1075) we provide only $M_{\rm gas,max}$, since the lower bound is constrained by the limits of the WB14 model grid, and at these very low masses photodissociation becomes important and self-shielding needs to be taken into account. For the 86 sources undetected in both lines, we give only $M_{\rm gas,max}$, which is set by the maximum WB14 model grid point consistent with the upper limits on both lines.

Alternatively, one can estimate $M_{\rm gas}$ by employing the grid of physical-chemical models presented in Miotello et al. (2016a) (hereafter M16). Using the code DALI (Bruderer et al., 2012), M16 investigated a range of realistic disk and stellar parameters, calculating self-consistently the thermal and chemical structure of the disk up to an age of 1 Myr. As shown by Miotello et al. (2016b) for Lupus disks, the simulated line fluxes can be used to derive $M_{\rm gas}$ from individual 13CO or C18O isotopologue line observations, assuming an abundance of volatile carbon. This is because the medians of the simulated line luminosities can be fit by simple functions of $M_{\rm gas}$ (see Equation 2 in M16). Namely, for low-mass disks the dependence of the line luminosities on $M_{\rm gas}$ is linear, while the trend becomes logarithmic for more massive disks due to the line emission becoming optically thick. Thus, the observed 13CO or C18O line luminosities directly trace $M_{\rm gas}$ only if they fall in the linear dependence regime found in M16.

Applying the methodology of M16 to our $\sigma$ Orionis sample, we find that the 13CO upper limits for the three sources detected only in 12CO (1153, 818, 1075) are in the linear dependence regime and thus provide upper limits on $M_{\rm gas}$. Although the three brightest disks (540, 1274, 1152) have 13CO line luminosities that fall in the logarithmic dependence regime, their C18O upper limits can still provide constraints on $M_{\rm gas}$. The derived $M_{\rm gas}$ upper limits and the CO isotopologue lines used for the calculations are reported in Table 3. These results place stronger constraints on $M_{\rm gas}$ than the WB14 model grid for all six sources detected in CO.

There are several caveats to our derived gas masses. Namely, these gas masses depend inversely on the assumed molecular CO abundance of ${\rm[CO]/[H_{2}]}=10^{-4}$ in the case of WB14, and the assumed volatile carbon abundance of ${\rm[C]/[H]}=1.35\times 10^{-4}$ and chemical age of 1 Myr in the case of M16. Additionally, both WB14 and M16 assume an isotopologue ratio of ${\rm[CO]/[^{13}CO]}=70$. The WB14 ${\rm[CO]/[H_{2}]}$ abundance is consistent with those measured in molecular clouds (Frerking et al., 1982; Lacy et al., 1994; Ripple et al., 2013; Shimajiri et al., 2014) as well as with a direct measurement in a disk (France et al., 2014). However, the strong HD (Bergin et al., 2013) but weak C18O emission toward the TW Hydra disk has been interpreted as resulting from significant carbon depletion of up to two orders of magnitude in this system (Favre et al., 2013; Kama et al., 2016; Schwarz et al., 2016).

If carbon depletion (rather than gas depletion) is the true cause of weak CO emission, the responsible physical mechanisms are not yet established. One hypothesis is that gas-phase reactions initiated by X-ray and cosmic ray ionization of He produce He+ atoms that react with gaseous CO to gradually extract carbon, which is then processed into more complex molecules that freeze onto cold dust grains at higher temperatures than CO (Aikawa et al., 1997; Bruderer et al., 2012; Favre et al., 2013; Bergin et al., 2014; Kama et al., 2016; Yu et al., 2016). Alternatively, CO can be turned into more complex organics such as CH3OH, or into CO2 and CH4, via ice chemistry reactions (e.g., see Figure 3c in Eistrup et al., 2016); these reactions have typical timescales of a few Myr (depending on the ionization rate) and thus could be more significant in older systems. Finally, volatile carbon may be locked up in large icy bodies in the disk midplane (Bergin et al., 2010; Ros & Johansen, 2013; Guidi et al., 2016). These large pebbles cannot diffuse upward and thus would no longer participate in gas-phase chemistry (Du et al., 2015; Kama et al., 2016). Such a process would “dry out” the CO from the warm molecular layer, analogous to what is proposed to explain the under-abundance of gas-phase water in disk atmospheres (Bergin et al., 2010; Hogerheijde et al., 2011). If any of these mechanisms significantly depletes carbon in the disk, our derived gas masses would be underestimated.

V.3. Stacking Analysis

We perform a stacking analysis to constrain the average dust and gas masses of the individually undetected sources in our sample. To stack the images, we average them in the image plane after centering them on their expected source locations. We then search for emission using the aperture photometry method described in Section IV.2.

We first stack the 55 sources undetected in the continuum (“Stack #1”), but do not find a significant mean signal in the continuum or any of the CO lines. The measured continuum mean signal is 0.05$\pm$0.03 mJy (1.7$\sigma$). We confirm this non-detection by calculating the mean and standard error on the mean from the continuum fluxes reported in Table 1, which similarly gives 0.03$\pm$0.02 mJy (1.5$\sigma$). This provides a 3$\sigma$ upper limit on the average dust mass of individually undetected continuum sources of 0.4 $M_{\oplus}$, which is 5$\times$ lower than the smallest dust mass among the continuum detected sources in $\sigma$ Orionis (see Figure 4). This striking difference in the dust masses of detected and undetected continuum sources was also seen in an ALMA survey of Lupus disks (see Figure 3 in Ansdell et al., 2016) and further supports theoretical models that predict protoplanetary disks dispersing rapidly once disk clearing begins (e.g., see review in Alexander et al., 2014).

We also stack the 31 sources that are detected in the continuum but undetected in 12CO (“Stack #2”), finding a mean continuum signal of $2.29\pm 0.09$ mJy as well as a significant mean 12CO signal of $36\pm 8$ mJy km s*-1* (4.5$\sigma$). No emission is detected in the 13CO or C18O lines with 3$\sigma$ upper limits of 14 and 11 mJy km s*-1*, respectively. The continuum flux corresponds to a dust mass of $\sim$10 $M_{\oplus}$, while the 12CO detection and 13CO upper limit correspond to a gas mass $<1.0~{}M_{\rm Jup}$ using the WB14 model grid (the 13CO and C18O upper limits also correspond to a gas mass $<1.0~{}M_{\rm Jup}$). This gives an average gas-to-dust ratio of $<30$ for sources detected in the continuum but not in CO, assuming standard CO abundances.

V.4. Transition Disks

Transition disks (TDs) are protoplanetary disks with large inner cavities in their dust distributions (see review in Espaillat et al., 2014). TDs can be identified by resolved (sub-)mm images, or by the mid-IR deficits in their SEDs, which indicate a lack of warm micron-sized dust grains close to the central star. Eight $\sigma$ Orionis members (1268, 1267, 908, 897, 818, 540, 411, 299) have been identified as TDs based on their SEDs (Hernández et al., 2007; Maucó et al., 2016) and all of these sources are included in our ALMA sample.

Six of these TDs (1267, 897, 818, 540, 411, 299) were detected by our ALMA continuum observations, and two of the continuum-detected sources (818, 540) were also detected in CO. Our ALMA observations did not resolve any dust cavities (see Figure 2), although this is not particularly surprising given the large beam size of our observations ($\sim$120$\times$95 AU at 385 pc; Section III).

Thus 16% (6/37) of our ALMA continuum detections are TDs (see Figure 4), which is interestingly similar to the 19% fraction found in Lupus (Ansdell et al., 2016) and consistent with previous findings that TDs tend to be among the brightest disks in a given star-forming region (Andrews et al., 2011; Ansdell et al., 2016). Moreover, the TD fraction across the entire protoplanetary disk population in $\sigma$ Orionis is 8% (8/92), which is consistent with TD fractions found in other star-forming regions (see Figure 11 in Espaillat et al., 2014).

VI. DISCUSSION

VI.1. External UV Photoevaporation from $\sigma$ Ori

OB associations host at their centers very massive stars with O and B type spectral classifications, surrounded by populations of several hundred to thousands of lower-mass stars. The massive OB stars emit large numbers of extreme-ultraviolet (EUV; $h\nu>13.6$ eV) and far-ultraviolet (FUV; 6 eV$~{}<h\nu<~{}$13.6 eV) photons, which can heat and photoevaporate circumstellar disks around the nearby low-mass stars. The resulting mass-loss rate for a given disk depends on the relative strengths of the incident FUV and EUV fluxes, and thus should depend on the distance of the source from the central OB stars (e.g., Johnstone et al., 1998; Störzer & Hollenbach, 1999).

Mann et al. (2014) used ALMA to search for observational evidence of these external photoevaporation effects on proplyds in the young ($\sim$1–2 Myr) Orion Nebula Cluster (ONC). They found a clear drop in disk mass at small projected separations from the central O6 star, $\theta^{1}$ Ori C. Namely, there was a lack of massive ($M_{\rm dust}\gtrsim 9~{}M_{\oplus}$) disks within $\sim$0.03 pc of $\theta^{1}$ Ori C, where models predict that EUV emission dominates the radiation field (Johnstone et al., 1998). At larger separations of $\sim$0.03–0.3 pc, where less energetic FUV emission is expected to dominate (Adams et al., 2004), they found a range of disk masses representative of typical low-mass star-forming regions, indicating that the lower mass-loss rates in FUV-dominated regions can preserve disk masses for up to a couple Myr. Mann et al. (2014) concluded that planet formation is likely inhibited for disks in the inner-most EUV-dominated regions of OB associations due to higher mass-loss rates, while disks in the FUV-dominated regions and beyond are relatively unaffected with planet formation proceeding as in isolated disks.

Mann et al. (2015) looked for similar effects in the very young NGC 2024 region, which at $\sim$0.5 Myr old (e.g., Levine et al., 2006) hosts the massive star IRS 2b of O8–B2 spectral type (Bik et al., 2003) as well as several hundred YSOs still heavily embedded in molecular cloud material. Mann et al. (2015) could not identify a distance-dependent disk mass distribution in NGC 2024, and instead found several massive ($M_{\rm dust}\gtrsim 17~{}M_{\oplus}$) disks located $<0.01$ pc from IRS 2b. They argued that this could be an evolutionary effect: the extremely young age of NGC 2024 simply means that processes like external photoevaporation have not yet had time to significantly reduce disk masses. Alternatively, they suggested this could be an environmental outcome: the significant cloud material in NGC 2024 may efficiently absorb the high-energy photons from IRS 2b, or the later spectral type of the star (compared to $\theta^{1}$ Ori C in the ONC) means that it does not produce sufficient photoionizing radiation.

Here we search for evidence of external photoevaporation in $\sigma$ Orionis, an OB association whose central trapezium system, $\sigma$ Ori, contains a massive O9 star. $\sigma$ Orionis is an interesting target for studying external photoevaporation, as its lack of cloud material and older age may both enhance the observable effects of external photoevaporation. Figure 6 (upper right panel) plots disk dust mass as a function of projected separation from $\sigma$ Ori ($\alpha$ = 05:38:44.779, $\delta=-$02:36:00.11). Similar to the ONC, we find a lack of massive ($M_{\rm dust}\gtrsim 3~{}M_{\oplus}$) disks close to the central OB system; however, the drop in occurrence is seen at $\sim$0.5 pc—a much larger projected distance compared to the $\sim$0.03 pc limit found for the ONC (Mann et al., 2014). Moreover, beyond $\sim$0.5 pc, we see a smooth distance-dependent dust mass distribution that extends out to several parsecs. We note that the smaller dust masses and larger projected distances found in $\sigma$ Orionis (compared to those found in the ONC) are both influenced by the older age of the region. Namely, dust mass distributions are known to decline with cluster age (Section VI.3) and typical intra-cluster velocity dispersions of several km s*-1* can result in cluster expansions of several parsecs by the age of $\sigma$ Orionis.

One concern is that, due to the $M_{\rm dust}$–$M_{\ast}$ relation (Section VI.2), mass segregation in clusters could produce these observed trends, if the least massive stars are preferentially located closer to the cluster centers. Mann et al. (2014) could not test this in the ONC, as the nature of proplyds complicates any estimates of stellar mass. Because we can estimate stellar masses in $\sigma$ Orionis (Section II), we also show in Figure 6 (lower right panel) the ratio of disk dust mass to stellar mass ($M_{\rm dust}/M_{\ast}$) as a function of projected separation, confirming that the distance-dependent trend still holds even when accounting for stellar mass differences. Moreover, Figure 6 shows our ALMA continuum detection fraction, illustrating a relatively constant detection rate of $\sim$30% out to $\sim$2 pc, after which the detection fraction more than doubles to $\sim$70%. Hernández et al. (2007) did not find a similar change in detection fraction with their Spitzer survey of $\sigma$ Orionis disks (see their Figure 16), however this may be because external photoevaporation does not remove the inner (i.e., more gravitationally bound) disk regions probed by Spitzer.

Interestingly, we also find that the CO detections in our sample (blue circles in Figure 6) only exist in the outer regions of the cluster. This is qualitatively consistent with the picture of external photoevaporation: for typical disks, the gas is generally more extended than the dust, and therefore less tightly bound to the star, making the gas more susceptible to external photoevaporation. However, our gas sample is small and Mann et al. (2014) were unable to reliably detect gas in ONC disks due to cloud confusion, making it important to confirm our finding with surveys of other OB associations. If external photoevaporation does have a more significant effect on gas relative to dust, this would impact the types of planets that can form in OB associations.

Some evidence for external photoevaporation has been previously found for $\sigma$ Orionis disks. Rigliaco et al. (2009) detected strong optical forbidden emission lines from SO 587, which they interpreted as an externally driven photoevaporative flow due to the very low stellar mass accretion rate for this source, the profile shapes and luminosities of the forbidden emission lines, and the small projected separation ($\sim$0.3 pc) of the disk from $\sigma$ Ori (we did not detect this source with our ALMA observations). Additionally, Maucó et al. (2016) fit irradiated accretion disk models to the SEDs of 18 sources in $\sigma$ Orionis to show decreased disk masses and sizes when compared to those in the younger ONC. They interpreted this as evidence for external photoevaporation, however their results were uncertain due to various model assumptions (e.g., constant $\alpha=0.01$) as well as the comparison of disk properties derived from disparate methods (e.g., they compared $\sigma$ Orionis disk radii derived from SED modeling, which probes the dust disk, to ONC disk radii derived from Hubble imagery, which probes the gas disk). Our ALMA observations therefore provide the clearest evidence to date that external photoevaporation is affecting disk masses throughout the $\sigma$ Orionis region.

Our findings also indicate that FUV (not just EUV) emission from OB stars is an important driver of external photoevaporation. Assuming a typical O9V FUV luminosity of $\log{(L_{\rm FUV}/L_{\odot})}=4.5$ for $\sigma$ Ori, the geometrically diluted FUV flux within the region can be expressed as $\sim 8000(d/{\rm pc})^{-2}\,G_{0}$, where $d$ represents the distance from the photoevaporative source in parsecs and $G_{0}=1.6\times 10^{-3}$ erg cm*-2* s*-1* (Habing, 1968). We note that although $\sigma$ Ori is a triple system, the FUV flux is usually dominated by the most massive star in the cluster (Fatuzzo & Adams, 2008; Holden et al., 2011). In this simple calculation, we also do not consider any extinction due to intra-cluster dust, which is observed to be at low densities in $\sigma$ Orionis (Walter et al., 2008), unlike in the ONC and NGC 2024. Figure 6 shows that external photoevaporation is affecting disk masses out to at least $\sim$2 pc, which when combined with the above equation corresponds to FUV fluxes $\gtrsim 2000\,G_{0}$.

This supports recent observations that suggest even moderate FUV fluxes can drive significant mass loss. Kim et al. (2016) observed 7 proplyds near a B star in NGC 1997, finding high mass-loss rates for an FUV flux of only \sim$$3000\,G_{0}. Haworth et al. (2017) also showed that the outer disk of IM Lup may be undergoing photoevaporation from an FUV flux of just \sim$$4\,G_{0}, where the high mass-loss rate can be explained by the large size of the disk (Cleeves et al., 2016), which causes gas in the outer regions to be only weakly gravitationally bound to the central star. Futhermore, Guarcello et al. (2016) found that disk frequency (as probed by near-IR excess) declines with smaller projected separation from the OB stars in Cygnus OB2 for FUV fluxes $\gtrsim 1000\,G_{0}$, similar to what we calculated for $\sigma$ Orionis.

Together, these observations support recent theoretical findings by Facchini et al. (2016) and Haworth et al. (2016), who predicted high mass-loss rates due to external photoevaporation for FUV fluxes $<3000\,G_{0}$. Moreover, these slow photoevaporative winds should be much more effective at removing gas and small ($\lesssim 1~{}\mu$m) dust particles compared to larger ($\gtrsim 1$ mm) solids (Facchini et al., 2016), which may help to explain our lack of gas detections at projected distances $\lesssim 1.5$ pc from $\sigma$ Ori.

Finally, we showed in Section V.3 that the average dust mass of the undetected sources in $\sigma$ Orionis is at least $\sim$5$\times$ lower than the smallest dust mass among the continuum detections, implying that disk dispersal occurs on short timescales. The rapid dispersal of disks impinged by intermediate external FUV fluxes has been predicted by Clarke (2007) and later by Anderson et al. (2013). Their models combined estimated mass-loss rates from external FUV photoevaporation with viscous disk evolution to show that disks should be dispersed from the outside in on timescales much shorter than the expected disk lifetime. The typical lifetime of a viscous disk impinged with a $\sim$3000$\,G_{0}$ FUV flux was predicted to be roughly a few Myr, in agreement with our observations.

In summary, our observations indicate that external photoevaporation due to FUV emission from OB stars is significantly affecting disk evolution throughout the $\sigma$ Orionis cluster, although other disk evolution mechanisms are also clearly at play (Section VI.2 & VI.3). This additional depletion of dust and gas for disks in OB clusters should have implications for planet formation, and detailed theoretical studies will be needed to quantify the impacts for different planet types and to identify any relations to trends seen in the exoplanet population.

VI.2. $M_{\rm dust}$–$M_{\ast}$ Relation

A clear correlation seen across protoplanetary disk populations is the positive relation between $M_{\rm dust}$ and $M_{\ast}$ (e.g., Figure 7). Evidence for the $M_{\rm dust}$–$M_{\ast}$ relation was first identified in pre-ALMA surveys of the Taurus star-forming region (Natta et al., 2000; Andrews et al., 2013) then later recovered with ALMA for other star-forming regions including the similarly aged Lupus clouds (Ansdell et al., 2016) and Chamaeleon I region (Pascucci et al., 2016) as well as the more evolved Upper Sco association (Barenfeld et al., 2016).

The $M_{\rm dust}$–$M_{\ast}$ relation is pertinent to understanding planet formation because it tells us how disks are related to the properties of their host stars, which can then be tied to similar trends seen in the exoplanet population. For example, as previously noted by Andrews et al. (2013), the $M_{\rm dust}$–$M_{\ast}$ relation may fundamentally explain the correlation between giant planet frequency and host star mass (Endl et al., 2006; Johnson et al., 2007; Bowler et al., 2010; Bonfils et al., 2013), as the cores of giant planets theoretically form more efficiently both in higher-mass disks (e.g., Thommes et al., 2008; Mordasini et al., 2012) and around higher-mass stars (e.g., Kennedy & Kenyon, 2008). Moreover, tracking the evolution of the $M_{\rm dust}$–$M_{\ast}$ relation with age can tell us whether disk evolution proceeds differently around low-mass stars compared to high-mass stars, which in turn can help to constrain the relative importance of different disk processes during planet formation (e.g., Pascucci et al., 2016).

However, parameterizing the $M_{\rm dust}$–$M_{\ast}$ relation is complicated by three main factors: measurement uncertainties on both variables, intrinsic scatter in the data, and upper limits. The procedure most often utilized in the disk survey literature is the Bayesian linear regression method of Kelly (2007), as it is capable of accounting for these three key factors simultaneously, unlike other linear regression methods (see Pascucci et al. 2016 for a detailed discussion). For a given dataset, the Kelly (2007) procedure fits a slope ($\beta$), intercept ($\alpha$), and intrinsic dispersion ($\delta$) with associated uncertainties on each parameter.

Using this method, Ansdell et al. (2016) showed that the fitted slopes to the protoplanetary disk populations in the young ($\sim$1–3 Myr) Taurus and Lupus regions were consistent with each other, but both shallower than that of the older ($\sim$5–10 Myr) Upper Sco association. Pascucci et al. (2016) then showed that Chamaeleon I had a slope consistent with the similarly aged Taurus and Lupus regions, further supporting a steepening of the $M_{\rm dust}$–$M_{\ast}$ relation with age. They also compared their results to theoretical models of grain growth, drift, and fragmentation to show that a steepening of the $M_{\rm dust}$–$M_{\ast}$ relation with age is consistent with the outer disk being in the fragmentation-limited regime. In this regime, grain sizes in the outer disk are limited by fragmenting collisions. When fragmentation sets the largest grain size, inward radial drift of dust occurs more rapidly around lower-mass stars, making their (sub-)mm continuum emission weaker and more compact with age compared to higher-mass stars.

Here we derive the $M_{\rm dust}$–$M_{\ast}$ relation for $\sigma$ Orionis disks, again using the Bayesian linear regression method of Kelly (2007). We only consider sources in our ALMA sample with $M_{\ast}\geq 0.1~{}M_{\odot}$, so that we can compare our results to the relations derived for other star-forming regions (see below; this only removes 5 sources from our sample and does not affect the fit results). Using the $M_{\ast}$ and $M_{\rm dust}$ values in Table 1, we derive a linear fit with $\alpha=1.0\pm 0.2$, $\beta=2.0\pm 0.4$, and $\delta=0.6\pm 0.1$ dex, as shown in Figure 7. To help illustrate that the fit is properly accounting for the numerous non-detections, we also show 3$\sigma$ upper limits from stacks of the non-detections in several stellar mass bins.

To compare our $M_{\rm dust}$–$M_{\ast}$ relation derived for $\sigma$ Orionis to those found for other star-forming regions in a consistent manner, we follow the procedure described in Ansdell et al. (2016). Namely, we calculate $M_{\rm dust}$ uniformly across each region by inputting the (sub-)mm continuum fluxes (or $3\sigma$ upper limits) from the literature into Equation 1, along with the cluster distances and observation wavelengths of the surveys. We assume $T_{\rm dust}=20$ K for all disks and adopt distances of 140 pc for Taurus (Kenyon et al., 2008), 150 pc or 200 pc for Lupus (Comerón, 2008), 160 pc for Chamaeleon I (Luhman, 2008), and 145 pc for Upper Sco (de Zeeuw et al., 1999). For Upper Sco, we only include the “full,” “evolved,” and “transitional” disks from the sample of Barenfeld et al. (2016), as these represent the “primordial” disks that do not yet show signs of disk clearing. Only sources with $M_{\ast}\geq 0.1~{}M_{\odot}$ were considered in order to exclude brown dwarfs, while also maintaining a common stellar mass limit among the surveys. Stellar masses taken from the literature were derived using the Siess et al. (2000) evolutionary tracks for all regions except Chamaeleon I, which used the Baraffe et al. (2015) models; the stellar masses derived using these two grids are generally consistent, thus any effects should be negligible.

The fitted linear regression parameters for each region are given in Table 4 and plotted in Figure 7. Our results further support the steepening of the $M_{\rm dust}$–$M_{\ast}$ relation with age, although the errors are large. Interestingly, we also find similarly large intrinsic dispersions for all five regions; as previously noted by Pascucci et al. (2016), this seems to be an inherent property of disk populations, reflecting a range of disk conditions (e.g., dust opacities, disk evolutionary states, dust temperatures) rather than the age and/or environment of the region, and may partially account for the diversity seen in the exoplanet population.

We note that our fitted values in Table 4 are mostly consistent with those found by Pascucci et al. (2016) (see their Table 4) despite differences in assumptions of grain opacity (e.g., they use $\beta=0.4$ in Equation 1 while we use $\beta=1.0$) and stellar mass cutoffs (e.g., they include sources with $M_{\ast}<0.1~{}M_{\odot}$, while we exclude brown dwarfs). Indeed, the main disagreement is the intercept estimate for Lupus, which differs because Pascucci et al. (2016) exclude the 20 sources in Lupus with unknown stellar masses, while we account for them using the MC approach described in Ansdell et al. (2016). The slope for Upper Sco is also noticeably different (although within errors) because Pascucci et al. (2016) only consider “full” and “transitional” disks from Upper Sco, while we also include “evolved” disks following the definition of “primordial” disks in Barenfeld et al. (2016).

Finally, we address three potential caveats to our Bayesian linear regression fit to $\sigma$ Orionis disks. First, Pascucci et al. (2016) found that shallower slopes can result when the sample is dominated by upper limits at low stellar masses. Although roughly two-thirds of our $\sigma$ Orionis sample was undetected in the continuum, the non-detections span a range of stellar masses (see Figure 7); moreover, even if the slope is actually steeper, this would only further distinguish $\sigma$ Orionis from the younger regions. Second, the source with the highest stellar mass in $\sigma$ Orionis is not detected in the continuum, which could potentially skew the fit. We could not find a valid reason for discarding this source (e.g., no evidence for binarity), although its spectral type is uncertain ($\pm$2.5 spectral types; Hernández et al., 2014). Nevertheless, when re-fitting the distribution without this source, we recover the same parameters. Third, we have shown in Section VI.1 that external photoevaporation is reducing disk dust masses in $\sigma$ Orionis. This may serve to steepen the $M_{\rm dust}$–$M_{\ast}$ relation, as external photoevaporation should be more effective at removing gas and small dust grains around lower-mass stars whose orbiting material is less gravitationally bound. Unfortunately this possible effect has not yet been tested with external photoevaporation models. Given these potential issues, a more sensitive (sub-)mm continuum survey of $\sigma$ Orionis, as well as additional observations and detailed modeling, will be needed to confirm the linear regression fit presented here. Furthermore, we note that several other linear regression methods for left-censored datasets are available in the statistical literature and should be tested for this specific astronomical problem.

VI.3. Dust Mass Distributions

Disk dispersal and grain growth should be reflected in a steady decline with age of the bulk dust mass probed by (sub-)mm continuum flux. Ansdell et al. (2016) showed that the average dust mass of disks in a given star-forming region, as measured from (sub-)mm continuum flux, does indeed decline with age. Specifically, they found that Lupus and Taurus have consistent mean dust masses ($15\pm 3$ $M_{\oplus}$ and $15\pm 2$ $M_{\oplus}$, respectively), while Upper Sco has a mean dust mass that is $\sim$5$\times$ lower ($5\pm 3$ $M_{\oplus}$). They calculated dust masses uniformly across each region (as described in Section VI.2) then derived the average dust mass using the Kaplan-Meier estimator in the ASURV package (Lavalley et al., 1992) to properly account for the upper limits using well-established techniques for left-censored datasets. We follow this procedure to now include results from our ALMA survey of $\sigma$ Orionis disks as well as the ALMA survey of Chamaeleon I disks (Pascucci et al., 2016). As shown in Figure 8, we find average disk dust masses of $7\pm 1$ $M_{\oplus}$ for $\sigma$ Orionis and $13\pm 4$ $M_{\oplus}$ for Chamaeleon I, further supporting the decline of this parameter with age.

When comparing the dust mass distributions of two regions, it is important to confirm that they have similar stellar mass distributions due to the $M_{\rm dust}$–$M_{\ast}$ relation discussed in Section VI.2. We therefore employ the two-sample tests in the ASURV package to determine the probabilities that the stellar masses of each region are drawn from the same parent population. We find that the stellar mass distributions are statistically indistinguishable for $\sigma$ Orionis and Lupus ($p=0.45$–0.96), $\sigma$ Orionis and Chamaeleon I ($p=0.06$–0.22), and $\sigma$ Orionis and Upper Sco ($p=0.21$–0.30). However we found statistically distinct stellar mass distributions for $\sigma$ Orionis and Taurus ($p=0.03$–0.04). Therefore, we can directly compare the dust mass distribution of $\sigma$ Orionis to those of Lupus, Chamaeleon I, and Upper Sco in Figure 8. However, caution should be taken when interpreting the comparison of $\sigma$ Orionis to Taurus in Figure 8 due to their potentially different stellar mass distributions. Nevertheless, our main conclusions in the previous paragraph remain the same.

Alternatively, to account for stellar mass differences between regions, one could instead compare the $M_{\rm dust}/M_{\ast}$ distributions (e.g., Barenfeld et al., 2016) or employ an MC approach that aims to normalize the stellar mass selection functions (e.g., Andrews et al., 2013; Williams et al., 2013; Ansdell et al., 2016). However, in this work we do not attempt these more detailed analyses given the larger uncertainties on our $M_{\ast}$ estimates, especially for the sources with photometrically derived stellar masses.

VI.4. Ingredients for Giant Planet Formation

Core accretion theory predicts that giant planets form when solid cores of a minimum critical mass assemble in the disk, enabling runaway accretion of the surrounding gaseous material (Pollack et al., 1996; Ida & Lin, 2004). The accretion of a gaseous envelope is expected to occur rapidly, where $\sim$10 $M_{\oplus}$ cores reach masses of $\sim$1 $M_{\rm Jup}$ within $\sim$0.1 Myr. Within the framework of this model, we can constrain the occurrence of giant planet formation by observing how quickly the dust and gas content in typical protoplanetary disks depletes to levels below what are thought to be needed to form a gas giant.

In particular, we can look at the fraction of protoplanetary disks in a region with dust masses above the \sim$$10~{}M_{\oplus} limit needed to form a giant planet core. For regions at $\sim$1–3 Myr, we see roughly a quarter of protoplanetary disks above this threshold (30% in Taurus, 26% in Lupus, and 23% in Chamaeleon I). At $\sim$3-5 Myr, we see this fraction cut in half (13% in $\sigma$ Orionis) and then halved again at $\sim$5–10 Myr (5% in Upper Sco). Although these are only rough estimates due to different survey completenesses, they seem to clearly reflect a sharp decline in the capacity of disks to form giant planets with age. Even in the youngest regions, the majority of disks appear to lack sufficient dust to form the solid cores needed to build giant planets, implying that giant planet formation is either rare or well on its way after just a few Myr.

Additionally, stacking (sub-)mm continuum non-detections allows us to put limits on the average amount of dust in the lowest-mass disks (e.g., Section V.3). In Lupus, Ansdell et al. (2016) found that the undetected disks had extremely low average dust masses of $\lesssim 6$ Lunar masses (0.03 $M_{\oplus}$), comparable to debris disk levels (Wyatt, 2008). Although the further distance of $\sigma$ Orionis results in looser constraints, we still find that undetected disks have $\lesssim 4$ Martian masses (0.4 $M_{\oplus}$) of dust on average (Section V.3). These findings support theoretical predictions that viscous disks evolve rapidly into debris disks once stellar accretion ceases and photoevaporation from the central star dominates, clearing the dust from the inside out and leaving behind larger solids such as pebbles and planetesimals (e.g., Clarke et al., 2001; Alexander et al., 2006). Previous studies of weak-lined T Tauri stars have provided observational evidence for rapid disk clearing (e.g., Cieza et al., 2013; Hardy et al., 2015), and Panić et al. (2013) used a compilation of sub-mm fluxes from the literature to show that debris disks are substantially less massive than disks around younger pre-main sequence stars. However, our larger and more homogeneous samples of Lupus and $\sigma$ Orionis disks confirm that rapid disk clearing is a uniform occurrence, even among young protoplanetary disk populations.

Another ingredient for giant planet formation is of course the gas. Although bulk gas masses are notoriously difficult to measure (Section V.2), the fact that we found only three disks in $\sigma$ Orionis that exhibit both 12CO and 13CO emission is telling; indeed, their line fluxes correspond to gas masses of just $\sim$2–7 $M_{\rm Jup}$ using the methodology described in Section V.2. For the remaining disks in $\sigma$ Orionis, we find upper limits on their individual gas masses of just $\sim$3 $M_{\rm Jup}$. Moreover, the average gas mass of disks detected in the continuum but undetected in 12CO is $<1.0~{}M_{\rm Jup}$ (Section V.3). These low gas masses again suggest that giant planet formation is either rare or nearly complete by the $\sim$3–5 Myr age of $\sigma$ Orionis (or that carbon is being significantly depleted in protoplanetary disks; see Section V.2).

VI.5. Relation between $M_{\rm disk}$ and $\dot{M}_{\rm acc}$

Viscously evolving protoplanetary disks should exhibit a direct relation between their total disk mass ($M_{\rm disk}$) and stellar mass accretion rate ($\dot{M}_{\rm acc}$), as shown in Hartmann et al. (1998). This theoretical prediction was only recently confirmed observationally by Manara et al. (2016), who combined the ALMA (Ansdell et al., 2016) and VLT/X-Shooter (Alcalá et al., 2014, 2016) surveys of protoplanetary disks in the Lupus clouds to reveal a linear relation between $\dot{M}_{\rm acc}$ and $M_{\rm disk}$. Manara et al. (2016) also showed that, when assuming an ISM gas-to-dust ratio of $\sim$100, the $M_{\rm disk}$/$\dot{M}_{\rm acc}$ ratio is consistent with the age of the region, as expected for viscously evolving disks (Rosotti et al., 2017). These observational findings suggest that mass accretion onto the stellar surface is indeed related to the properties of the outer disk.

To search for a similar correlation in $\sigma$ Orionis, we combine our ALMA observations with the $U$–band survey of Rigliaco et al. (2011). The latter used their photometric data to estimate $\dot{M}_{\rm acc}$, making these accretion rates more uncertain than those obtained with spectroscopy for Lupus disks. As shown in Figure 9, only 20 sources (the most massive disks in $\sigma$ Orionis) have both $M_{\rm dust}$ estimates and constraints on $\dot{M}_{\rm acc}$ (the latter of which are mostly upper or lower limits). The sparseness of this sub-sample, combined with the smaller $M_{\rm dust}$ range when compared to the Lupus sample, limits our ability to fit a relation similar to that in Manara et al. (2016). However, we note that, as for Lupus disks, the $M_{\rm disk}$/$\dot{M}_{\rm acc}$ ratios are consistent with the age of $\sigma$ Orionis for an ISM gas-to-dust ratio of $\sim$100.

The disks undetected with ALMA (i.e., the least massive disks in $\sigma$ Orionis) span over 2 dex in $\dot{M}_{\rm acc}$, and thus a correspondingly large range of $M_{\rm disk}$/$\dot{M}_{\rm acc}$ ratios, as shown in Figure 9. The undetected disks with low $\dot{M}_{\rm acc}$ values or upper limits are consistent with expectations from viscous evolution. However, the undetected disks with moderate-to-high $\dot{M}_{\rm acc}$ values ($\gtrsim 2\times 10^{-10}M_{\odot}$ yr*-1*) are unexpected, as these should have lifetimes shorter than the age of the region. The sources $\lesssim 0.5$ pc from the central OB star (outlined in red in Figure 9) are most readily explained by external photoevaporation (see Section VI.1), which would serve to reduce $M_{\rm disk}$ by removing mass from the outer disk, thereby increasing the $M_{\rm disk}/\dot{M}_{\rm acc}$ ratio (Rosotti et al., 2017). Although this accounts for only a handful of objects, there is evidence that external photoevaporation is occurring throughout the region (Section VI.1), thus may apply to more of the undetected disks. Additionally, the $\dot{M}_{\rm acc}$ values estimated from $U$-band photometry are uncertain and need to be confirmed with spectroscopy. Nevertheless, there are a sufficient number of sources with low disk masses and significant accretion rates to warrant further investigation. These objects might have strongly variable accretion rates, or alternatively the accreting gas may come from the evaporation of ice-coated dust grains.

VII. Summary

We have used ALMA to conduct a high-sensitivity mm survey of protoplanetary disks in the $\sigma$ Orionis cluster. This region is particularly interesting for studying disk evolution as its intermediate age (\sim$$3–5 Myr) is comparable to the median disk lifetime, and therefore corresponds to a potentially important phase of disk evolution and planet formation.

We used ALMA to survey the dust and gas in 92 protoplanetary disks around $\sigma$ Orionis members with $M_{\ast}\gtrsim 0.1~{}M_{\odot}$. Our observations cover the 1.33 mm continuum as well as the 12CO, 13CO, and C18O $J=2$–1 lines. Out of the 92 sources, we detected only 37 in the mm continuum and six in 12CO, three in 13CO, and none in C18O. 2. 2.

The continuum emission constrained dust masses to $\sim$2 $M_{\oplus}$, while the CO line emission constrained gas masses to $\sim$3 $M_{\rm Jup}$. Only 11 disks had $M_{\rm dust}\gtrsim 10~{}M_{\oplus}$, indicating that after a few Myr of evolution the vast majority of disks lack sufficient dust to form giant planet cores. The low gas masses also indicate that giant planet formation must be rapid or rare, but may also reflect significant carbon depletion in protoplanetary disks. Moreover, stacking the individually undetected continuum sources limited their average dust mass to $\sim$5$\times$ lower than that of the faintest detected disk, supporting theoretical models that predict disks dissipating rapidly once accretion stops and photoevaporation dominates. 3. 3.

We found that external photoevaporation from the central OB stars is influencing disk evolution throughout the region. Namely, disk dust masses clearly decline with decreasing projected separation from the photoionizing source, and the handful of cluster members with detected CO emission exist only at projected separations $>1.5$ pc. This indicates that even moderate external FUV fluxes can result in significant mass-loss rates, and future theoretical studies will be needed to quantify the implications for planet formation in OB clusters. 4. 4.

Comparing the protoplanetary disk population in $\sigma$ Orionis to those of other star-forming regions provided continuing support of the steady decline in average disk dust mass and steepening of the $M_{\rm dust}$–$M_{\ast}$ relation with age. Quantifying these evolutionary trends can help to determine the relative importance of different disk processes during key eras of planet formation. However, for $\sigma$ Orionis, these trends may also be influenced by the effects of external photoevaporation from the central OB stars. 5. 5.

Collectively, our findings indicate that giant planet formation is inherently rare and/or well underway by a few Myr of age. However, due to the abundance of upper limits in our ALMA sample, and the need for better constraints on stellar properties, a higher-sensitivity (sub-)mm survey as well as a complete spectroscopic survey of the members of $\sigma$ Orionis analyzed in this work should be conducted to confirm these results.

MCA and JPW were supported by NSF and NASA grants AST-1208911 and NNX15AC92G, respectively. MCA acknowledges student observing support from NRAO. CFM acknowledges an ESA Research Fellowship. NM is supported in part by the Beatrice W. Parrent Fellowship in Astronomy at the University of Hawaii. Leiden is supported by the European Union A-ERC grant 291141 CHEMPLAN, by the Netherlands Research School for Astronomy (NOVA), and by grant 614.001.352 from the Netherlands Organization for Scientific Research (NWO). This paper makes use of the following ALMA data: ADS/JAO.ALMA2015.1.00089.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Re- public of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.