GAS INSIDE THE 97 AU CAVITY AROUND THE TRANSITION DISK Sz 91

Sz 91 was observed with ALMA band 6 (231 GHz) and band 3 (110 GHz) during Cycle 0 (program 2011.0.00733, PI: M.R. Schreiber). The observations are summarized in Table 1.

Sz 91 was observed in five different epochs with different weather conditions and antenna configurations for the band 6 data. The system temperature ranged from 53 to 96 K. For the band 3 data, the ALMA correlator was configured to provide two spectral windows centered on the continuum and two spectral windows centered on the $^{13}$CO(1–0) (110.20135 GHz) and the C¹⁸O(1−0) (109.78218 GHz) lines. The total bandwidth of the band 3 observations was 7.5 GHz ($4\times 1.875$ GHz). For the band 6 observations, the ALMA correlator was initially configured to provide one spectral window centered on the continuum, and three spectral windows centered on the ¹²CO(2−1), 230.538 GHz (hereafter ¹²CO), ¹³CO(2−1), and C¹⁸O(2−1) lines. Unfortunately, because of technical problems only one sideband of the correlator could be configured, and only the spectral windows centered on the continuum and on the ¹²CO line could be produced. The total bandwidth of the band 6 observations was $2\times 1.875$ GHz. Both ALMA bands were sampled at 0.488 MHz (e.g., 0.635 km s⁻¹ in the ¹²CO line).

The quasar QSO B1730-130 and the active galactic nucleus ICRF J160431.0-444131 were used as bandpass and primary phase calibrators in both bands, respectively. Neptune and Titan were used to calibrate amplitude in band 3 and band 6, respectively. The observations of the calibrators were alternated with the science observations. The individual exposure times for the calibrators and the science targets was 6.05 s, amounting to a total exposure time for the science observations of 226 s (3.77 minutes) in band 6 and 24 s (0.4 minutes) in band 3. Phase correction based on WVR measurements was performed in offline mode as part of the basic ALMA Cycle 0 corrections. We used the dispersion of the band 6 flux calibrations to estimate a flux uncertainty of $\approx 15\%$. This value is a lower limit since it does not include uncertainties from the amplitude calibration.¹²

The observations were processed using the Common Astronomy Software Application (CASA) package (McMullin et al. 2007). The visibilities were Fourier transformed and deconvolved, using natural weighting, with the CLEAN algorithm (Högbom 1974). For the band 6 data, we combined the four data sets to increase the UV coverage before cleaning. The cleaned image produced from this concatenated data set had a lower rms than the images created from individually cleaned data sets. After combining the band 6 data, the two bands contained 276 baselines, ranging from 21 to 402 m (16 to 309 kλ) in band 6 and from 21 to 452 m (7 to 167 kλ) in band 3.

In the continuum, we reach an rms of 0.1 mJy beam⁻¹ in band 6 and 0.1 mJy beam⁻¹ in band 3. The median rms on the channels associated with ¹²CO is 6.2 mJy beam⁻¹. In the continuum, the beam is $0\buildrel{\prime\prime}\over{.} 86\times 0\buildrel{\prime\prime}\over{.} 60$ with a position angle (PA) of 889 in band 6 and $2\buildrel{\prime\prime}\over{.} 08\times 1\buildrel{\prime\prime}\over{.} 51$ with a PA of 994 in band 3. No emission was detected in either the ¹³CO(1−0) or the C¹⁸O(1−0) lines in band 3.

We add two WISE fluxes at 12 and 22 μm; two Herschel/PACS fluxes at 100 and 160 μm; three Herschel/Spire fluxes at 250, 350, and 500 μm; and the two millimeter fluxes discussed here to the photometry presented by Romero et al. (2012). For the Herschel fluxes we adopt the most recent values presented by Bustamante et al. (2015). The fluxes at wavelengths shorter than 24 μm were corrected for extinction by applying the dereddening relations listed in Cieza et al. (2007). We assume calibration errors of 20% for the optical photometry and 15% for the photometry up to 24 μm. The photometry between 100 and 500 μm is affected by background contamination from the cloud in which Sz 91 is embedded, particularly at 250, 350, and 500 μm (Bustamante et al. 2015). We therefore consider the Herschel fluxes at 250, 350, and 500 to be upper limits. For the 870 μm flux we used the 8.5% uncertainty listed by Romero et al. (2012). To derive the fluxes from our ALMA observations we perform a fit in the uv plane. Because of the lack of data at short baselines (i.e., large spatial scales), deriving the total flux from the cleaned images would result in underestimated values. Instead, we fit different profiles to the visibilities to estimate the flux at the center of the uv plane. In band 6 Sz 91 is resolved (see the next section), so we fit a Gaussian profile; in band 3 Sz 91 is not resolved, so we use a point-source profile. The fluxes comprising the SED, their errors, and references are listed in Table 2.

