A HERSCHEL VIEW OF PROTOPLANETARY DISKS IN THE σ ORI CLUSTER

As shown in Figure 2, 142 TTSs in σ Ori (H07a) fall in the PACS FOV but only 23% of them are detected. Here we compare the detections with the entire TTS population in the field, aiming to understand what makes them different.

The left panel of Figure 4 shows the PACS 70 μm detections in the [3.6]–[4.5] versus [5.8]–[8.0] diagram, as well as all the σ Ori sources in the PACS field. Sources around [0, 0] are referred to as diskless stars while PACS detections are identified following the disk classification of H07a. Note that most of the PACS-detected sources are inside the rectangle that encompasses the optically thick disk region (Hartmann et al. 2005; D'Alessio et al. 2006) and are considered as CTTSs. This is consistent with the shape of their SEDs and their similarities with the Taurus median (Figure 3). Excluding the TDs/PTDs and the two evolved disks in our sample, which have a decrease of IR excess in the IRAC bands and therefore will have lower colors than those predicted by optically thick disks, four sources fall outside the loci of CTTSs: SO823, SO462, SO774, and SO927. The source SO823 has a low [5.8]–[8.0] color. The others (SO462, SO774, and SO927) present higher [5.8]–[8.0] colors. The stars SO462 and SO774 seem to have a slight decay in the first three IRAC bands (3.6, 4.5, 5.8 μm), which causes an apparent redder [5.8]–[8.0] color (note their high 8 μm fluxes resembling those of TDs/PTDs in our sample). The source SO927 has a peculiar SED with very large near-IR and mid-IR excesses.

Zoom In Zoom Out Reset image size

The right panel in Figure 4 displays the K–[5.8] versus [8.0]–[24] SED slopes. The slope of the SED is defined as

$$\begin{eqnarray}&&n=\displaystyle \frac{\mathrm{log}({\lambda }_{1}{F}_{{\lambda }_{1}})-\mathrm{log}({\lambda }_{2}{F}_{{\lambda }_{2}})}{\mathrm{log}({\lambda }_{1})-\mathrm{log}({\lambda }_{2})},\end{eqnarray} \tag{ 1 }$$

where ${\lambda }_{i}$ and ${F}_{{\lambda }_{{i}}}$ are the wavelength and the observed flux at that particular wavelength, respectively. This parameter has been used to classify young stellar objects (Lada 1987, 2006; H07a). We also indicate the lower limit of primordial disks in Taurus (Luhman et al. 2010, dashed line). Photospheric limits are shown by dotted lines. H07a have classified the object SO540 as a class II star based on its slopes between 3.6 and 8.0 μm; however, using the K–[5.8] and the [8]–[24] slopes indicates that it could be a TD candidate. Another interesting point to notice is the fact that one of our two evolved disks, the source SO615, appears just in the limit of primordial disks and has the lowest [8.0]–[24] color of the PACS sample. This may be the result of chromospheric contamination, which is more significant in evolved stars with low mass accretion rates, where magnetic activity from the chromosphere of the star may produce an extra heating of dust, which is reflected at 5.8 μm (Ingleby et al. 2013).

Figure 5 shows the [24] versus K–[24] color–magnitude diagram of σ Ori members, highlighting our PACS disks. Also shown are the histograms of MIPS 24 μm (right) and of K–[24] color (top). As shown in the figure, our detections have the largest excesses at 24 μm, despite their similar near- and mid-IR colors to those of non-detections (Figure 4), with an average magnitude of 5.65 that corresponds to a flux of 67.27 mJy. The lowest 24 μm magnitude detected corresponds to the source SO514 with a value of 8.3 mag. Non-detections, on the other hand, present an average magnitude of 8.16. The source SO566 is a variable star with variability amplitude in the J band larger than 2.7 mag (H14), as can be seen in its non-contemporaneous SED in Figure 3; thus, its large K–[24] color is not real.

Zoom In Zoom Out Reset image size

Figure 6 shows a histogram of spectral types for σ Ori members (H14). PACS 70 μm detections cover spectral types from K2.0 to M5.0. Objects with spectral types later than M5 are unlikely to be detected because the average flux at 70 μm would be below the detection limit by more than a factor of 2 (Section 2). On the other hand, disk-bearing stars with spectral types earlier than K0 have evolved disks with 24 μm excesses below the medians of Taurus and σ Ori, and therefore 70 μm fluxes below the detection limit, or are located out of the PACS FOV.

Zoom In Zoom Out Reset image size

4.1.1. Statistics of Transitional and Pre-transitional Disks

We estimated stellar and accretion properties of all the σ Ori members reported in H07a and H14 that lie inside the PACS FOV, for which we have the necessary spectra and photometry. This information complements the spectroscopic census of σ Ori sources made by H14. To characterize the stellar properties of the sources we located them in the H–R diagram. For this, we estimated the luminosity of our sample using J and V photometry, visual extinctions ${A}_{{\rm{v}}}$, and spectral types from H14. For stars of F0 or later, bolometric correction and effective temperatures were obtained from the standard table for 5–30 Myr old pre-main-sequence (PMS) stars from Pecaut & Mamajek (2013, Table 6). For stars earlier than F0 we used the standard table of main-sequence stars reported by Pecaut & Mamajek (2013, Table 5). Using these luminosities and effective temperatures we estimated stellar radii for σ Ori members. We adopted the Mathis reddening law (Mathis 1990, R = 3.1). We assumed a distance to the σ Ori cluster of 385 pc (Caballero 2008c).

The H–R diagram for σ Ori members inside the FOV of PACS is displayed in Figure 7. Objects detected with PACS are, in general, late-type PMS stars, as discussed in Section 4.1. Using the isochrones of Siess stellar models (Siess et al. 2000), we were able to calculate the age and mass of all our stars in the sample. Stellar properties are listed in Tables 2 and 3 for detected and undetected sources, respectively.

Zoom In Zoom Out Reset image size

Mass accretion rates ($\dot{M}$) were estimated from the ${\rm{H}}\alpha$ line luminosity (Muzerolle et al. 2003, 2005; Natta et al. 2004; Rigliaco et al. 2012; Ingleby et al. 2013). The Hα luminosity was estimated by approximating the flux of the line (${F}_{{\rm{H}}\alpha }$) as EW(Hα) × ${F}_{{\rm{cont}}}$, where ${F}_{{\rm{cont}}}$ and EW(Hα) are the continuum flux around the line and the equivalent width, respectively. We calculated ${F}_{{\rm{cont}}}$ from the ${R}_{{\rm{c}}}$ magnitude of each source reported in H07a corrected for extinction, and the equivalent widths EW(Hα) from low-resolution spectra reported in H14. Once the line fluxes are estimated, the line luminosities are given by

$$\begin{eqnarray}&&{L}_{{\rm{H}}\alpha }=4\pi {d}^{2}{F}_{{\rm{H}}\alpha },\end{eqnarray} \tag{ 2 }$$

where d is the distance to the σ Ori cluster.

In order to obtain the accretion luminosities we used the relation between the Hα luminosity (${L}_{{\rm{H}}\alpha }$) and the accretion luminosity (${L}_{{\rm{acc}}}$) from Ingleby et al. (2013):

$$\begin{eqnarray}&&\mathrm{log}({L}_{\mathrm{acc}})=1.0(\pm 0.2)\mathrm{log}({L}_{{\rm{H}}\alpha })+1.3(\pm 0.7).\end{eqnarray} \tag{ 3 }$$

Finally, ${L}_{{\rm{acc}}}$ can be converted into $\dot{M}$ using

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{{L}_{\mathrm{acc}}{R}_{* }}{{{GM}}_{* }},\end{eqnarray} \tag{ 4 }$$

where ${M}_{* }$ and R_* are the stellar mass and radius respectively. Three sources have no estimation of mass accretion rates since they do not have photometry in the R_c band (SO540, SO566, and SO638).

Figure 8 shows a histogram of the mass accretion rates for the σ Ori cluster. PACS (70 μm) detections are plotted as gray bars and undetected sources as light blue bars. PACS detections exhibit $\dot{M}\leqslant {10}^{-8}$ ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. Note, however, that values of $\dot{M}\lt {10}^{-9}$ ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ are unreliable because of chromospheric contamination, as demonstrated by Ingleby et al. (2013), where an active chromosphere can mask all evidence of an accretion shock excess for the lowest accretors. Additionally, chromospheric activity adds uncertainty to the estimation of the luminosities of the lines associated with accretion for low-mass stars at later evolutionary stages (Manara et al. 2013). As noted in the figure, 23 undetected sources have $\dot{M}\geqslant {10}^{-9}$ ${M}_{\odot }\,{\mathrm{yr}}^{-1}$; this is consistent with the discussion of Figures 4 and 5 where the undetected CTTSs in the PACS FOV are those with the lowest excesses at 24 μm. Most of these sources are the less massive and intrinsically weaker CTTSs in the cluster.

Zoom In Zoom Out Reset image size

We compared our mass accretion rates with those reported by Rigliaco et al. (2011), who estimated $\dot{M}$ from U-band excesses. In general, our mass accretion rates are higher than Rigliaco's. For the PACS sample we found that 10 out of 19 sources are consistent with the values reported by Rigliaco et al. (2011) within a factor of less than 3, while for the undetected sources we found that 10 out of 16 sources agree within the same factor. Rigliaco et al. (2011) neglected extinction, but H07a found extinctions up to 2.94 mag, which can be significant at the U band. This can explain some of the discrepancies between their values and ours. Another source of uncertainty is due to the fact that mass accretion rates determined from the U-band excess are not adequate for sources with $\dot{M}$ <10⁻⁹ ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ because of chromospheric contamination (Ingleby et al. 2013), which is the case for several of our undetected sources. Therefore, we also compared our mass accretion rates with those reported by Rigliaco et al. (2012), who estimated $\dot{M}$ from the Hα luminosities. Even though they also neglected extinction, accretion rates estimated from Hα line luminosities are less affected by extinction than $\dot{M}$ determined from U-band excesses. We found that, for the six sources we have in common, all but one (SO490) agree within a factor of ≲3. Gatti et al. (2008) used near-IR hydrogen recombination lines (Paγ) to measure mass accretion rates of 35 objects in σ Ori. If we compare our mass accretion rates with Gatti's we find that, for the 25 objects we have in common, all but three agree within a factor of ≤3. The main uncertainties are due to the stellar mass adopted in each case, which in turn depends on the luminosity and the extinction of the sources, and on the spectral type.

All the accretion parameters of PACS sources are shown in Table 2. Additionally, Table 3 lists the accretion parameters for σ Ori members in the PACS FOV but without detections.

We follow the methods of D'Alessio et al. (2006) to calculate the structure and emission of accretion disk models. The main input parameters are the stellar properties (M_*, R_*, L_*), the mass accretion rate ($\dot{M}$), the viscosity parameter (α), the cosine of the inclination angle (μ), the disk outer radius (${R}_{{\rm{d}}}$), the maximum grain size at the disk midplane (${a}_{\mathrm{maxb}}$), and the dust settling, assumed to be constant throughout the disk.

The composition of the disk consists of amorphous silicates (pyroxenes) with a mass fraction of 0.004 and carbonates in the form of graphite with a mass fraction of 0.0025 (all these mass fractions are relative to gas). The PACS fluxes did not show a strong dependence on H₂O abundance, so we kept this parameter fixed at 10⁻⁵. The emission from the inner edge or wall of the dust disk is calculated from the stellar properties, the maximum grain size (${a}_{\mathrm{maxw}}$), and the temperature (${T}_{{\rm{wall}}}$), which is assumed to be the sublimation temperature of silicate grains (1400 K) for the diffuse ISM (Draine & Lee 1984; D'Alessio et al. 2006), with a dust composition formed entirely by pyroxenes.

To simulate the settling of dust, D'Alessio et al. (2006) considered the ideal case of two dust populations, which differ in grain size. Both populations follow a size distribution with a function of the form $n(a)\propto {a}^{-p}$, where a is the radius of the grain, whose minimum value is ${a}_{{\rm{\min }}}=0.005$ μm, and p is an exponent equal to 3.5, resembling the size distribution found in the ISM. The first population corresponds to small grains with ${a}_{{\rm{\max }}}=0.25$ μm—characteristic of ISM dust—which is expected in unevolved disks. These grains are mostly located in the upper layers and are assumed to be well mixed with the gas throughout the disk. The second population consists of larger grains, ${a}_{{\rm{\max }}}=1$ mm, and lives in the disk midplane. As the settling process takes place on the disk, dust particles will leave the upper layer and leave behind a depleted population of small grains that has a dust-to-gas mass ratio, ${\zeta }_{{\rm{small}}}$, lower than the initial value, ${\zeta }_{{\rm{std}}}$. On the other hand, the population of big grains will have a dust-to-gas mass ratio, ${\zeta }_{{\rm{big}}}$, larger than the initial value, since small dust particles that settle from upper layers will add mass to this population. Therefore, a decrease in ${\zeta }_{{\rm{small}}}$ automatically implies an increase (of the same proportion) in ${\zeta }_{{\rm{big}}}$. We parameterize settling with the parameter , referred to as the dust-to-gas mass ratio and defined as $\epsilon ={\zeta }_{{\rm{small}}}/{\zeta }_{{\rm{std}}}$, i.e., the mass fraction of the small grains relative to the standard value. Therefore, lower values of represent more settled disks.

Since evolutionary effects, like dust settling and grain growth, are more apparent at longer wavelengths, we made plots of spectral indices using the PACS photometry reported here and IRAC and MIPS fluxes from H07a, in order to determine the overall degree of dust settling in our disks.

The left panel of Figure 9 displays a diagram of n_4.5–24 versus n_24–70 for σ Ori members (solid circles). We compare the observed spectral indices with theoretical indices calculated for a set of models with the following parameters: ${M}_{* }=0.5$ ${M}_{\odot }$, ${R}_{* }=2$ ${R}_{\odot }$, ${T}_{{\rm{eff}}}=4000$ K (all these parameters constitute mean values in the cluster as shown in Tables 2 and 3) and two values of mass accretion rate, $\dot{M}={10}^{-8}\,{M}_{\odot }{\mathrm{yr}}^{-1}$ (light blue symbols) and $\dot{M}={10}^{-9}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (red symbols), where the latter is a common value for σ Ori sources (Figure 8). We fixed the viscosity parameter α to a typical value of 0.01 (Hartmann et al. 1998). As we already said in Section 4.2, we used a distance to the σ Ori cluster of 385 pc. Predicted spectral indices are indicated by open circles for models with ${R}_{{\rm{d}}}=250$ au, and plus signs for models with ${R}_{{\rm{d}}}=10$ au, where the size of the symbol varies according to the degree of dust settling, with larger symbols representing more settled disks (lower values of ). Squares and triangles account for TD/PTD disks and evolved disks, respectively. The scatter observed in each case represents the variety of properties such as maximum grain size at the disk midplane, ${a}_{{\rm{maxb}}}$, the maximum grain size at the disk wall, ${a}_{{\rm{maxw}}}$, disk inclination angle, μ, disk external radius, ${R}_{{\rm{d}}}$, and ice abundances adopted in the models. A similar plot is shown in the right panel of Figure 9 for the spectral index n_24–160 on the x axis. As shown, models with different degrees of settling populate different regions on this diagram, the more settled disks being bluer than the less settled ones. Note as well the dependence of mass accretion rate on the n_4.5–24 slope. Models with $\dot{M}={10}^{-9}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ have n_4.5–24 spectral indices as low as −1.5, while models with $\dot{M}={10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ have indices higher than −1.0. As can be inferred from these diagrams, the majority of our objects fall in regions of  ≤ 0.01, indicating that most of the σ Ori disks have a significant degree of settling. In both cases (left and right panels of Figure 9) the indices n_24–70 and n_24–160 for a set of objects lie outside the region populated by the models. However, in Section 5.3 we modeled almost all of these outliers.

Zoom In Zoom Out Reset image size

Figure 10 displays the n_4.5–24 versus n_24–70 spectral indices of σ Ori sources (black dots) compared to disks in Taurus (1–3 Myr, diamonds) from Howard et al. (2013) and Luhman et al. (2010). Also shown are the medians and quartiles for the n_24–70 (top) and n_4.5–24 (right) indices of each region for comparison. The σ Ori sample has spectral indices similar to the population of Taurus.

Zoom In Zoom Out Reset image size

In order to inquire deeply into this subject as well as to understand the degree of dust settling in the sample, we compared in Figure 11 histograms of observed n_24–70 and n_24–160 indices for σ Ori (black) and Taurus (light blue) disks. The ranges of the spectral indices for the disk models shown in Figure 9 were plotted as bands on top of these figures for values of 1.0, 0.1, 0.01, 0.001, and 0.0001 from top to bottom. The Taurus histograms were made with MIPS 24 μm photometry from Luhman et al. (2010) while the PACS fluxes were taken from Howard et al. (2013), giving an overall average of 21 sources. As shown in the left panel, most of our objects have a considerable degree of dust settling, since the peak of the histogram, labeled as σ Ori, falls in the same range as models with  ≤ 0.01. A similar situation happens for the n_24–160 index, shown in the right panel. In this case, however, there is more overlap between ranges occupied by models with different degrees of settling. These results imply that disks in our sample have experienced some evolution over time. Note, as well, how disks in Taurus have a lower degree of dust settling, since the peaks of the Taurus histograms fall slightly to the right (bigger values of ) compared to a significant number of σ Ori sources with  ≤ 0.01 in the figures. This is consistent with previous studies of dust evolution where older disks are more settled and emit smaller IR excesses than young ones, which causes a decrease in the disk fraction as a function of age for a given stellar group (H07a). It is important to note that the fact that disks on Taurus already exhibit signs of dust settling indicates that this process must happen at an early evolutionary stage.

Zoom In Zoom Out Reset image size

To better constrain the dust content in our PACS disks classified as class II stars and EV disks (the TDs will be studied in a future paper), we modeled the SEDs of individual objects with irradiated accretion disk models, using as input the stellar properties and accretion rates reported in Table 2, assuming a viscosity parameter α = 0.01 and a distance to σ Ori of 385 pc (Section 4.2). For each object we calculated a total of 1800 disk models, varying the outer radius (${R}_{{\rm{d}}}$), the degree of dust settling (), the cosine of the inclination angle (μ), and the maximum grain size at the disk midplane (${a}_{{\rm{maxb}}}$) and at the wall (${a}_{{\rm{maxw}}}$). Table 6 lists all the relevant parameters. We assumed a maximum grain size in the upper layer of the disk of 0.25 μm. We selected as the best fit the model that yielded the minimum value of ${\chi }^{2}$. Upper limits in the observed fluxes were excluded from the calculation of the ${\chi }^{2}$ value.

Table 6.  Model Parameters

Parameter	Value
${R}_{{\rm{d}}}$ (au)a	10, 30, 50, 70, 100, 150, 200, 250, 300, 350
	1, 0.1, 0.01, 0.001, 0.0001
μ	0.3, 0.6, 0.9
${a}_{{\rm{maxb}}}$	100 μm, 1 mm, 10 cm
${a}_{{\rm{maxw}}}$	0.25 μm, 10μm

Note.

^aThe source SO697 needed a smaller ${R}_{{\rm{d}}}$ (Table 7).

Download table as:  ASCIITypeset image

We modeled 18 sources out of 27 in the sample considered here (SO540 is taken as a TD, see Section 4.1): 17 class II stars and one EV disk. One of these objects, SO844, has millimeter detections taken with SCUBA-2 at 850 μm, and the rest have only upper limits (Williams et al. 2013). Of the nine stars not modeled in this work, the source SO984 also has observations at 1300 μm taken with the Submillimeter Array (Williams et al. 2013). This object exhibits a MIPS 24 μm flux lower than that expected for class II stars (see Figure 3) and resembling that of PTDs, so it will be modeled, separately, in a future work. The rest are objects with significant uncertainties in their estimations of extinction and spectral type or objects without mass accretion rates (Section 4.2), so they were not modeled here.

Figure 12 shows the SEDs of our PACS disks (open circles) with the resulting fit (solid lines). Dashed lines indicate the σ Ori median (Table 4). Photospheric fluxes, using the colors of Kenyon & Hartmann (1995), are indicated by dotted lines. Table 7 gives the parameters of the best-fit model for each object with its corresponding reduced ${\chi }^{2}$ value (${\chi }_{\mathrm{red}}^{2}$) in the following order: target name, disk outer radius (R_d), confidence interval of disk outer radius (${R}_{{\rm{d}}}^{-}-{R}_{{\rm{d}}}^{+}$), degree of dust settling (), cosine of the inclination angle (μ), maximum grain size at the disk midplane (${a}_{{\rm{maxb}}}$), maximum grain size at the wall (${a}_{{\rm{maxw}}}$), disk mass (${M}_{{\rm{d}}}$), and confidence interval of disk mass (${M}_{{\rm{d}}}^{-}-{M}_{{\rm{d}}}^{+}$).

Zoom In Zoom Out Reset image size

Table 7.  Estimated Disk Parameters

Target	${R}_{{\rm{d}}}$	${R}_{{\rm{d}}}^{-}$–${R}_{{\rm{d}}}^{+}$		μ	${a}_{{\rm{maxb}}}$	${a}_{{\rm{maxw}}}$	${M}_{{\rm{d}}}$	${M}_{{\rm{d}}}^{-}$–${M}_{{\rm{d}}}^{+}$	${\chi }_{\mathrm{red}}^{2}$
	(au)	(au)				(μm)	(${10}^{-3}\,{M}_{\odot }$)	(${10}^{-3}\,{M}_{\odot }$)
SO1036	70	[69–94]	0.01	0.3	10 cm	10	5.19	[4.78–6.10]	6.116
SO1260	250	[180–250)	0.0001	0.9	100 μm	10	3.11	[2.57–3.11)	14.33
SO341	200	[158–208]	0.01	0.6	10 cm	10	1.24	[1.05–1.34]	24.784
SO362	150	[126–156]	0.0001	0.9	100 μm	0.25	1.01	[0.91–1.11]	31.489
SO396	10	(10–29]	0.01	0.3	100 μm	10	1.15	(1.15–2.46]	2.040
SO462	10	(10–16]	0.01	0.9	10 cm	0.25	0.098	(0.09–0.14]	265.259
SO562	200	[142–229]	0.0001	0.9	100 μm	10	2.70	[2.24–3.26]	55.817
SO615	200	[100–238]	0.0001	0.9	100 μm	10	1.25	[0.68–1.49]	51.825
SO662	200	[79–250]	0.001	0.9	10 cm	10	2.44	[1.32–2.91]	52.388
SO682	10	(10–12]	0.01	0.6	1 mm	0.25	0.031	(0.028–0.035]	14.469
SO697	7	[6–11]	0.01	0.6	10 cm	10	0.11	[0.10–0.16]	110.766
SO710	100	[63–126]	0.0001	0.9	100 μm	10	0.69	[0.45–0.85]	44.761
SO736	10	(10–12]	0.001	0.3	10 cm	0.25	0.22	(0.22–0.27]	27.975
SO774	250	[203–250)	0.01	0.9	10 cm	10	1.68	[1.43–1.68)	47.192
SO827	250	[185–250)	0.01	0.6	100 μm	10	4.94	[3.90–4.60)	34.874
SO844a	150	[149–158]	0.1	0.6	10 cm	10	39.40	[39.40–43.70]	48.159
SO859	30	[10–36]	0.01	0.3	100 μm	10	0.47	[0.23–0.57]	16.131
SO865	250	[237–250)	0.01	0.9	10 cm	0.25	0.85	[0.77–0.85)	3.450

Note. Column 1: ID following H07a; Column 2: disk outer radius; Column 3: confidence interval of disk outer radius; Column 4: degree of dust settling; Column 5: cosine of the inclination angle; Column 6: maximum grain size at the disk midplane; Column 7: maximum grain size at the disk wall; Column 8: disk mass; Column 9: confidence interval of disk mass; Column 10: ${\chi }_{\mathrm{red}}^{2}$ value of the fit.

^aObject with SCUBA-2 850 μm photometry from Williams et al. (2013).

Download table as:  ASCIITypeset image

In order to calculate confidence intervals of the outer radius of the disks, we first estimated the likelihood function, ${ \mathcal L }$, which is related to the ${\chi }^{2}$ values through the expression ${ \mathcal L }=\exp (-{\chi }^{2}/2)$. Since ${\chi }^{2}$ is a multidimensional function, at every radius we have several values of ${\chi }^{2}$, one for each one of the calculated models. Thus, the likelihood ${ \mathcal L }$ is computed using the minimum ${\chi }^{2}$ value at each radius. In Figure 13 we show plots of the normalized likelihood as a function of radius for all the 18 sources that were successfully modeled. We have restrained the x axis around the maximum peak in each case in order to better visualize the likelihood function. As usual, the confidence intervals of R_d are given as those radii ${R}_{{\rm{d}},1}$ and ${R}_{{\rm{d}},2}$ at which the area below the likelihood curve is 95% of its total area (Sivia & Skilling 2012). These intervals are indicated by light blue shaded regions in each panel and are reported in column 3 of Table 7. For those cases where the best radius falls on one of the edges of the range of radii used in the models, we have considered these values as upper or lower limits, and they are indicated by parenthesis instead of square brackets in Table 7. On the other hand, since the mass of the disk is a parameter that correlates with size, the confidence intervals of the mass are given as the values of the mass given by the best-fit model at ${R}_{{\rm{d}},1}$ and ${R}_{{\rm{d}},2}$. These values are reported in column 9 of Table 7.

Zoom In Zoom Out Reset image size

Some sources have peculiar SEDs: SO844 shows strong excess emission from the mid-IR to the far-IR, particularly in all four IRAC bands. This is probably due to variability. SO462, on the other hand, has a poor fit from H to IRAC 5.8 μm, which significantly increases the value of ${\chi }^{2};$ it is worth noting that this source has the largest ${A}_{{\rm{v}}}$ (=2.94, see Table 2). Finally, SO697 needed a smaller radius than the other sources (R_d = 7 au). This is due to its low flux at 70 μm combined with one of the lowest 160 μm upper limits. Therefore, any disk model with an outer radius greater than 7 au will overestimate the flux at 70 μm and will produce a flux greater than the upper limit at 160 μm.

We point out that we used the values of ${\chi }^{2}$ as a guide in order to obtain the model that provides the best fit to the photometry of each source. However, the values of ${\chi }^{2}$ are expected to be large. One of the most problematic points to fit is the IRAC 8 μm band. In a real spectrum, this band contains a contribution from the 10 μm silicate feature. The total flux of this feature is highly dependent on the adopted composition of the silicates, as detailed modeling of IRS spectra has shown (see McClure et al. 2012, 2013b). However, we felt that one photometric point did not give enough constraints to justify having the silicate composition as an additional parameter. Another reason for the large values of ${\chi }^{2}$ is that we are including only the errors in the photometry, which are small, and not other sources of uncertainty such as the inherent uncertainties in the distance, luminosity, spectral types, and mass accretion rates. In a forthcoming paper we will model these sources including their IRS spectra, in order to characterize their silicate features, which will improve the current fit. Through this study we will characterize the inner parts of the disk in the sense of grain growth, dust processing, and types of silicates present in the disks.

The flux in the mid-IR is mostly dictated by the mass accretion rate and the inclination, since for settled disks a large fraction of the flux arises in the innermost disk regions. For instance, 80% of the flux at 24 μm and shorter wavelengths arises inside ∼2 au for a disk with  = 0.001, instead of ∼20 au for a well mixed disk (D'Alessio et al. 2006; McClure et al. 2013a). In addition, as settling increases, the depleted upper layers become optically thin, exposing the disk midplane. Therefore, highly settled disks are colder and radiate less than disks with a low degree of settling (D'Alessio et al. 2006). We found that most of our objects (60%) can be explained by a significant degree of dust settling, $\epsilon =0.01$, consistent with studies in young regions like Taurus and Ophiuchus (Furlan et al. 2006; McClure et al. 2010). We also found that some of our objects have an even higher degree of dust settling ($\epsilon =0.0001$), which, in turn, is consistent with previous studies of disk fraction as a function of age in different young stellar populations (H07a).

The dependence of the flux on ${a}_{{\rm{maxb}}}$ has been discussed by D'Alessio et al. (2001). Even though we found nine sources to have ${a}_{\mathrm{maxb}}=100$ μm (Table 7), the left panel of Figure 14 shows that models with  = 0.01 and ${a}_{\mathrm{maxb}}$ $\leqslant$ 1 mm only differ by less than 10%, a factor that is even smaller for  = 0.0001, so we cannot discriminate sizes of ${a}_{{\rm{maxb}}}\leqslant 1$ mm. The right panel, on the other hand, shows the effects of changing the disk radius on the SED. As expected, smaller disks produced steeper mid- and far-IR slopes and lower fluxes. Differences in disks radii are more apparent for ${R}_{{\rm{d}}}\lt 70\,{\rm{a}}{\rm{u}}$. We found that 65% of our objects can be modeled with large sizes, ${R}_{{\rm{d}}}\geqslant 100\,{\rm{a}}{\rm{u}}$. The rest have dust disk radii less than 80 au. Note how changes in R_d affect the SED from MIPS 24 μm to millimeter wavelengths while changes in ${a}_{\mathrm{maxb}}$ are apparent only beyond 24 μm.

Zoom In Zoom Out Reset image size

We found no correlation between the mass accretion rate and the degree of dust settling. The disks with highest and lowest mass accretion rate in our sample have $\epsilon =0.01$. This seems to imply that the turbulence expected in the high accretors is not influencing the vertical distribution of dust in disks, which is contrary to expectations.

We have obtained detailed disk structures for 18 disks in the σ Ori cluster using irradiated accretion disk models (D'Alessio et al. 2006), constrained by mass accretion rates estimated independently from optical spectra. As a result of our modeling, we inferred disk masses and radii, assuming a gas-to-dust ratio close to that of the ISM. The masses of the best-fit models are given in Table 7; most of the masses are of the order of ∼10⁻³ to 10⁻⁴ ${M}_{\odot }$. By analyzing possible degeneracies in SED fitting using a set of 1800 generic models with different parameters, we estimated that the uncertainty in the inferred disk masses is a factor of 3, if the mass accretion rates of the sources are known.

Our mass determinations are consistent with those inferred by Williams et al. (2013) from 850 μm SCUBA-2 observations. Williams et al. (2013) detected nine disks, but only five of these lie inside the FOV of our PACS observations. Four of these sources (SO540, SO844, SO984, and SO1153) were detected at 70 and 160 μm. The only source not detected, SO609, is a class III star (H07a). Of the detections, SO1153 is a class I star, as seen in Figure 3, and it was not modeled here. Similarly, the source SO984 will be modeled in a future work together with other TDs such as SO540 (see Section 5.3). The SED and best-fit model for the source SO844 are shown in Figure 12 with inferred properties in Table 7. The disk mass for this object, 39 ${M}_{{\rm{Jup}}}$, is a factor of ∼8 higher than the mass estimated by Williams et al. (2013). The reason for this is that our dust grain opacity (${\kappa }_{\nu }$) at 850 μm for this object is significantly lower than theirs. Unlike Williams et al. (2013), which assumed a dust opacity of ${\kappa }_{\nu }=0.1(\nu /1200\,\mathrm{GHz})$ ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, we used a consistent opacity law (estimated through detailed modeling of the SEDs) for each one of our objects. This opacity depends on the mix of materials assumed in our disk models (silicates, graphite, water, etc.), their abundances, and their grain size distributions. Therefore, each object has a distinct disk dust opacity. With our dust opacity we can reproduce not only the flux at 850 μm but the entire SED (see Figure 12).

Figure 15 shows the disk mass distribution for our σ Ori sources in black compared to the mass distribution found by Williams et al. (2013) in red. Our detections include significantly lower values than the survey of Williams et al. (2013). The fact that we have detected 24 more sources than Williams et al. (2013) inside the FOV of their observations indicates that the PACS photometry was more suitable for detecting even small young disks than the shallow submillimeter observations. However, since our longest wavelength is 160 μm, we may be missing emission from the largest grains and underestimating our masses. Our models are consistent with the upper limits of the 850 μm SCUBA-2 observations (Figure 12), which suggests that the mass deficit is not large. Moreover, Williams et al. (2013) reported a 3σ limit of $4.3\times {10}^{-3}$ ${M}_{\odot }\sim 4.5$ ${M}_{\mathrm{Jup}}$, so their non-detections imply disk masses <5 ${M}_{\mathrm{Jup}}$, a value that is very similar to the upper limit of our disk masses for sources not detected with millimeter observations (∼5.5 ${M}_{\mathrm{Jup}}$, source SO1036; see Table 7). This supports the fact that we are truly looking at disks with very low masses. More sensitive millimeter observations are still required to determine this mass deficit. In any event, the masses of the σ Ori disks apart from SO844 range between 0.03 and ∼6 ${M}_{{\rm{Jup}}}$, with 35% being lower than 1 ${M}_{\mathrm{Jup}}$. Like Williams et al. (2013), we conclude that Jupiter-scale giant planet formation must be complete in these objects, which indicates either that giant planets form in less than 3 Myr or that it is difficult to make them in clustered environments.

Zoom In Zoom Out Reset image size

Figure 16 shows the disk mass distribution of σ Ori (our sample and Williams') compared to that of Taurus. An overall decrease in disk mass is clearly seen. This behavior of disk mass with age is expected from viscous evolution (Hartmann et al. 1998). For instance, between 1 and 3 Myr, the ages of Taurus and σ Ori, the disk mass should decrease by a factor of 1.5. However, the disk sizes should correspondingly increase by a factor of 2.5, which we do not see (see next section and Table 7). One possibility is that these disks have been subject to radial drift, according to which the millimeter-size grains have drifted inward and have piled up into pressure bumps, resulting in an observationally smaller dust disk (Pinilla et al. 2012b; Birnstiel & Andrews 2014), an effect that may be at play in TW Hya (Andrews et al. 2012) and other disks (Panić et al. 2009; Laibe et al. 2012; Pinilla et al. 2012a; Rosenfeld et al. 2013; Zhang et al. 2014). Another possible explanation for the presence of small disks is due to interactions with stellar companions. This effect has been observed in Taurus for close binaries with separations <40 au (Kraus & Hillenbrand 2012). Multiplicity in σ Ori has been reviewed by Caballero (2014), where he notes that outside the central arcminute, only 10 close binaries have been reported with angular separations between 04 and 30 (∼150–1200 au). Of these, only two have been imaged with adaptive optics. In short, more high-resolution imaging surveys of close binaries in σ Ori are needed in order to determine whether multiplicity is a major factor affecting the evolution of disks. Alternatively, these disks may have suffered the effects of photoevaporation, in which case both the dust and the gas have dissipated from the outer disk. We explore this possibility in the next section.

Zoom In Zoom Out Reset image size

The age of the σ Ori cluster is estimated to be ∼3 Myr (H07a), which is older than Taurus (Kenyon & Hartmann 1995, ∼1–2 Myr) or the ONC (Hillenbrand 1997, <1 Myr). By applying Kolmogorov–Smirnov tests for the n_24–70 and n_24–160 spectral indices between disks in Taurus and σ Ori (see Figure 11), we found that the levels of significance p for these indices are not small enough (compared to the statistics D) to say that the distributions are intrinsically different (Table 8). However, testing differences in the inner parts of the disks between the two distributions (using n_K–24, n_K–70, and n_K–160 spectral indices) shows that, as the wavelength increases, the distributions become substantially different, such that for the case of the last index, n_K–160, we can reject the null hypothesis that the samples are drawn from the same distribution ($D\gt p$, Table 8). This may imply that dissipation due to nearby massive stars dominates in the σ Ori cluster, photoevaporating the outer disks, in contrast to the inside-out dissipation in regions without massive stars, like Taurus.

Since disks in σ Ori are accretion disks, with substantial mass accretion rates (Section 4.2), viscous evolution must also be at play. Two alternatives are then possible: isolated viscous evolution or viscous evolution mediated by external photoevaporation. To consider the latter as a possibility, we must first consider its feasibility. To do so, we estimated the field intensity, G₀, produced by the three brightest objects in the system, the σ Ori Aa, Ab (O9.5, B0.5) pair and σ Ori B (∼B0–B2) with a total mass of ∼44 ${M}_{\odot }$ (Caballero 2014; Simón-Díaz et al. 2015). Using the reported stellar luminosities for these stars, and assuming that most of the stellar radiation is in the form of FUV photons, we estimated the intensity of the incident radiation at a given distance r from the source:

$$\begin{eqnarray}&&{G}_{0}=\displaystyle \frac{1}{{F}_{0}}\displaystyle \frac{{L}_{{\rm{FUV}}}}{4\pi {r}^{2}},\end{eqnarray} \tag{ 5 }$$

where F₀ is the typical interstellar flux level with a value of 1.6 × 10⁻³ erg s⁻¹ cm⁻² (Habing 1968), L_FUV is the FUV luminosity of the σ Ori multiple system and r the true distance to the ionizing sources. We used Monte Carlo simulations to calculate the distribution of true distances to the center for a given impact parameter, and from this the expected distribution of G₀, following Anderson et al. (2013).

We obtained that 80% of our disks are exposed to flux levels of $300\lesssim {G}_{0}\lesssim 1000$ (see Figure 17). Even though these are modest values of G₀, Anderson et al. (2013) show that even with values of external FUV radiation fields as low as ${G}_{0}=300$, external photoevaporation can play a significant role by greatly reducing the lifetime of the disks, as well as by truncating their outer edges. This is consistent with the fact that a third of our disks have radii ≤30 au (see Table 7). As we showed in Section 5.3 small radii can reliably be determined from the SED. These radii refer to the extent of the dust disk since we lack resolved gas observations; nevertheless, this result indicates that a possible explanation for the presence of small disks in populated regions is due to external photoevaporation by massive OB stars.

Zoom In Zoom Out Reset image size

Even though there is no correlation between disk radii and projected distance—a result also found in the ONC proplyds (Vicente & Alves 2005; Mann & Williams 2010; Mann et al. 2014)—our argument is reinforced by the fact that (a) two of our smallest disks, the sources SO697 (${R}_{d}=7$ au) and SO682 (${R}_{d}=10$ au), have the closest projected distances to the σ Ori multiple system (0.45 and 0.14 pc, respectively), (b) 81% of our largest disks (e.g., the sources SO774, SO865, and SO1260) are located at projected distances greater than 1 pc, and (c) most of the 850 μm SCUBA-2 detections of Williams et al. (2013) are outside the PACS field, at distances $\geqslant 2.23$ pc from the central ionizing sources.

In Figure 18 we show our determinations of disk masses and radii in σ Ori compared to those of proplyds in the ONC. We have included in the figure the confidence intervals of ${M}_{{\rm{d}}}$ and R_d estimated in Section 5.3. The disk masses for the ONC are taken from Mann & Williams (2010) and radii from Vicente & Alves (2005). We also show the expected evolution of viscous disks subject to photoevaporation from Anderson et al. (2013). Unlike isolated viscous evolution, in which disks expand as the disk mass increases, viscous disks subject to external photoevaporation first expand until they reach the evaporation radius, where photoevaporation begins to act and to rip out the external disk regions. As a result, both the disk mass and radius decrease with age (Anderson et al. 2013). The σ Ori disks fit this trend. Their masses are lower than those of the ONC disks, since they are older, while their radii are comparable or smaller. The two models shown in Figure 18 correspond to two values of G₀, 300 and 3000, with the lower value corresponding to overall larger disk radii. This suggests that the shift between the ONC disks and the σ Ori disks may be due to the latter being subject to comparatively lower levels of FUV flux (which explains the slightly larger radii) for a longer period of time (which may explain the smaller radii). However, given that disk radii of proplyds were measured from Hubble Space Telescope/WFPC2 Hα images, and hence these values are related to the gaseous disks, and we estimated disk radii from modeling the SEDs, which are associated with the dusty disks, there are likely systematic differences between the two measurements of disk radius made using very different methods. We need more sensitive millimeter observations of gas in these disks in order to make a more robust comparison.

Zoom In Zoom Out Reset image size

The evolutionary models shown in Figure 18 assume a constant value of G₀ throughout the disk evolution. This in only an approximation since stars are not static—they orbit around the cluster center and therefore will experience different FUV flux levels during their lifetime. For very eccentric orbits the stars will stay longer in outer regions, where the FUV field is low, and retain their disk masses longer. Furthermore, these models assume a constant value of α; adopting a non-uniform radial profile for the viscosity alters the disk structure (Kalyaan et al. 2015) and therefore the mass–radius evolutionary tracks. We adopted a viscosity parameter of α = 0.01. With this viscosity the disks are depleted down to 10^–5${M}_{\odot }$ very quickly (<1 Myr). Using a lower α will increase the lifetime of the disks but will not produce disks as large as observed. Note that the timescale for disk dispersal depends strongly on the assumed disk viscosity, and therefore a relatively small range of α can lead to a wide range in the disk mass and radius at a given time. The models also assume an initial disk radius of 30 au. In reality, disks will exhibit a wide range of radii depending on the initial conditions of the core where the stars are formed. However, Anderson et al. (2013) showed that varying the initial radius by a factor of two does not change the evolutionary tracks. Another approximation of the models is the assumption of a stellar mass of ${M}_{* }=1\,{M}_{\odot }$. Stars with smaller masses will produce shallower gravitational potential wells, allowing the unbound material to be located closer to the star, and possibly producing smaller disks. All these approximations may be responsible for the differences in radii observed between the models and disks in σ Ori. Finally, as mentioned above, the radial distribution of dust may differ from that of the gas. However, even with these simplifying assumptions, the models reproduce the observations of the ONC remarkably well and those of the σ Ori disks sample quite well.

Additional support for the photoevaporative hypothesis comes from the presence of forbidden lines in the optical spectra of the disk sources. These lines are expected to form in a photoevaporating wind (Hartigan et al. 1995; Acke et al. 2005; Pascucci & Sterzik 2009; Rigliaco et al. 2009, 2012, 2013; Gorti et al. 2011; Natta et al. 2014). Rigliaco et al. (2009) found forbidden lines of [S ii] and [N ii] in high-resolution spectra of SO587, located at d_p = 0.35 pc. They associated these emission lines with a photoevaporating wind estimated to give ${\dot{M}}_{\mathrm{loss}}\sim {10}^{-9}$ ${M}_{\odot }\,{\mathrm{yr}}^{-1}$. Unfortunately, we did not model this object because it was not detected with PACS.

The survey of H14 includes high-resolution spectra of some of the PACS sources; these spectra cover the region around Hα and include the [N ii] λ6548.05 and λ6583.45 forbidden lines. The top panel of Figure 19 displays the spectrum of SO662 (R_d = 200 au, M_d ∼ 2 M_Jup; Table 7), showing the Hα line in emission, characteristic of TTSs, and the forbidden lines of [N ii]. This object is located at d_p = 0.6 pc. In the bottom panel we also show the lines in velocity space. The blueshift of these lines, if any, is smaller than 3 km s⁻¹ (our error in velocity determination), consistent with the lowest blueshifts found in TTSs for the [O i] line (Rigliaco et al. 2013; Natta et al. 2014). The low velocity value indicates that the forbidden line emission does not come from a jet, for which velocities reach hundreds of km s⁻¹ (Hartigan et al. 1995). Figure 20, on the other hand, shows the spectrum of SO697, our smallest disk (R_d = 7 au, M_d ∼ 0.1 M_Jup; Table 7) located at d_p = 0.45 pc. This object does not exhibit the [N ii] forbidden line at 6548.05 Å and has a very weak emission at 6583.45 Å (bottom right panel). Similarly, none of our small disks with high-resolution spectra (SO396, SO462, and SO859) shows any emission lines of [N ii].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Strong optical forbidden emission lines are expected to form in highly ionized regions, such as in the bright cusps characteristic of proplyds in the ONC, so the fact that within the first 0.6 pc from σ Ori, one of our largest disks exhibits these lines while none of our small disks has any features suggests that photoevaporation is in action in the cluster. If this is indeed the case, then SO662 is a photoevaporating disk in its initial phase while the small disks are those that moved close enough to the hot stars during the lifetime of the cluster to be stripped of their outermost regions and hence cannot produce the forbidden lines. We need a greater sample of disks with high-resolution spectra covering a range of projected distances from the central stars to confirm this hypothesis and characterize the region in the cluster that is affected by photoevaporation.