First estimations of the Sz 91 stellar parameters were obtained by Hughes et al. (1994), who used a distance of 140 ± 20 AU, an extinction A_v of 2, and a spectral type of M0.5. These authors derived a stellar mass range of ${{M}_{\star }}=0.49{\rm --}0.69\;{{M}_{\odot }}$, an effective temperature of ${{T}_{{\rm eff}}}=3723$ K, and an age ranging from 5 to 7 Myr. Romero et al. (2012) classified Sz 91 as an M1.5 star using high-resolution spectra, in general agreement with the previous result. Using the latter spectral type, the most recent distance estimate of ∼200 pc, and an extinction of A_v = 2, we derive the stellar parameters using the stellar evolutionary models of Siess et al. (2000). We obtain a stellar radius of ${{R}_{\star }}=1.46\;{{R}_{\odot }}$, a temperature of ${{T}_{{\rm eff}}}=3720$ K, a mass of ${{M}_{\star }}=0.47\;{{M}_{\odot }}$, and an age below 1 Myr.¹³ In what follows we use these values.

At millimeter wavelengths, we can use the Rayleigh–Jeans approximation to express the flux as ${{F}_{\nu }}\propto {{\nu }^{{{\alpha }_{{\rm mm}}}}}$. Assuming optically thin emission, the millimeter slope is a function of the dust opacity index β, i.e., ${{\alpha }_{{\rm mm}}}=2+\beta$. In particular, $\beta \sim 2$ for the interstellar medium (ISM; see Williams & Cieza 2011 and references therein). Trying to fit the fluxes at 0.85, 1.3, and 2.7 mm in this way, it is not possible to match the observed fluxes, as shown in Figure 3. The slope derived from a fit to the three fluxes results in ${{\alpha }_{0.8{\rm --}2.7\,{\rm mm}}}=3.36\pm 0.14$ (solid line). Using the two shorter wavelengths we obtain ${{\alpha }_{0.8{\rm --}1.3\,{\rm mm}}}=2.34\pm 0.40$ (dotted line), whereas using the fluxes at 1.3 and 2.7 mm yields ${{\alpha }_{1.3{\rm --}2.7\,{\rm mm}}}=4.07\pm 0.29$ (dashed line). This large ${{\alpha }_{1.3{\rm --}2.7\,{\rm mm}}}$ would imply moderate or very low grain growth. On the other hand, our derived ${{\alpha }_{0.8{\rm --}1.3\,{\rm mm}}}$ implies ${{\beta }_{0.8{\rm --}1.3\,{\rm mm}}}=0.34\pm 0.40$, which suggests the opposite. The latter value agrees with the result of T2014 ($\beta =0.5\pm 0.1$) and with the average values for a set of transition disks and protoplanetary disks derived by Pinilla et al. (2014). In addition, we conclude from our models (see Section 5) that grains of at least 1 mm in size are needed to match the SED at millimeter wavelengths, implying that grain growth is happening in Sz 91 (see Section 5). One possible explanation of the different slopes is an optical depth effect. Enhancements in the optical thickness could reduce the slope, explaining the reduction in ${{\alpha }_{0.8{\rm --}1.2{\rm mm}}}$ compared to ${{\alpha }_{1.3{\rm --}2.7\,{\rm mm}}}$. However, using a sample of 50 disks, Ricci et al. (2012) find that low values of ${{\alpha }_{{\rm mm}}}$ can be explained by emission from optically thick regions only in the most massive disks of their sample. Given the relatively low mass of Sz 91 (Section 5), the good agreement with the results by T2014, and the significant difference between ${{\alpha }_{0.8{\rm --}1.3{\rm mm}}}$ and ${{\alpha }_{1.3{\rm --}2.7{\rm mm}}}$, we consider that an underestimation of the error bars represents a more likely explanation of the broken ${{\alpha }_{{\rm mm}}}$ slope and that ${{\alpha }_{0.8{\rm --}1.3{\rm mm}}}=2.34\pm 0.40$ probably represents a realistic estimate of the true millimeter slope.

Zoom In Zoom Out Reset image size

Given the strong absorption observed toward the redshifted channels, deriving the disk's center and geometry by means of a Keplerian fit to the line profile as in, e.g., Mathews et al. (2012) would result in large uncertainties. Instead, we determine the disk inclination (i), position angle (PA), and center by fitting different disk profiles to the continuum visibilities using the uvfit/MIRIAD task.

A ring morphology well matches our continuum observations and the results of T2014. Unfortunately, uvfit does not allow us to adjust the thickness of the fitted ring. To test for the effect of using different disk profiles, we performed three fits using a ring, a disk, and a Gaussian profile. Combining these results, we derived a PA of 182 ± 3.5 (east of north) and an inclination (i) of 495 ± 3.5. These results are in agreement with those derived by T2014 (PA = 175 ± 177, $i=40{}^\circ \pm 15{}^\circ$). We find the disk's center to have an offset of ${\Delta }\alpha =-0\buildrel{\prime\prime}\over{.} 29\pm 0\buildrel{\prime\prime}\over{.} 01$ and ${\Delta }\delta =-0\buildrel{\prime\prime}\over{.} 06\pm 0\buildrel{\prime\prime}\over{.} 02$ with respect to the center of the ALMA field. This offset is consistent with the measured proper motion of (−20.3, −18.4) mas ${\rm y}{{{\rm r}}^{-1}}$ (UCAC3, Zacharias et al. 2010) when including the ∼01 astrometric uncertainty from our ALMA Cycle 0 observations and the 006 uncertainty from our input coordinates (2MASS catalog, Cutri et al. 2003). We corrected for this offset in our analysis of the disk.

A quick inspection of the continuum (Figure 2, left panel) and the ¹²CO moment 0 and moment 1 (Figure 2, center panel) images shows that the gaseous disk extends farther out than the continuum disk. The resolution of our images along the major axis of the disk is set by the beam size in that direction, i.e., $0\buildrel{\prime\prime}\over{.} 61$. Keeping this in mind, we produced a 3 pixel wide cut along the blueshifted (i.e., less affected by cloud emission) semimajor axis of the disk for the continuum and ¹²CO moment 0 images (Figure 4). We compute the fluxes every $0\buildrel{\prime\prime}\over{.} 16$, which means that they are correlated (i.e., their separation is smaller than the beam size). Still, we can use them to obtain information about the extension of the gaseous and dusty disks. In the plot, the vertical dotted lines indicate the position at which the emission becomes significant (i.e., above $3\sigma$). The continuum emission becomes undetectable beyond $\sim 1\buildrel{\prime\prime}\over{.} 1$ (i.e., ∼220 AU) from the center, whereas the ¹²CO emission is still detected at $\sim 2^{\prime\prime}$ (∼400 AU), i.e., almost 1.5 beams farther out than the dusty disk. We therefore conclude that, although we are limited by the moderate resolution of the beam size, the gaseous disk is undoubtedly detected at a larger distance from the star than the dusty disk. The outer radius of the dusty and gaseous disk estimated by T2014 at 345 GHz (170 ± 20 AU for the continuum, 420 AU for the gas) matches particularly well with our results. Our limited spatial resolution does not allow us to probe for the sharp outer edge in the continuum disk observed in other disks, e.g., TW Hya (Andrews et al. 2011) and HD 163296 (de Gregorio-Monsalvo et al. 2013), that could be a consequence of radial drift (Birnstiel et al. 2013).

Zoom In Zoom Out Reset image size

The shape of the continuum cleaned image is indicative of the presence of an inner cavity. To explore the partially resolved inner disk we here focus on the radially averaged deprojected ALMA visibilities. We split the continuum and ¹²CO visibilities into 30 radial bins, ensuring that all bins contain the same amount of visibilities. This results in minimum and maximum bin sizes of 3.7 and 22 kλ, respectively, and represents a good compromise between UV sampling and error bars. The real part of the binned visibilities of our band 6 observations as a function of UV distance are shown in Figure 5. The steeper slope of the ¹²CO visibilities shows that the detected gaseous disk is indeed more extended than the dusty disk, as already noted from the cleaned images. The null in the continuum visibilities at $\sim 150\;{\rm k}\lambda$ indicates a discontinuity caused by the lack of millimeter-sized particles in the inner disk regions, or, in other words, the presence of a cavity (see Hughes et al. 2007, 2009, Andrews et al. 2011). With our observations we cannot estimate the sharpness of the inner wall (see, e.g., Mulders et al. 2013). In the next section we derive the cavity size in the continuum by fitting the uv plane.

Zoom In Zoom Out Reset image size

The ¹²CO visibilities have lower S/Ns than the continuum and they are severely affected by the cloud. This complicates the description of the gaseous disk in terms of the analysis of the visibility profile. However, if the gaseous disk velocities are governed by Keplerian rotation, we can use the ¹²CO line profile to estimate how close to the central star we detect gas using the projected observed velocity $v{{(r)}_{{\rm Kep}}}=\sqrt{G{{M}_{\star }}/r}{\rm sin} \;i.$ Assuming that the line center is at ${{v}_{{\rm LSR}}}=3.4$ km s⁻¹, we detect ¹²CO emission down to ${{v}_{{\rm LSR}}}=0.5$ km s⁻¹ with an S/N of ∼6 (Figure 2). This corresponds to −2.9 km s⁻¹ in the star's reference system. Using this velocity, the stellar mass and inclination previously derived, and a distance of 200 pc, we conclude that our observations probe the ¹²CO gas down to ∼28 AU from the central star.

To account for the SED shape at mid-infrared wavelengths, we define two inner disk regions: the innermost region extending from the dust sublimation radius to ${{R}_{{\rm out},1}}$ which is severely dust-depleted ("Region 1") and accounts for the fluxes at 12 μm and below, and the intermediate region extending from ${{R}_{{\rm in},2}}$ to ${{R}_{{\rm out},2}}$ ("Region 2") where some dust is emitting at 22 and 24 μm. There is no observational constraint on the surface density distribution of these two regions, and therefore we adopt a simple power law of the form

$$\begin{eqnarray}&&{\Sigma }(r)={{{\Sigma }}_{100}}{{\left( \frac{r}{100\;{\rm AU}} \right)}^{-{{p}_{s}}}},\end{eqnarray} \tag{ 1 }$$

where ${{{\Sigma }}_{100}}$ is the surface density distribution at 100 AU and p_s is the power-law index. We define the outer disk, "Region 3," as the part of the disk extending from the cavity radius (${{R}_{{\rm cav}}}$) outward. For a number of disks where the gaseous disk is detected beyond the dusty disk, as in Sz 91, a tapered-edge profile can naturally explain the observations (e.g., Hughes et al. 2008; Andrews et al. 2011; Williams & Cieza 2011; Cieza et al. 2012; Mathews et al. 2012; Zhang et al. 2014). There are three exceptions to this behavior. To explain the high-resolution observations of IM Lup, TW Hya, and V4046 Sgr, strong variations of the gas-to-dust ratio at large radii are required (Panić et al. 2009; Andrews et al. 2012; Rosenfeld et al. 2013). Our data do not have the high spatial resolution and sensitivity needed to test for radial changes in the gas-to-dust ratio and therefore we adopt the tapered-edge prescription to describe the outer disk. The surface density distribution in Region 3 is then defined as

$$\begin{eqnarray}&&{\Sigma }(r)={{{\Sigma }}_{C}}\;{{r}^{-\gamma }}\;{\rm exp} \left[ -{{\left( \frac{r}{{{R}_{C}}} \right)}^{2-\gamma }} \right],\end{eqnarray} \tag{ 2 }$$

where R_C is the characteristic radius that defines where the tapered edge begins to dominate over the power-law component, ${{{\Sigma }}_{C}}$ is the surface density at R_C, and γ corresponds to the viscosity power-law index in accretion disk theory (ν ∝ R^γ; Lynden-Bell & Pringle 1974; Hartmann et al. 1998). For each region we define the scale height H(r) assuming that the dust follows a Gaussian vertical density profile

$$\begin{eqnarray}&&H(r)={{H}_{0}}{{(r/100\;{\rm AU})}^{\psi }},\end{eqnarray} \tag{ 3 }$$

where ψ is the flaring parameter of the disk and H₀ is the scale height at r = 100 AU.

We assume homogeneous spherical dust particles. The scattering opacities and phase functions, extinctions, and Mueller matrices are computed using the Mie theory. For the composition we assumed compact astro-silicates (Mathis & Whiffen 1989) with a power-law size distribution ($dn(a)\propto {{a}^{-p}}da$, with p = 3.5; e.g., Draine 2006). This assumption requires minimum and maximum grain sizes of ${{a}_{{\rm min} }}=2\;\mu {\rm m},{{a}_{{\rm max} }}=15\;\mu {\rm m}$ in Region 1 to not overproduce the flux at 8 and 12 μm due to the silicate resonance bands. We note that using a different grain composition, e.g., carbon grains, it is also possible to match the SED using smaller values of ${{a}_{{\rm min} }}$ inside the cavity. In any case, the lack of emission above photospheric levels below $8\;\mu {\rm m}$ indicates that the innermost region of the cavity must be severely depleted of dust particles. Region 2 is mainly constrained by the fluxes at 22, 24, and 70 μm and therefore many combinations of a_min, a_max and the dust mass can reproduce the observed fluxes. The minimum grain size, however, cannot exceed ${{a}_{{\rm min} }}=0.1\;\mu {\rm m}$ in order to match the observed fluxes.

We use the stellar parameters derived in Section 4.1 and the inclination (i) and PA derived in Section 4.3. We focus on identifying a model that can describe the thermal structure of the disk (constrained by SED and continuum observations) using the disk structure and dust parameters described in the previous sections. Regions 1 and 2 are poorly constrained and therefore we fix the power-law indexes of their surface density distribution to the standard value of p_s = 1 following Andrews & Williams (2007). In Region 1 we fix the minimum and maximum grain sizes to the values discussed in the previous section. Using this prescription we find that we need a maximum dust mass of ${{M}_{{\rm dust},1}}=1\times {{10}^{-4}}\;{{M}_{\oplus }}$ to reproduce the SED shape below 12 μm. We adopt this value to describe the dust content inside Region 1. A first exploration of the H₀ parameter in this region reveals that it is completely unconstrained, and we adopt an arbitrary value of ${{H}_{{{0}_{1}}}}=5{\rm AU}.$ In Regions 1 and 2 we fix the minimum grain size to ${{a}_{{\rm min} }}=0.05\;\mu {\rm m}$, as in the ISM. For Region 3 we explore the effects of varying the surface density parameter γ within the range of $[-1,1]$. We find that our data are equally well fitted with values ranging from $\gamma =[-0.2,0.8]$ and use $\gamma =0.3$. Similarly, our data do not provide observational constraints on the flaring index (ψ) in any of the three regions. For all three regions we adopt a typical value of $\psi =1.15$ from model results in, e.g., T Cha (Cieza et al. 2011), IM Lup (Pinte et al. 2008), and HD 163296 (de Gregorio-Monsalvo et al. 2013). The free parameters in our model approach are the outer radius of Regions 1 and 2 (${{R}_{{\rm ou}{{{\rm t}}_{1,2}}}}$), the inner radius of Region 2 (${{R}_{{\rm in},2}}$), and the cavity and characteristic radius of Region 3 (${{R}_{{\rm cav}}},{{R}_{C}}$). We also fit the dust masses ${{M}_{{\rm dus}{{{\rm t}}_{2,3}}}}$, the maximum grain sizes ${{a}_{{{{\rm max} }_{2,3}}}}$, and the scale heights ${{H}_{{{0}_{2,3}}}}$.

We explore the parameter space by means of a genetic algorithm. The procedure starts by randomly selecting a number of models. The parameters of these models are defined as their "genes." Models producing a better match to our observational data set have the highest chance of reproducing their parameter values (genes) into the next generation of models. The quality of a fit is measured with the combined reduced ${{\chi }^{2}}$ given by the sum of reduced ${{\chi }^{2}}$ from SED and the continuum visibilities, i.e., ${{\chi }^{2}}=\chi _{{\rm SED}}^{2}+\chi _{{\rm UV}}^{2}$. We use a population of 50 "parent" models and let them evolve through 50 generations, resulting in a total of 2500 models. We find that after 35 generations the models quickly converge, and after 45 generations the variations of each free parameter are within 5% of the average value of the parameter.

We find a best-fit model that agrees well with our ALMA observations and with the SED of Sz 91. Its parameters are summarized in Table 3. The most significant result of our modeling is that a large cavity of ∼97 AU in the dust grain distribution is required to mach our observations. This is clearly illustrated by the null in the real part of the continuum deprojected visibilities as shown in Figure 5 (left panel). For direct comparison with our ALMA band 6 images we created synthetic ALMA observations using the simalma/casa ¹⁴ package. We created an antenna configuration file that mimics our observations using the buildConfigurationFile/Analysis Utilities ¹⁵ task. We then run the simulation using the same PWV, exposure time, and hour angle as in our observations. The synthetic ALMA image of our model for Sz 91 and the residuals obtained by subtracting the synthetic image from the real observations are shown in Figure 6 (center and right panels). The residual image has an rms of 0.1 mJy beam⁻¹, similar to the rms of the observations. Interestingly, our simulated image also shows a brightness difference between the northern and southern lobes of the disk, with the northern lobe being the brighter one (although our MCFOST model image is axisymmetric). To investigate this feature we repeated the SIMALMA/CASA simulations using different UV coverages and PWV. We find that the asymmetry is a consequence of the low S/N of the visibilities at longer baselines, as increasing the UV coverage at larger baselines or decreasing the PWV of the simulated ALMA observations removes the apparent difference in brightness.

Zoom In Zoom Out Reset image size

Table 3.  Model Results

Parameter	Value
${{T}_{\star }}$	3720 K
${{R}_{\star }}$	1.46 ${{R}_{\odot }}$
${{M}_{\star }}$	0.47 ${{M}_{\odot }}$
d	200 pc
i	495
PA	182
p	3.5
${{a}_{{{{\rm min} }_{1}}}},{{a}_{{{{\rm max} }_{1}}}}$	$2,15$ μm
${{a}_{{{{\rm min} }_{2,3}}}}$	$0.05\;\mu {\rm m}$
${{H}_{{{0}_{1}}}}$	5 AU
${{p}_{{{s}_{1,2}}}}$	1
γ	0.3
${{\psi }_{1,2,3}}$	1.15
${{M}_{{\rm dus}{{{\rm t}}_{1}}}}$	1 × 10−4 ${{M}_{\oplus }}$
${{R}_{1}}$	0.025 AU
${{a}_{{{{\rm max} }_{2}}}}$	$5\;\mu {\rm m}$
${{a}_{{{{\rm max} }_{3}}}}$	1000 μm
${{R}_{{\rm ou}{{{\rm t}}_{1}}}}$	85 AU
${{R}_{{\rm i}{{{\rm n}}_{2}}}}$	85 AU
${{R}_{{\rm ou}{{{\rm t}}_{2}}}}$	97 AU
Rcav	97 AU
RC	100 AU
${{M}_{{\rm dus}{{{\rm t}}_{2}}}}$	0.7 ${{M}_{\oplus }}$
${{M}_{{\rm dus}{{{\rm t}}_{3}}}}$	14.3 ${{M}_{\oplus }}$
${{H}_{{{0}_{2}}}}$	5 AU
${{H}_{{{0}_{3}}}}$	10 AU
${\rm log} ({{{\Sigma }}_{{{100}_{1}}}})$	−8.20 (${\rm gr}\;{\rm c}{{{\rm m}}^{-2}}$)a
${\rm log} ({{{\Sigma }}_{{{100}_{2}}}})$	−2.95 (${\rm gr}\;{\rm c}{{{\rm m}}^{-2}}$)a
${\rm log} ({{{\Sigma }}_{C}})$	−1.98 (${\rm gr}\;{\rm c}{{{\rm m}}^{-2}}$)a

Notes. Fixed parameters are listed above the line.

^aDerived values.

Download table as:  ASCIITypeset image

The SED and surface density distribution of the model are shown in Figures 7 and 8, respectively. In Figure 7 we additionally show the individual contributions of Regions 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To estimate which model parameters are well constrained from our observational data set we proceed as follows. We fix all but one parameter and run a grid of models exploring the local parameter space. We then derive the reduced ${{\chi }^{2}}$ of the models forming this local grid. This procedure is repeated for all free parameters. Keeping in mind that this is equivalent to assuming that all parameters are uncorrelated (which is clearly not the case), this approach allow us to get a first idea of which parameters are clearly constrained by our observations without spending unreasonable amounts of computing time. The results of this exploration are shown in Figure 9. These plots illustrate the severe degeneracy of most of the parameters describing Regions 1 and 2 in our model disk. In contrast, we find that the parameters such as M_dust, a_max, and the cavity size (R_cav) in Region 3 are more constrained by our observations. In particular, values of R_cav below 90 AU are very difficult to reconcile with the ALMA observations. According to our model, most of the disk's mass (10 ${{M}_{\oplus }},\sim 70\%$ of ${{M}_{{\rm dust},3}}$) is concentrated in a ring ranging from 97 to 140 AU.

Zoom In Zoom Out Reset image size

Our observations trace the optically thick ¹²CO line, which is severely affected by the cloud. We therefore only present a simple model of the gas component of the disk around Sz 91 that is able to match the ¹²CO line profile reasonably well using standard assumptions. We only consider the blueshifted (less absorbed) side of the line. Given the sparse observational constraints on the gaseous disk, we build this model upon the best-fit model for the dusty disk obtained in the previous section. The lack of measurements of the isotopologues ¹³CO and C¹⁸O prevents us from estimating the gas-to-dust ratio and total gas mass as in Williams & Best (2014). We fix the gas-to-dust ratio of the outer region (Region 3) to the standard value of 100, and we ensure that the ¹²CO freezes out at 20 K (Qi et al. 2004; de Gregorio-Monsalvo et al. 2013). We assume that the gas and the dust have the same temperature (local thermal equilibrium, LTE). For the ¹²CO abundance we use the standard ${{{\rm H}}_{2}}$/¹²CO ratio (10⁴), which is a reasonable value for protoplanetary disks as confirmed by France et al. (2014). The gas is assumed to be in Keplerian rotation with an inner radius of at least 28 AU (in agreement with our observations). Assuming that the gas is in hydrostatic equilibrium in this direction, the gas profile is described by

$$\begin{eqnarray}&&\rho (r,z)=\rho (r,0){\rm exp} \left[ -\frac{{{z}^{2}}}{2H{{(r)}^{2}}} \right].\end{eqnarray} \tag{ 4 }$$

The scale height of the gas component H(r) is given by

$$\begin{eqnarray}&&H(r)=\sqrt{\frac{{{r}^{3}}{{k}_{{\rm B}}}T(r)}{G{{M}_{\star }}\mu {{m}_{H}}}},\end{eqnarray} \tag{ 5 }$$

where ${{k}_{{\rm B}}}$ is the Boltzmann constant, G is the gravitational constant, and m_H is the hydrogen mass. We also take into account the thermal broadening of the velocities (defined as ${{V}_{{\rm th}}}=\sqrt{2{{k}_{{\rm b}}}{{T}_{{\rm CO}}}/{{m}_{{\rm CO}}}}$, where m_CO is the molecular mass of ¹²CO) and the effect of turbulence (v_turb) on the ¹²CO emission along the line of sight.

Under these assumptions, we manually explore different models by varying the turbulence velocity (v_turb) within a range of [10,300] m s⁻¹ (in accordance with the observations discussed by Hughes et al. 2011) and the gas mass inside the dust-depleted regions (Regions 1 and 2). We find a model that agrees relatively well with our observations using a continuous gaseous disk (i.e., without depletion or an inner cavity) with ${{v}_{{\rm turb}}}=120$ m s⁻¹. The radially binned deprojected visibilities and the surface density of this model are shown in Figure 5 (right panel) and Figure 8 (red line), respectively. The integrated ¹²CO line profile of the model, along with the real observations, is shown in Figure 10.

Zoom In Zoom Out Reset image size

Assuming three different disk regions, we find a model that reproduces our observational data set reasonably well. As a separate exercise we also tried to fit a disk model composed of only two regions (the inner cavity and the outer disk). After running the genetic search we could not find a good match to our observations with such a disk structure. This result is in agreement with the findings of T2014, who need an extra component, either a circumplanetary disk or a narrow optically thick inner ring, to reproduce the fluxes at 22 and 24 $\mu {\rm m}.$

The local parameter exploration (Figure 9) represents a useful approach to estimate the sensitivity of the model parameters to our observations. We find that most of the parameters describing Regions 1 and 2 remain unconstrained from our observations. For instance, in our representative model there are no empty gaps, as we find ${{R}_{{\rm out},1}}={{R}_{{\rm in},2}}$ and ${{R}_{{\rm out},2}}={{R}_{{\rm cav}}}$, although other values for these parameters (i.e., an empty gap between Regions 1 and 2 or a dust empty innermost region) could equally reproduce our observations (as indicated by the ${\Delta }{{\chi }^{2}}$ profile in the R_in and R_out panels in Figure 9). For all three regions we find that the the exponents of the surface density distributions are completely unconstrained by our observations. The same happens with the flaring index ψ, as expected if the emission is optically thin at these wavelengths. However, the exercise we performed does not allow us to derive proper statistical errors of the fitted parameters, as we did not investigate correlations between the parameters. The latter would require us to map the entire multidimensional parameter space, which is beyond the scope of this paper and would require spending an unreasonable amount of computing time. We reserve this exercise for the future, when data with more sensitivity and much better angular resolution will be available. Regarding the gas model, we find that it is difficult to reproduce the observed ¹²CO profile with our model. Overall, we see that compared to the observations the model shows a lack of emission at higher velocities and an excess at lower velocities. Even using a non-depleted gaseous disk, our model produces less flux at higher velocities than the observations. We note that this does not exclude the possibility of a gas-depleted cavity (as found for HD 142527; Perez et al. 2015), since we also find that it is possible to produce a similar ¹²CO profile by modeling a gas-depleted cavity using higher values of turbulence or using different surface density profiles for the gaseous and dusty disks (i.e., different γ for the gas and dust disks).

With a size of 97 AU in radius, Sz 91 has by far the largest the inner cavity observed in a transition disk around a T Tauri star. This cavity size and accretion rate ($7.9\times {{10}^{-11}}\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}}$) exclude photoevaporation alone as the mechanism responsible for the cavity (e.g., Owen et al. 2012; Rosotti et al. 2013). Similarly, studies of grain growth conclude that grain growth alone cannot produce the large cavities resolved by interferometric millimeter observations (Birnstiel et al. 2013). Romero et al. (2012) excluded a stellar companion down to separations of ∼30 AU with a contrast of ${\Delta }K\sim 3.3$. Using NextGen models (e.g., Allard et al. 1997), an age of 1 Myr and the stellar mass derived in Section 4.1 translate to a sensitivity close to the brown dwarf limit of $\sim 80\;{{M}_{{\rm J}}}$. In addition, Melo (2003) found no evidence for a close-in binary companion around Sz 91 in their 3 yr radial velocity survey (down to masses of $\approx 0.2\;{{M}_{\odot }}$). However, we cannot yet exclude that a low-mass stellar companion is carving the cavity in the disk around Sz 91.The other possible mechanism is planet formation. In fact, our direct detection of gas inside of the large dust-depleted cavity, the presence of μm-sized dust inside the cavity, and the large cavity size indicated by the continuum visibilities matches several predictions of models for multiple (Jupiter-sized) planet formation (e.g., Lubow and D'Angelo 2006; Dodson-Robinson & Salyk 2011; Zhu et al. 2011).

There are at least four transition objects with detectable (but different) accretion rates that can be explained by an inner disk structure similar to the one derived for Sz 91. DM Tau and RX J1615.3-3255 have resolved cavities and require an empty disk very close to the star because of the lack of near-IR excess, but also some dust in the outer regions of the cavity to account for the MIR emission (Andrews et al. 2011). However, a direct comparison with Sz 91 is difficult because the modeling approach adopted by Andrews et al. (2011) is different from ours. RX J1633.9-2442 and [PZ99] J160421.7-213028 have been analyzed using tools similar to those used in this work for Sz 91. In both systems, an empty inner disk and some dust inside the main cavity (Cieza et al. 2012; Mathews et al. 2012; Zhang et al. 2014) are needed to explain the millimeter observations and SED. Interestingly, the millimeter-sized dust distribution in [PZ99] J160421.7-213028 is concentrated in a ring similar to what we find for Sz 91 (Mathews et al. 2012; Zhang et al. 2014). Pinilla et al. (2012) conclude that this ring-shaped dusty structure is a product of the planet–disk interaction.

While in all five cases strong parameter degeneracies of the models do not allow us to constrain the exact configuration of the dust inside the cavity, photoevaporation and grain growth can be ruled out as the mechanisms carving the inner hole of these disks. Furthermore, for DM Tau, RX J1633.9-2442, and [PZ99] close binarity can also be excluded, and planet formation has been suggested to explain both the cavities and the presence of dust inside the cavity (e.g., Rice et al. 2006; Cieza et al. 2012; Mathews et al. 2012). As suggested by Cieza et al. (2012), these inner structures are likely simple approximations to the complex structures predicted by hydrodynamical models of planet formation (Dodson-Robinson & Salyk 2011), which have now been directly observed (e.g., Casassus et al. 2013; Quanz et al. 2013) in a few systems.