Hubble Space Telescope Scattered-light Imaging and Modeling of the Edge-on Protoplanetary Disk ESO-Hα 569

ESO-Hα 569 (2MASS J11111083−7641574) was first identified as a target of interest in the Comerón et al. (2004) European Southern Observatory survey of young stars with strong Hα emission in Cha I SFR. Comerón et al. (2004) classified the central star as K7 using ground-based spectroscopy. The authors noted that this object is severely underluminous for a K7 star (by ∼2 orders of magnitude), which made it a prime candidate for our edge-on disks survey. The Luhman (2007) survey of the stellar population in Chamaeleon obtained an $R\approx 5000$ spectrum from 0.6 to 0.9 $\mu {\rm{m}}$, which gave a spectral type of M2.5, an effective temperature of 3488 K, and an apparent bolometric luminosity of L_bol = 0.0030 L_⊙. More recently, broadband spectroscopy with VLT/X-Shooter provides a spectral type of M1 ± 2 subtypes (Manara et al. 2017, their Table 3), and confirms that the target appears underluminous. Because ESO-Hα 569 is heavily extincted by the disk, the apparent luminosity is an unreliable estimator for the true bolometric luminosity of the central star. For stars of the same spectral type in the Luhman (2007) survey, the average bolometric luminosity is $0.34\pm 0.08{L}_{\odot }$. The current study adopts this luminosity and a ∼3500 K effective temperature. Using the theoretical evolutionary models of Baraffe et al. (1998) for low-mass stars with solar metallicity gives a mass for the central star of 0.35 M_⊙. The associated stellar radius is 1.13 R_⊙.

Prior attempts have been made to infer the disk properties of ESO-Hα 569 based on its SED. Luhman (2007) noted that the X-ray non-detection of this star indicates an extinction of ${A}_{K}\geqslant 60$, consistent with obscuration by an edge-on disk, assuming its X-ray luminosity is that of a typical T Tauri star. Robberto et al. (2012) combine published 2MASS and Spitzer photometry, with unresolved HST fluxes to fit properties of the disk and central star using the online library of 20,000 models of young circumstellar systems compiled by Robitaille et al. (2006). These models include the central star, a diffuse envelope, and an accreting disk (Whitney et al. 2003a, 2003b, 2004). The authors find that the disk is best fit by an inclination of ∼871, ${L}_{\mathrm{bol}}=0.8\pm 0.4{L}_{\odot }$, ${M}_{\mathrm{star}}=0.33\pm 0.03{M}_{\odot }$, and ${R}_{\mathrm{star}}=2.5\pm 0.6{R}_{\odot }$, and give an upper limit for the submillimeter disk mass of $0.005{M}_{\odot }$. Rodgers-Lee et al. (2014) included Herschel data and found a best-fit inclination of 814. More recently, Pascucci et al. (2016) provided 1.3 mm continuum data, which correspond to a disk mass estimate of $0.0046\,{M}_{\odot }$, using an assumed opacity of $\kappa =2.3\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, a gas-to-dust ratio of 100, and a disk temperature of 20 K.

Scattered-light images of the ESO-Hα 569 disk were obtained using the HST Advanced Camera for Surveys (ACS)/Wide-field Camera (WFC) in both the F814W and F606W broadband filters on 2012 March 9 as part of program GO 12514. The total exposure times were 1440 s for F606W and 960 s for F814, with each filter's exposure split as two integrations for cosmic-ray rejection. The reduced and calibrated data produced by the HST pipeline were retrieved from the Mikulski Archive for Space Telescopes (MAST).

Figure 1 provides the reduced images, rotated to place the disk major axis horizontally. The bipolar appearance unequivocally demonstrates the edge-on nature of ESO-Hα 569. The western side is much brighter than the eastern side (by $\sim 20\times$ comparing their peak surface brightnesses) and, along with the curvature of the nebula, indicates that this side is tilted slightly toward us. There is no sign of starlight directly peeking through as an unresolved point source. The position angle of the disk's minor axis was evaluated to be 65° ± 1°. This was computed as the position angle for which mirroring the image across the minor axis minimized the flux difference between the left and right sides. The disk is close to left/right symmetric, although the southern side (right side as shown in Figure 1) is very slightly brighter.

Zoom In Zoom Out Reset image size

The disk is very red (much brighter in F814W than in F606W). The flux density of the disk was measured in both filters using a 50 pixel aperture, which corresponds to a spatial scale of 2'' × 2'' and was chosen to encompass all disk flux with surface brightness $\geqslant 3\sigma$ above the background noise. The measured flux density is 0.058 ± 0.001 mJy for F606W and 0.21 ± 0.01 mJy for F814W, which gives a color [F814W] – [F606W] = 1.4 AB magnitudes.

The disk has an apparent outer radius of 0.80 ± 005 which corresponds to 125 ± 8 au at a distance of 160 pc. Here, the outer radius is inferred as the offset at which the flux declines to less than 10% of the peak value for the widest part of the disk.

The strong Hα emission in the spectrum of this young object indicates ongoing accretion onto the central star, which is often associated with the launching of outflow jets. Bally et al. (2006) suggested that ESO-Hα 569 as the possible source for the Herbig Haro object 919. HH 919 is an arcminute-long filament with a PA of ∼60°–75° and is located $22^{\prime\prime}$ (0.05 pc) southwest of ESO-Hα 569. A jet is visible in the F606W scattered-light image extending vertically from the disk and is ∼$0\buildrel{\prime\prime}\over{.} 25$ wide. This is consistent with the emission lines of Hα and S ii as are commonly seen in such outflows. A line connecting ESO-Hα 569 with HH 919 would have a PA of ∼63^◦, giving an orientation consistent with the ESO-Hα 569 jet serving as the culprit for the HH 919 filament.

Figure 2 presents a wider field of view showing the interaction of this disk with the surrounding ISM. Diffuse nebulosity is visible extending outward from the disk. An image of the jet was created by subtracting the F814W image (scaled by a factor of 2.5) from the F606W image (Figure 3). The flux from the jet is difficult to decouple from the disk flux, but the jet accounts for roughly 50% of the local surface brightness from the disk. This value is taken from an average of the flux over 9 pixels with the jet superimposed on the disk and compared to the flux in nine neighboring pixels with no jet signature. The peak surface brightness of the jet is ∼0.19 mJy arcsec⁻². The ability to measure color variations in the shape of the disk and width of the dark lane between the F606W and F814W bands is hindered by the presence of this bright jet.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An SED for the disk was compiled from the literature, including data from HST, 2MASS, Spitzer, WISE, Herschel, ALMA, and the LABOCA instrument on the APEX telescope (see Figure 4). Table 1 provides the SED values with photometric errors and references for each value. The SED shows the characteristic double-peaked shape of edge-on disks with contributions from both the scattered light from the central star peaking at about 1.5 μm and the thermal emission from the surrounding optically thick disk peaking at roughly 70 μm. Data at similar wavelengths from different epochs show variability at the 10%–20% level, consistent with the variability seen in other young disks (Muzerolle et al. 2009; Espaillat et al. 2011; Flaherty et al. 2012).

Zoom In Zoom Out Reset image size

Table 1.  Spectral Energy Distribution Photometry and References

$\lambda (\mu {\rm{m}})$	Flux (mJy)	Source	Instrument	Bandwidth ($\mu {\rm{m}}$)	Angular Resolution	Date
0.551	0.030 ± 0.004	Robberto et al. (2012)	HST WFPC2	0.14	$0\buildrel{\prime\prime}\over{.} 0996$	2009 Apr 27
0.606	0.058 ± 0.001	This work	HST ACS	0.27	$0\buildrel{\prime\prime}\over{.} 05$	2012 Mar 09
0.814	0.21 ± 0.01	This work	HST ACS	0.31	$0\buildrel{\prime\prime}\over{.} 05$	2012 Mar 09
1.235	0.66 ± 0.05	Skrutskie et al. (2006)	2MASS	0.16	∼5''	2000 Jan 25
1.662	0.97 ± 0.08	Skrutskie et al. (2006)	2MASS	0.25	∼5''	2000 Jan 25
2.15	0.98 ± 0.09	Skrutskie et al. (2006)	2MASS	0.26	∼5''	2000 Jan 25
3.6	0.58 ± 0.03	Luhman et al. (2008)	Spitzer IRAC	0.75	∼2''	2004 Jul 04
4.5	0.57 ± 0.05	Luhman et al. (2008)	Spitzer IRAC	1.02	∼2''	2004 Jul 04
5.8	0.58 ± 0.05	Luhman et al. (2008)	Spitzer IRAC	1.43	∼2''	2004 Jul 04
8.0	0.67 ± 0.05	Luhman et al. (2008)	Spitzer IRAC	2.91	∼2''	2004 Jul 04
3.4	0.63 ± 0.02	Cutri et al. (2012)	WISE	0.66	$6\buildrel{\prime\prime}\over{.} 1$	2010 Feb 13, 20
4.6	0.71 ± 0.02	Cutri et al. (2012)	WISE	1.04	$6\buildrel{\prime\prime}\over{.} 4$	2010 Feb 13, 20
12	0.65 ± 0.07	Cutri et al. (2012)	WISE	5.51	$6\buildrel{\prime\prime}\over{.} 5$	2010 Feb 13, 20
22	7.5 ± 0.89	Cutri et al. (2012)	WISE	4.10	$12\buildrel{\prime\prime}\over{.} 0$	2010 Feb 13, 20
24	8.36 ± 0.77	Luhman et al. (2008)	Spitzer MIPS	5.3	$6^{\prime\prime}$	2004 Apr 08
70	107 ± 10.8	Luhman et al. (2008)	Spitzer MIPS	19	$18^{\prime\prime}$	2004 Apr 08
70	200 ± 100	Winston et al. (2012)	Herschel PACS	25	$5\buildrel{\prime\prime}\over{.} 8$	2011 Jun 23
160*	200 ± 200	Winston et al. (2012)	Herschel PACS	85	$12\buildrel{\prime\prime}\over{.} 0$	2011 Jun 23
250*	100 ± 100	Winston et al. (2012)	Herschel SPIRE	25	$18^{\prime\prime}$	2011 Jun 23
350*	50 ± 50	Winston et al. (2012)	Herschel SPIRE	25	$25^{\prime\prime}$	2011 Jun 23
500*	50 ± 50	Winston et al. (2012)	Herschel SPIRE	25	$37^{\prime\prime}$	2011 Jun 23
870	72 ± 14	Belloche et al. (2011)	APEX/LABOCA	150	$19\buildrel{\prime\prime}\over{.} 2$	2008 May
2830	3.2 ± 0.1	Dunham et al. (2016)	ALMA	55	∼2''	2013 Nov 29–2014 Mar 08

Note. Photometry at wavelengths marked with an "*" represent only upper limits and are not included in the spectral energy distribution modeling.

Download table as:  ASCIITypeset image

ESO-Hα 569 was imaged with Herschel as part of the Gould Belt survey in the PACS 70 and 160 μm bands and the SPIRE in the 250, 350, and 500 μm bands (Winston et al. 2012). The source is barely detected in the PACS bands, hence the large uncertainties reported by Winston et al. (2012). There seems to be a very marginally detected point source in the SPIRE bands (100 ± 100 mJy at 250 μm and 50 ± 50 mJy at 350 and 500 μm), but given the coarse angular resolution, it is hard to exclude contamination from dust emission from the surrounding cloud itself. Given the low significance of these detections, the Herschel fluxes are not included in our SED fits.

An ALMA Band 3 continuum (2.8 mm, 106 GHz) measurement was obtained by Dunham et al. (2016) with ALMA in a compact configuration that achieved a ∼2'' beam and does not resolve the disk. Because the disk is not resolved, the ALMA continuum flux could be contaminated with flux from a remnant envelope. However, there cannot be too much non-disk material present, or it would be too opaque to see the central disk in scattered light in the visible as we do. This measurement was published after our initial rounds of disk SED fitting as described below (the grid fit described in Sections 3.2–3.3, and the ${\chi }^{2}$-based MCMC fit in Section 3.4.1), but this data point has been included in our final SED model fitting (the covariance-based MCMC fit described in Section 3.4.2).

In addition to the continuum measurement at 0.55 μm included in Table 1, Robberto et al. (2012) provide fluxes for ESO-Hα 569 in HST WFPC2's F631N, F656N, and F673N narrowband filters associated with [O i], Hα, and [S ii] emission, respectively. The disk is not resolved and the measured fluxes are near the detection limits: 0.21 ± 0.15 (F631N), 2.4 ± 0.7 (F656N), and 0.35 ± 0.12 (F673N) $\times \,{10}^{-16}$ erg s⁻¹ cm⁻² Å⁻¹. These emission lines are all consistent with the spectrum of Luhman (2007), which shows strong Hα emission and [S ii] emission. Given the large uncertainties, and the fact that the model is not set up to simulate line emission, these emission lines were not included in the SED fits.

Circumstellar disks are complex objects: mixtures of gas and dust, containing solid bodies from the smallest planetisimals to giant Jovian planets, shaped by many dynamical forces across evolutionary states from the youngest protoplanetary disks through transitional regimes to second-generation debris disks. This complexity can now be probed by powerful observational capabilities across the entire electromagnetic spectrum, with especially detailed views provided in the visible by the HST, in the infrared by 8–10 m telescopes with adaptive optics and soon by JWST, and in the millimeter and submillimeter by ALMA and other interferometers. In some cases, a particular physical property of interest can be directly measured from a given observation, but more typically forward modeling of the data must be performed to derive constraints on the underlying physics. This is particularly necessary for observations of disks at wavelengths where they are optically thick, which is the case for observations of protoplanetary disks at visual and near-IR wavelengths.

The general outline of such inference is well known: start from a model of the system's properties and physics with some number of free parameters. Construct synthetic observables using that model, for instance through Monte Carlo radiative transfer (MCRT) calculations. Then compare the synthetic observables to data in order to constrain the free parameters and draw conclusions about their most likely values and the ranges of uncertainty. This process sounds simple enough in theory, but is often a practical challenge due to several confounding factors, among them the complexity of the underlying physics (which inevitably requires simplifications in the models), the nonetheless high dimensionality of the model parameter space, and the need to confront heterogeneous and multiwavelength observations in order to resolve model degeneracies.

The current work makes use of the MCFOST radiative transfer code (Pinte et al. 2006), one of a broad class of class of MCRT programs designed to study circumstellar disks (for a review of radiative transfer codes, see Steinacker et al. 2013). In short, such a code begins with a numerical model of the physical properties within the disk, such as the density of dust in each grid cell, and the mineralogical composition and size distribution of dust particles. It then computes the temperature and scattering source function everywhere in the disk via a Monte Carlo method: photon packets are propagated stochastically through the model volume following the equations of radiative transfer, and information on their properties is retained along their path. The radiation field and the quantities derived from it (for instance, temperature, radiation pressure, etc.) are obtained by averaging this Monte Carlo information. Observable quantities (SEDs and images) are then obtained via a ray-tracing method, which calculates the output intensities by formally integrating the source function estimated by the Monte Carlo calculations. This approach naturally allows simulation of disk images, which are dominated by scattered starlight, thermal emission from the dust, or a combination thereof.

Comparison of the simulated images and SEDs against observations then allows inference about which ranges of model parameters are compatible. There are a couple different approaches to performing such comparisons. One option is to compute a grid of models spanning the parameter space of interest (e.g., Robitaille et al. 2006; Pinte et al 2008; Woitke et al. 2010). Bayesian techniques allow the derivation of uncertainty ranges around the best-fit grid point (e.g., Chiang et al. 2012). However, even with hundreds of thousands of models computed, given the high dimensionality of the parameter spaces, each parameter must often be quite coarsely sampled at just a few discrete values, which can limit the results achieved. The grid technique is also computationally inefficient because it blindly allocates equal effort to both the best- and worst-fitting portions of parameter space. As is well known, the MCMC paradigm improves on this; the MCMC algorithm allows the efficient exploration of parameter space and yields detailed information on parameter posterior probability distributions and correlations.

However, most disk model-fitting efforts to date have concentrated on fitting either SEDs alone (e.g., Huélamo et al. 2010; Ribas et al. 2016) or images or interferometric visibilities alone (e.g., Millar-Blanchaer et al. 2015; Ricci et al. 2015; Pohl et al. 2017). This is broadly the case independent of the choice of either grid fitting or MCMC fitting. But fits to SEDs alone are notoriously degenerate (Chiang et al. 2001; Woitke 2015), and spatially resolved image data or interferometric visibilities are required in order to place robust constraints on many properties of interest. Only a handful of disk studies have successfully and rigorously fit models to heterogeneous observables, including SEDs and images or interferometric visibilities, but when this has been achieved it has often yielded particularly powerful constraints and detailed insights into disk structures (e.g., Pinte et al 2008; Duchêne et al. 2010; Lebreton et al. 2012; Carmona et al. 2014; Milli et al. 2015; Cleeves et al. 2016).

Such works have most often used the grid-fitting approach rather than MCMC, perhaps due to the increased technical complexity of integrating the MCMC framework with heterogeneous observables. One detail—but an important one in this context—is that the MCMC approach necessarily assesses a single goodness-of-fit metric which must combine both SED and image data together, such as a sum of ${\chi }^{2}$ values from the SED and image (or more generally from any combination of distinct observables). In the case where the best-fitting ${\chi }^{2}$ for one observable is systematically much higher than that for the other observable(s), the model fitting will be driven by that first observable and will likely not deliver an adequate simultaneous fit to the others. Models must necessarily simplify, and imperfect models lead to correlated systematic residuals that increase the minimum ${\chi }^{2}$. Consider, for instance, attempting to fit a simple axisymmetric model to an eccentric disk. This problem is generally worse for images than for SEDs, because the one-dimensional nature of SEDs collapses much of the parameter space. In other words, the well-known degeneracies of SEDs can hide disk offsets, eccentricities, spiral arms, and other asymmetries that are immediately apparent in sufficiently high-resolution images. As a result, it becomes difficult to develop a good metric that combines both images and SEDs in a well-balanced manner for the purposes of a simultaneous MCMC fit.

To address this difficulty in fitting disk observations, a new method has been developed that explicitly takes into account the covariant and correlated residuals in the image fitting. Czekala et al. (2015) introduced this approach in the context of 1D spectral fitting. That approach has been extended to work on heterogeneous disk data sets, including two-dimensional images, and use that to implement an MCMC fitting process that balances both the image and SED data for ESO-Hα 569.

For this work, the MCFOST radiative transfer code (Pinte et al. 2006, 2009) was used to construct SEDs and 0.8 $\mu {\rm{m}}$ scattered-light images for each of the models. The 0.6 $\mu {\rm{m}}$ scattered-light images were not modelled because the strong jet signature required masking ≥50% of the integrated disk flux.

The selected model assumes an axisymmetric disk with a surface density, Σ, described by a power-law distribution in radius given by ${\rm{\Sigma }}={{\rm{\Sigma }}}_{0}{(R/{R}_{0})}^{\alpha }$, where α is termed the surface density exponent and R₀ is the reference radius of 100 au. In this "sharp-edged" model, the disk is abruptly truncated at an outer radius R_out. In order to achieve a good fit to the diffuse emission above the disk and the disk mass and inclination simultaneously, a "tapered-edged" disk model was tested, in which the density Σ falls off exponentially with some critical radius R_c of material outside of the disk:¹⁰

$$\begin{eqnarray}&&{\rm{\Sigma }}={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{\alpha }\exp \left[{\left(-\displaystyle \frac{R}{{R}_{c}}\right)}^{2+\alpha }\right].\end{eqnarray} \tag{ 1 }$$

For this work, ${R}_{c}={R}_{\mathrm{out}}$. This exponential taper is predicted by physical models of viscous accretion disks (Hartmann et al. 1998), but observations were not sensitive enough to detect this outer gradual falloff until Hughes et al. (2008) used this form to model both gas and dust continuum observations in the millimeter. It is expected that the small dust grains seen in scattered light should be well-coupled with the gas for young disks, suggesting that the use of this surface density distribution is justified here. (See also the recent work by Guidi et al. 2016 and Pohl et al. 2017 for HD 163296 and T Cha, respectively). The scale height is also defined as a power law in radius by $H{(R)={H}_{0}(R/{R}_{0})}^{\beta }$, where β is the flaring exponent describing the curvature of the disk and again R₀ = 100 au.

Several model parameters were held fixed to minimize the degrees of freedom and to save computation time. Values for these parameters were either measured directly from the HST images or taken from the literature. The disk is within the SFR Cha I; therefore, we fix the distance to the disk at 160 pc (Whittet et al. 1997). From the angular size of the disk measured above and the distance, we calculate an outer radius of 125 au. The inner radius was defined by a conservative estimate of the sublimation radius, ${R}_{\mathrm{sub}}={R}_{\mathrm{star}}{({T}_{\mathrm{star}}/{T}_{\mathrm{sub}})}^{2.1}\sim 0.1\,\mathrm{au}$, where ${T}_{\mathrm{sub}}=1600\,{\rm{K}}$ (Robitaille et al. 2006).

The free parameters in the model are inclination (with 90° as edge-on), scale height, dust mass, maximum dust particle size, dust porosity, disk vertical flaring exponent (β), surface density exponent (α), and disk edge type (sharp or tapered). For the maximum particle size, the grain population is described by a single species of amorphous dust of olivine composition (Dorschner et al. 1995) with a particle size distribution following a −3.5 power law extending from 0.03 μm up to the free parameter a_max. We assume that the dust is well-mixed with the gas, irrespective of the particle size. This combination of dust properties (with ${a}_{\max }=100\,\mu {\rm{m}}$) results in a mean scattering phase function asymmetry factor of g = 0.54. Dust porosity is modeled simply as a fraction between 0 and 1 of vacuum that is mixed with the silicates following the Bruggeman effective mixing rule.

For comparison with the observed 0.8 $\mu {\rm{m}}$ scattered-light images, each model image was convolved with a Tiny Tim simulated PSF (Krist 1995). The 0.8 $\mu {\rm{m}}$ observations were masked to select only the pixels with flux values $\geqslant 3\sigma$ above the background noise level. A 2D map of the noise was generated by converting the observed image to electrons and assigning $\sigma =\sqrt{{N}_{e-}}$ for the ${\chi }^{2}$ values. The model images were aligned with the observations via a cross-correlation and normalized to the total observed flux. The models were then compared to the data via an error-weighted pixel-by-pixel ${\chi }^{2}$ calculation. For a similar work, see Duchêne et al. (2010) and McCabe et al. (2011). For the SEDs, when fitting each model point, the foreground extinction is allowed to vary from ${A}_{V}=0\,-\,10$ with ${R}_{V}=3.1$, and the extinction value that minimizes the observed–model residuals is chosen.

Although the robust treatment of radiative transfer provided by MCFOST is essential for modeling optically thick disks, it is computationally intensive. Generating a single model SED requires ∼three minutes of desktop CPU time, with an additional ∼minute to generate synthetic images at each desired wavelength. MCFOST allows the user to parallelize the computation; however, systematic explorations of parameter space can quickly become very time consuming.

This complex parameter space was explored in two stages using two different techniques. First, a coarse model grid was computed with a wide range of allowed model parameter values to get a handle on reasonable regions of parameter space. Section 3.3 describes the initial exploration of parameter space via a grid search, with results in Section 3.4. This work was used to inform a more robust Markov Chain Monte Carlo exploration for finer sampling of allowed parameter values, with the methods described in Section 3.5 and the results given in Section 3.6.

Our initial modeling used a uniform grid sampling, with the explored parameter space shown in Table 2. For each set of disk model parameters, 15 disk inclinations were sampled uniformly in $\cos i$ between 60° and 90°. This resulted in a grid of over 200,000 models. Comparison with data was performed using a custom IDL software. A benefit of the grid search approach is that multiple goodness-of-fit metrics may be evaluated across all sampled points. ${\chi }^{2}$ values were computed separately for the 0.8 $\mu {\rm{m}}$ image and SED for each model along with the combined total ${\chi }_{\mathrm{tot}}^{2}={\chi }_{0.8\,\mu {\rm{m}}}^{2}+{\chi }_{\mathrm{SED}}^{2}$. Bayesian probabilities are derived from the likelihood function wherein the ${\chi }^{2}$ value for a given model with unique parameter values is related to a probability $\exp (-{\chi }^{2}/2)$, and the sum of all probabilities is normalized to unity (e.g., Pinte et al 2008).

Table 2.  Modeled Disk Parameters

Parameter	Grid Values	MCMC Values
Distance (pc)	160 (Fixed)	160 (Fixed)
Outer Radius (au)	125 (Fixed)	125 (Fixed)
Minimum Particle Size ($\mu {\rm{m}}$)	0.03 (Fixed)	0.03 (Fixed)
Inclination (°)	60–90	65–90
Scale Height (H in au at R = 100 au)	10, 15, 20, 25	5–25
Dust Mass (M in M⊙)	10−4, $3\times {10}^{-4}$, 10−3	10−5–10−3
Surface Density (α)	−2.0, −1.5, −1.0, −0.5, 0.0	−2.0 to 0
Flaring Exponent (β)	1.1, 1.2, 1.3, 1.4, 1.5	1.0–1.5
Maximum Grain Sizea ($\mu {\rm{m}}$)	100, 1000, 3000	100–3000
Weightb	⋯	0.3–0.7
Grain Porosity	0.0, 0.25, 0.5	⋯
Structure	Disk, Tapered-edge Disk	⋯

Notes.

^aGrain size was kept at a constant value of 100 $\mu {\rm{m}}$ for the covariance-based MCMC run. ^bDuring the ${\chi }^{2}$-based MCMC run, a weighting term was used to describe the relative contribution of the image and SED fits to the log-likelihood value of each model.

Download table as:  ASCIITypeset image

The grid sampling is a simple way to explore the parameter space initially, but its sampling of parameters proved to be inadequate for several reasons. First, it is too sparse to provide clear insight into degeneracies between the various parameters.Second, the discrete sampling limits the precision with which best-fit values can be determined, and does not allow a rigorous computation of uncertainties. These factors motivated the later development of our MCMC model-fitting toolkit described below. Nonetheless, the results of the grid search helped clarify relevant portions of parameter space and informed our understanding of the disk.

For the grid search approach, the best-fit model for the disk was found using a tapered-edged disk with non-porous grains, an inclination of 755, and a scale height of 20 au at a reference radius of 100 au. The preferred maximum particle size is 3000 $\mu {\rm{m}}$, the dust mass is $3\,\times {10}^{-4}$ M_⊙, the flaring exponent β is 1.3, and the surface density exponent α is −0.5. The separate SED and image fits for the α and β exponents favor opposing extremes of parameter space, but the combined ${\chi }_{\mathrm{tot}}^{2}$ likelihood distribution peaks in the middle at physically reasonable values.

Figure 5 illustrates the likelihood distributions for the inclination and scale height. The sparse sampling and disagreement between the model parameters preferred by the image and SED (most pronounced in the scale height) demonstrate the limitations of the grid-fitting approach.

Zoom In Zoom Out Reset image size

3.4.1. Porosity

3.4.2. Disk Structure: Sharp versus Tapered Outer Edge

When modeling the disk with a sharp outer edge, the SED and image fits preferred very different regions of parameter space. Specifically, it was difficult to simultaneously fit the flux ratio between the top and bottom nebulae of the disk, the diffuse emission above the plane of the disk, and the shape of the disk. Because the disk is not precisely edge-on, the scattering angles differ between the upper and lower disk nebulae. Therefore, changes in the scattering phase function of the grains will change the peak-to-peak flux ratio. Any parameter that would increase the flux ratio and emission above the disk (for example increasing the inclination or porosity of the grains) caused too much forward-scattering and allowed too much of the light from the central star to appear in the peak. Similarly, the diffuse emission above the disk could not be described well by a low-mass spherical envelope.

The tapered-edged disk did much better in accounting for both the emission above the disk and matching the flux ratio between the top and bottom sides of the disk. This is clearly demonstrated in Figure 6, which compares the observations to the best-fit tapered-edged disk model and the corresponding sharp-edged model. The right panel shows the surface brightness profiles through several vertical cuts across the disk for both the sharp- and tapered-edged models.

Zoom In Zoom Out Reset image size

To more efficiently the sample parameter space and gain a better understanding of the uncertainties, an MCMC approach was applied to the model optimization. We used the Python package emcee (Foreman-Mackey et al. 2013) which implements the Affine Invariant MCMC algorithm by Goodman & Weare (2010). Specifically, we selected the parallel-tempered MCMC ensemble sampler designed to improve convergence in degenerate parameter spaces. The MCMC samples the posterior distribution given by

$$\begin{eqnarray}&&P({\rm{\Theta }}| D)\propto P(D| {\rm{\Theta }})P({\rm{\Theta }}),\end{eqnarray} \tag{ 2 }$$

where D represents the observations and Θ the free parameters in the model. Here, P(D $| {\rm{\Theta }}$) is the likelihood of the data given the model and P(Θ) is the prior distribution. Uniform priors were adopted for each parameter over the allowable range.

In order to implement this code in conjunction with the MCFOST radiative transfer code, we developed a suite of software tools in Python to interact with the observations, generate models, and calculate goodness-of-fit metrics to inform the MCMC iterations. The toolkit is general enough to be usable with any disk image, provided a PSF and uncertainty map are available. By combining the detailed modeling capabilities of MCFOST with the efficient parameter space sampling of the emcee package, the goal was to self-consistently and simultaneously fit a wide variety of observables in order to place constraints on the physical properties of a given disk, while also rigorously assessing the uncertainties in the derived properties. The mcfost-python package is publicly available on GitHub,¹¹ and the authors encourage its use by the disk modeling community (Wolff et al. 2017).

The mcfost-python package was designed to be modular, with different components to read in the observables, interact with the MCFOST parameter files, generate model SEDs and images, compare them to data, and set up and control the overall MCMC run. To validate the functions for comparing models to data, benchmark cross-checks were performed to compare the new Python fitting code to existing ${\chi }^{2}$ routines in IDL and Yorick. Although this code was originally designed to work with HST data and the MCFOST modeling package as described in this paper, it has also been expanded to work with data from different instruments, including polarimetry data, and can be used with other radiative transfer modeling codes.

3.5.1.  ${\chi }^{2}$-based Log-likelihood Estimation

The mcfost-python package allows the user to choose between two goodness-of-fit metrics. This section discusses the first of those, the ${\chi }^{2}$ metric. A simple benchmark comparison of the ${\chi }^{2}$ and covariance likelihood methods is provided in the Appendix. At each step in the MCMC iteration, a model image and SED are created for the chosen parameter values and a ${\chi }^{2}$ value is calculated using the same methodology as the grid sampling approach. The emcee code requires a log-likelihood distribution which is computed from the ${\chi }^{2}$ assuming a multidimensional Gaussian likelihood function:

$$\begin{eqnarray}&&\mathrm{ln}[P(D| {\rm{\Theta }})]=-0.5N\mathrm{ln}2\pi +\displaystyle \sum _{i=1}^{N}(-\mathrm{ln}| {\sigma }_{i}| )-\displaystyle \frac{1}{2}{\chi }_{i}^{2}.\end{eqnarray} \tag{ 3 }$$

Here, N is the number of data points and σ is our uncertainty. The MCMC approach inherently requires a single goodness-of-fit metric, so it is essential to combine the SED and image metrics into a single log-likelihood function for use by emcee. The log-likelihood distribution is computed separately for the images and SEDs, and a weighted average is used to determine the goodness of fit. During initial tests using the ${\chi }^{2}$-based log-likelihood goodness-of-fit metric, we chose to allow the relative weighting between the image and SED to vary. The best way to handle the relative weighting between different types of observations for a single disk model was not well-understood and is a nuisance parameter that does not, itself, inform us about any inherent physical properties of the disk. By marginalizing over it in this way, the intent was to produce a best-fit model that was informed by both the SED and image data without a bias toward one or the other. The weighting was allowed to vary between 0.3 and 0.7 for a minimum of 30% weighting to either the image or SED fits. We found that the image likelihood values were downweighted due to their systematically higher ${\chi }^{2}$ values, and the MCMC chains worked to improve the images while largely ignoring the better SED fits. In our first round of MCMC calculations, the image reduced ${\chi }^{2}$ values tended to be more than an order of magnitude above the SED reduced ${\chi }^{2}$ values (best ${\chi }_{\mathrm{SED}}^{2}=1.3$, ${\chi }_{0.8\,\mu {\rm{m}}}^{2}=66$), due to the larger number of measurements in the images presumably with underestimated uncertainties.

3.5.2. Covariance-based Log-likelihood Estimation

The imbalance between the image and SED ${\chi }^{2}$ values served as the impetus for the development of the covariance matrix likelihood estimation method, which ultimately provided a much better relative weighting of the different observables. Given that each model image is convolved with an instrumental PSF, neighboring pixels must be covariant. Furthermore, this approach lets us correct for the global limitations of the disk model to fit the data set. Model systematics present as correlated uncertainties. For a more complete estimate of the errors in our HST images, we adopt and extend the covariance-based method for the log-likelihood estimation presented by Czekala et al. (2015) in the context of 1D spectral fitting. That approach must be extended to work in the context of 2D images. In this case, we convert Equation (3), which describes the likelihood of the data given the model assuming a Gaussian likelihood distribution, into a matrix formalism in Equation (4):

$$\begin{eqnarray}\mathrm{ln}[P(D| {\rm{\Theta }})] & = & -\displaystyle \frac{1}{2}({R}^{T}{C}^{-1}R+\mathrm{ln}[\det (C)]\\ & & +N\mathrm{ln}[2\pi ]),\end{eqnarray} \tag{ 4 }$$

where R represents the residuals of the observations subtracted by the model, C is the covariance matrix defined below, and N is the total number of pixels in the image (not the number of pixels along a given dimension of the array).

To apply this approach, each 2D image must first be "unwrapped" into a 1D array. In practice, not all pixels in an image may have a sufficient signal-to-noise ratio disk detection to justify fitting. Excluding such pixels from the unwrapping improves the overall computational efficiency, particularly for the matrix inversion calculation, at the cost of somewhat more complex bookkeeping between the 2D and 1D versions of the image.

The covariance matrix C (of size ${N}_{\mathrm{pix}}\times {N}_{\mathrm{pix}}$) incorporates both the noise in each individual pixel and global covariances between adjacent pixels (represented by K^G): ${C}_{i,j}\,={\delta }_{i,j}{\sigma }_{i,j}^{2}+{K}_{i,j}^{G}$. An example source of global covariance is the FWHM of a telescope PSF. For a non-zero PSF FWHM, neighboring pixels cannot be treated as individual measurements of the disk surface brightness. Additionally, any global limitations of the model to fit the data can be implicitly included in the covariance structure. For example, when using a symmetric disk model, any asymmetries in the observed image of the disk will necessarily lead to higher correlated residuals even for the best-fitting model parameters. These residuals will in general be spatially correlated on one or more scales from the angular resolution to the size of the observed asymmetry. Incorporating our knowledge of these residuals in the covariance matrix improves our ability to draw conclusions given such necessarily imperfect models. Likewise, the choice of incomplete or simplified parameterizations of the disk physics/structure in the model can be handled the same way. For instance, if there exists an additional unmodelled component, such as a more vertically extended disk atmosphere or significant residual jet emission on the top/bottom on the disk, or if the functional form of the power law adopted for the disk surface density is an oversimplified description of the true disk properties, such systematics would lead to correlated residuals in data–model comparisons. This covariance framework allows the downweighting of these contributions within the correlated residuals without masking them altogether.

The field of Gaussian processes has developed several useful analytic models for convolution kernels that can be used to parameterize the covariant structure. For instance, Czekala et al. (2015) adopt the Matérn kernel truncated by a Hann window function. This kernel has several free parameters, which can be solved for as nuisance parameters as part of the MCMC fit. Of course, this increases the dimensionality of the parameter space that must be explored, which can in practice increase computation time by an order of magnitude or more. Czekala et al. (2015) note that, because the best-fit model parameters are relatively insensitive to the precise values of the covariance parameters (i.e., a reasonably good but perhaps not optimal covariance model often suffices), one can first roughly optimize the covariance model and then perform the MCMC fit with that model fixed. Given the computational demands of disk radiative transfer model fitting, a variant of that approach is adopted here.

The global covariance is estimated empirically by computing the average autocorrelation of the residuals from a subtraction of our 0.8 $\mu {\rm{m}}$ image and a subset of 1000 randomly chosen model disk images from a uniform sampling of the parameter space within the limits of our priors (Figure 7). This provides, in a computationally tractable way, a reasonable model for the covariant structure found in residuals for the parameter space of interest and allows us to hold the covariance model fixed in subsequent MCMC runs. The 2D autocorrelation is collapsed along the horizontal axis to generate a 1D autocorrelation function (Figure 8). The horizontal axis was chosen to generate the covariance matrix because it provided the most conservative estimate, with a wider tail similar to the Matérn kernel, and did not exhibit the anticorrelation found in the vertical axis due to the dark lane of the disk. For comparison, several $\nu =3/2$ Matérn kernels are shown, following the chosen formalism from Czekala et al. (2015). To compute the covariance matrix K^G, for each pair of pixels $i,j$, the distance between them given by ${r}_{i,j}=\sqrt{{({x}_{i}-{x}_{j})}^{2}+{({y}_{i}-{y}_{j})}^{2}}$ is calculated. For each entry of ${K}_{i,j}^{G}$, the analytic autocorrelation function is interpolated to the value for ${r}_{i,j}$, with a cutoff outside of 20 pixels to make computations of ${C}_{i,j}$ manageable. The resulting covariance matrix is shown in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For consistency, the likelihood of each model SED is computed using the covariance matrix framework from Equation (4). In this case, the covariance matrix contains only the individual uncertainties for each point multiplied by an identity matrix. Any global limitations of the model SEDs to fit the data set are neglected. Given the low ${\chi }^{2}$ values achieved for the SED fitting in the grid search described above (lowest SED ${\chi }^{2}$ ∼1.3), the uncertainties in the SED are presumed to be well-estimated. For this data set, covariances are not anticipated between neighboring photometric points, but this formalism would naturally handle any such correlations and makes it straightforward to include continuous spectra as part of a unified fit alongside broadband photometry.

The covariance framework is also capable of including model terms for additional regions of the locally covariant structure (K^L), as discussed in Czekala et al. (2015). We leave the application of such local covariances to disk image fitting for future work, along with the exploration of how best to explicitly model covariances between the SED and image portions of the overall fit.

3.5.3. Choice of Parameter Values for MCMC

The best-fit parameter values are shown in Table 3. The data are best fit by a tapered-edged disk with an inclination of ${83.0}_{-4.8}^{+2.6}$ degrees, a scale height of ${16.2}_{-2.0}^{+1.7}$ au at a reference radius of 100 au, a disk dust mass of ${0.00057}_{-0.00022}^{+0.00017}$ M_⊙ (assuming a gas-to-dust ratio of 100), a surface density exponent (α) of $-{1.77}_{-0.14}^{+0.94}$, and a flaring exponent (β) of ${1.19}_{-0.08}^{+0.09}$. The image and SED combined best-fit model is illustrated in Figure 10 (single best fit in red, along with an ensemble of well-fitting models in gray) and together provide a close fit to the observations. Parameter distributions are shown in Figure 11. The results for individual parameters are discussed in more detail below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  MCMC Best-fit Parameters

Parameters	Best-fit Values
Inclination (°)	${83.0}_{-4.8}^{+2.6}$
Scale Height (au)	${16.2}_{-2.0}^{+1.7}$
Dust Mass (M⊙)	${0.00057}_{-0.00022}^{+0.00017}$
Surface Density α	$-{1.77}_{-0.14}^{+0.94}$
Flaring β	${1.19}_{-0.08}^{+0.09}$

Note. Best-fit values for the covariance likelihood estimation mode of the MCMC.

Download table as:  ASCIITypeset image

The best-fit parameters provide a compromise between the image and SED fit. Therefore, this combined fit to the SED and image is not as favorable as if the fits had been performed separately on each individual data set. For example, the best-fit model SED underpredicts the flux in the 20–100 μm region of the SED (by a factor of 20 at 20 μm and 1.5 at 70 μm), while the best-fit model image underpredicts the flux ratio between the top and bottom nebulae by a factor of ∼4. Either of these could have been improved individually if we had only optimized the fit for that metric alone. Models that best fit the image tend to overpredict the disk flux at all wavelengths, while the models that best fit the SED tend to produce images that have very steep surface density profiles, which remove the diffuse material on the outer edges of the disk provided by the tapered edge.

The apparent disagreement is likely a result of some limitations in the disk model. If the opacity of the dust grains in the disk was decreased, the optically thick/thin boundary would move to shorter wavelengths, recovering some of the flux in the several tens of micron range of the SED. However, to improve the flux ratio between the top/bottom nebulae in the modeled image, we would need to move the inclination farther from edge-on and/or change the scattering properties of the grains (i.e., increase the forward-scattering or decrease the dust albedo), which would most likely necessitate an increase in the dust opacity.

The best-fit disk dust mass is ${0.00057}_{-0.00022}^{+0.00017}\,{M}_{\odot }$, which corresponds to a disk mass of 0.057 M_⊙ (assuming the standard ISM gas-to-dust ratio of 100). This is 16% of the stellar mass for a 0.35 M_⊙ star like ESO-Hα 569, a surprisingly high disk-to-star mass ratio.

In the grid fit and the initial MCMC run using the ${\chi }^{2}$ estimator, the fit to the disk mass relied heavily on the 870 μm measurement. Given the surprisingly high mass estimate, we speculated that the 870 μm photometry might be in some way compromised (for instance, contaminated by excess flux from a background source). To test the dependence of the derived mass on this measurement, another MCMC run was tested excluding this data point, but the overall fit still preferred high disk masses. Subsequent to these initial MCMC runs, Dunham et al. (2016) published their ALMA 2.8 mm continuum observations, from which they found a total disk mass (gas+dust assuming a gas-to-dust ratio of 100) of $0.057\pm 0.002\,{M}_{\odot }$. The excellent agreement between these independent results (their estimate and our result from fits without including their 2.8 mm data point) provides increased confidence in the apparently high mass of this disk. Our final MCMC fit using the covariance framework included the 2.8 mm measurement as well as the other photometry.

A key aspect here is the assumed dust opacity. The mass estimate by Dunham et al. (2016) was made under the assumption that the disk is optically thin, in which case the mass may be directly computed via $M(\mathrm{gas}\,+\,\mathrm{dust})=\tfrac{{F}_{\nu }{D}^{2}}{{\kappa }_{\nu }{B}_{\nu }(T)}$. The derived mass thus depends on both the opacity and the disk temperature. Dunham et al. (2016) assumed a disk average temperature of T = 10 K. Our best-fit MCFOST model yields the disk internal temperature as a byproduct of the MCMC radiative transfer calculation, and the results are fairly consistent: a calculated disk midplane temperature of 10 K at 100 au, increasing inwards to 30 K at 5 au. The larger source of potential systematic bias in the disk mass is thus the assumed dust opacity. Dunham et al. (2016) use a dust opacity characteristic of coagulated dust grains with thin icy mantles ($\kappa =0.23$, cm² g⁻¹; Ossenkopf & Henning 1994). Our model uses olivine dust as described in Section 3.2, which yields a similar opacity at 2.8 mm within a factor of 2. But other results can easily be obtained. For instance, if we instead adopt the dust opacity law from Beckwith et al. (1990) (${\kappa }_{\nu }=0.03$ cm² g⁻¹ at 870 μm) and use T = 20 K, then that yields an estimated disk mass $M(\mathrm{gas}+\mathrm{dust})=\tfrac{{F}_{\nu }{D}^{2}}{{\kappa }_{\nu }{B}_{\nu }(T)}=0.0055{M}_{\odot }$, a factor of 10 lower (again assuming the standard gas-to-dust ratio of 100). Better constraints on the dust particle properties and thus the millimeter opacities could help clarify the true mass of this disk. Some disagreement between predictions from the 2.8 mm and 870 μm continuum measurements may be unsurprising since derived dust masses from submillimeter data may be biased downwards in the case of EODs if they begin to become optically thick at that wavelength. That said, the good agreement between the Dunham et al. (2016) millimeter-continuum-derived mass and the result from our fit to the full SED and the HST scattered-light image seems to indicate that the relative importance of absorption/emission and scattering of the dust model (which includes, but is not limited to, the dust albedo) used here is a reasonable approximation.

Lastly, we note that a reduced gas-to-dust ratio would of course directly affect the inferred total gas+dust to star mass ratio, but the available observations do not provide any evidence toward (or against) such a hypothesis.

The best-=fit scale height of ${16.2}_{-2.0}^{+1.7}$ au (at 100 au) is consistent with the low mass of the central star. For a disk that is pressure-supported and vertically isothermal with temperature, the Gaussian vertical density distribution is described by Equation (5) (Burrows et al. 1996):

$$\begin{eqnarray}&&H(r)=\sqrt{\displaystyle \frac{{k}_{B}T(r){r}^{3}}{{{GM}}_{\mathrm{star}}\mu {m}_{p}}},\end{eqnarray} \tag{ 5 }$$

where we assume a reduced mass (μ) of 2.3 If we adopt the best-fit scale height value of 16.18 au at a reference radius of 100 au and calculate the temperature of the disk at this radius, we obtain T ∼ 23 K. This disk temperature agrees well with observations of other edge-on disks (e.g., HH 30; Burrows et al. 1996).

Additionally, MCFOST is capable of producing the temperature structure within the disk along with the images and SEDs. This can be used as a cross-check on the physical self-consistency of our best-fit model parameters. The mass-averaged temperature (across the vertical direction) for our best-fit model at the reference radius (100 au) is T = 29 K. Surface effects that are exacerbated in scattered light could account for the slight discrepancy between the analytically and numerically estimated disk temperatures, as the surface gas is superheated by stellar radiation. The agreement between the dust scale height inferred from the image and the gas scale height computed from the model suggests that the dust grains are well-mixed vertically with very little dust settling, at least for the small dust particles (≲10 μm) that dominate the opacity at visible wavelengths.

The best-fit flaring exponent was $\beta ={1.19}_{-0.08}^{+0.09}$. Kenyon & Hartmann (1987) provide an analytical model for the temperature profile of a flared disk wherein the surface layers are heated by the direct stellar radiation and the energy is re-radiated thermally. Assuming the gas and dust are well-mixed vertically, and that the incident angle of the stellar radiation on the flared surface is small, $T(R)=T({R}_{0}){\left(\tfrac{R}{{R}_{0}}\right)}^{2\beta -3}$. We fit the modeled mass-averaged disk temperature profile to this analytic solution, revealing that a flaring exponent of $\beta =1.29$ is preferred. This value is consistent within ∼1σ of the model preferred value. The best-fit value is slightly shallower than has been predicted for other young, flared disks with $\beta =1.3\mbox{--}1.5$ (e.g., Chiang & Goldreich 1997). This could be an indication of early dust settling in the disk, decoupling the dust and gas, and changing the disk thermal pressure profile. In this model, dust particles of all sizes are assumed to be evenly distributed vertically throughout the disk. An investigation into the effect of settling of larger grains to the disk midplane is left for future work.

The surface density exponent is best fit by $\alpha =-{1.77}_{-0.14}^{+0.94}$, which is near the lower edge of the allowed parameter space. However, allowing for steeper surface density profiles would push the models into a highly unphysical range. The SED favors a very steep surface density profile (also seen for HV Tau C; Duchêne et al. 2010), while the images favor a shallow profile with a more gradual taper at the disk edge. It is possible that the steep best-fit surface density profile is a reaction to the large disk masses required to fit the millimeter data in the SED, whereby mass is being concentrated in the center of the disk, where our data set is poorly equipped to constrain the disk properties. The SED was not expected to have a strong dependence on the surface density slope. The disk is presumed to be very optically thick across most of the IR portion of the SED. Consequently, the surface density profile would not impact the location of the disk scattering surface, which is intercepting and re-radiating light from the central star. It is possible that a degeneracy between the surface density exponent and some other star/disk property is influencing this fit (e.g., stellar luminosity, dust albedo, etc.).

It is unexpected that a disk surface density power law would be steeper than the $\alpha =-1.5$ value for the minimum mass solar nebula (Weidenschilling 1977). Indeed, Andrews & Williams (2007) conducted a resolved submillimeter continuum survey of circumstellar disks and find a mean value of $\alpha =-0.5$. Instead, this steep profile is probably indicative of some shortcoming in our model parameterization. Invoking separate power laws for the inner and outer regions of the disk may provide a solution, but is beyond the scope of this paper. Although we have spatially resolved images at optical wavelengths, the disk is highly optically thick, causing any effects of radial density gradients to be undetectable. Characterizing these would require resolved images at wavelengths where the disk is optically thin (e.g., resolved millimeter-continuum images, though it is uncertain if the disk is truly optically thin at these wavelengths). Scattered-light imaging alone simply does not constrain the surface density exponent in the innermost regions of the disk. Previous studies of the radial structure of protoplanetary disks observed in millimeter continuum find surface density profiles that are generally shallower than presented here, although there are a few exceptions (e.g., DG Tau, GM Aur; Guilloteau et al. 2011). Estimates of the surface density distributions inferred from resolved millimeter data at different wavelengths vary widely (Isella et al. 2010), suggesting that these disks are not optically thin even in the 1–3 mm range.

The best-fit dust mass ($0.00057\,{M}_{\odot }$ or $190\,{M}_{\oplus }$) and the associated total (gas + dust) disk mass ($60\,{M}_{\mathrm{Jup}}$, assuming a gas-to-dust ratio of 100) imply a disk mass to star mass ratio (${M}_{D}/{M}_{\mathrm{star}}$) significantly higher than expected for its age and spectral type. Williams & Cieza (2011) provide a review of protoplanetary disks and report a relatively flat distribution of disk masses when spaced logarithmically, with a sharp drop outside of ∼50 M_Jup and an average disk mass to host a stellar mass ratio of 0.01 albeit with a large scatter. The median mass (assuming a gas-to-dust ratio of 100) of disks around GKM spectral type hosts is 5 M_Jup (implying a dust mass of ∼16 ${M}_{\oplus }$).

This trend of low ${M}_{D}/{M}_{\mathrm{star}}$ mass ratios seems to continue for low-mass stars. Van der Plas et al. (2016) conducted a survey of disk masses for low-mass stars with ALMA, finding a range of masses between 0.1 and 1 ${M}_{\oplus }$ for their eight targets. One target in their sample, Allers 8 (an M3 star with a mass of 0.34 M_⊙), has similar stellar parameters to ESO-Hα 569, but a significantly lower dust mass of 1.05 ${M}_{\oplus }$. However, the bulk of the disks in their sample are located in the Upper Scorpius SFR (∼10 Myr; David et al. 2016) and are older than our target. The authors have an additional data set for the younger Taurus SFR, with preliminary estimates for the dust mass upper limit of 25 ${M}_{\oplus }$ for a sample of stars with an earliest spectral type of M4 (K. Ward-Duong 2017, private communication). Additionally, Andrews et al. (2013) conduct a survey of the protoplanetary disks with low-mass hosts (spectral types earlier than M8.5) in the Taurus SFR and find slightly higher disk masses. The authors estimate the disk masses from their millimeter-wave continuum luminosity, and find that the median disk mass to stellar mass ratio is 0.3%, with very few disks having a ratio of ≥10%. Targets in their sample in the M3–M4 spectral type range have disk dust masses of 2–17 ${M}_{\oplus }$, with an average of 9 ${M}_{\oplus }$.

While uncommon, protoplanetary disks with large disk masses are not unprecedented. Duchêne et al. (2010) model scattered-light images and SEDs for the HV Tau C system and find a best-fit dust mass of ${M}_{\mathrm{dust}}\geqslant {10}^{-3}\,{M}_{\odot }$, which gives ${M}_{D}/{M}_{\mathrm{star}}\sim 0.2$ (meaning the disk is 20% of the mass of the central star), assuming a gas-to-dust ratio of 100. Likewise, Duchêne et al. (2003) model a millimeter image of the HK Tau B protoplanetary disk and get a best-fit total disk mass of ${M}_{\mathrm{disk}}\simeq 2\times {10}^{-2}\,{M}_{\odot }$, which gives ${M}_{D}/{M}_{\mathrm{star}}\sim 0.04$ (4% of the stellar mass). Glauser et al. (2008) present an in-depth study of the IRAS 04158 + 2805 disk using images in the optical and NIR, polarization maps in the optical and mid-IR, and X-ray spectra. The dust mass is constrained to be ${M}_{\mathrm{dust}}=1.0\mbox{--}1.75\times {10}^{-4}\,{M}_{\odot }$, which also gives ${M}_{D}/{M}_{\mathrm{star}}\,=\sim 0.04$ (4%). All three disks are in the Taurus SFR and the first two disks above are in multiple systems. Likewise, all of these sources are viewed edge-on. It is possible that the large inferred disk masses could be the result of a selection effect (observations of edge-on disks are only sensitive to the most massive disks), or some artifact of our fitting method which compensates for missing physics by placing more mass in the disk. The fit for the dust mass is driven by the SED, but the spectral coverage is poor in the millimeter. The mass estimates could be reduced by including larger opacities in the millimeter, for instance by adding amorphous carbon into the mixture, or by using a more complex, nonuniform particle distribution.

Throughout this paper, we assumed a gas-to-dust mass ratio of 100, as is typical of other young disks and the ISM. However, very recent work by Long et al. (2017) estimate the gas mass around ESO-Hα 569 from ALMA ¹³CO line emission and find only $\sim 1.3\,{M}_{\mathrm{Jup}}$ of gas mass in the disk, although optical depth effects and details of the CO freeze-out are likely to introduce major sources of uncertainty. Combined with our own dust mass estimate, this gives an uncharacteristically low gas-to-dust ratio of only ∼2. Although gas depletion in the disk would lower the unusually high best-fit total disk mass, the flared appearance strongly confirms this is a young pressure-supported gas+dust disk. A gas-to-dust ratio of 2:1 is suggestive of a later evolutionary stage. This disagreement highlights the challenges of measuring disk masses for EODs, which are generally optically thick even at millimeter wavelengths.

The MCMC radiative transfer fit prefers a disk with an abnormally large disk mass that is ∼16% of the mass of the central star. We investigate the stability of the disk via the Toomre Q parameter,

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{s}\kappa }{\pi G{\rm{\Sigma }}},\end{eqnarray} \tag{ 6 }$$

where c_s is the sound speed in the disk, κ is the epicyclic frequency, and Σ is the surface density profile of the disk. For a vertically isothermal disk with a Keplerian velocity, $\kappa ={\rm{\Omega }}=\sqrt{\tfrac{{{GM}}_{\mathrm{Star}}}{{R}^{3}}}$ and a sound speed ${c}_{s}=\sqrt{\tfrac{{k}_{B}T}{\mu {m}_{p}}}$, where a reduced mass (μ) of 2.3 was assumed. Figure 12 shows the radial profile of the Toomre Q parameter. It shows that the disk appears to be stable at all radii.

Zoom In Zoom Out Reset image size

It is worth noting here that a good constraint is not expected on the properties of the inner regions of the disk from scattered-light imaging and the SED alone. Any change in the interior structure (e.g., an inner wall, spiral structure, or a broken surface density power law) of the disk would affect stability. Each of these mechanisms would increase the variability of the system, possibly accounting for the observed variability in several of the photometric points included in the SED.

Rodgers-Lee et al. (2014) conducted a survey of disks in Cha I as identified from IR excesses in the SEDs. For 34 objects, disk masses were estimated. The median of the distribution of disk masses is 0.005 M_⊙, which corresponds to 0.5% of the stellar mass, while the tail of the distribution stretches to 0.1 M_⊙ for more massive central stars. ESO-Hα 569 is a clear outlier with 10 times more mass than the median value.

The Luhman (2007) survey of Cha I names six members as likely edge-on disk candidates because they are underluminous for their spectral type and are seen in scattered light (CHSM 15991, T14A, ISO 225, ESO-Hα 569 and 574, and Cha J11081938–7731522). One of those objects, Cha J11081938–7731522, appears extended in their survey with a butterfly morphology, providing further support that these targets are all likely edge-on disks. Two members of that list, ESO-Hα 569 and 574, were observed in our HST campaign, which confirmed that both are edge-on protoplanetary disks.

Luhman et al. (2008) use Spitzer colors to estimate the disk fraction as a function of stellar mass. For stars of spectral type between K6 and M3.5, the disk fraction in Cha I is 0.64 ± 0.06 disks per star. If we multiply this fraction by the fraction of disks expected to have inclinations between 75° and 90°, roughly 17% of stars with young disks should host edge-on disks. However, a recent survey of 44 YSOs hosting circumstellar disks detectable with Herschel found only two edge-on disks (Rodgers-Lee et al. 2014) as classified from the SEDs. Although the sample size of this survey is small, this surprising lack of known EODs is a common phenomena seen for many SFRs (Stapelfeldt et al. 2014) and was one of the key motivating factors for our HST survey. While that program doubled the number of known EODs, the increased sample remains smaller than would be predicted from purely geometrical grounds. This suggests that many disks must be near edge-on but with insufficient material and/or vertical extent to block the direct light of the star. Flatter disks, with lower H/R values than ESO-Hα 569, would only appear edge-on for a narrower inclination range. For instance, if the "typical" young disk is flared enough to only occult its star within 5°, considering a range from 85° to 90° would give an edge-on disk fraction per star of 4%, more in line with what is observed.

Alternatively, this could suggest that the "typical" double-peaked SED assumed for edge-on disks may only be present for disks with an unusually high disk mass. The targets for this edge-on disk survey were selected based on the shape of the SEDs—specifically, targets with a doubled-peaked SED, where the stellar peak flux was of the same order as the dust peak flux in the IR. Figure 13 shows the effect of changing dust mass on the structure of the SED and the scattered-light image for a fixed inclination. Disk masses shown are for the best-fit disk mass divided by factors of 3, 10, 30, and 100. After dividing by a factor of 10 (for a more reasonable ${M}_{D}/{M}_{\odot }\sim 0.016$), the double-peaked structure disappeared, and we did not include this target in our sample. It is possible that this could account for the relative lack of known edge-on disks; the selection metrics used are biased toward detecting only the most massive disks, as they require fairly large line-of-sight opacities.

Zoom In Zoom Out Reset image size

It is therefore possible that the edge-on disk detections thus far are outliers in the population of young disks. Double-peaked SEDs alone are an insufficient indicator of the edge-on disk fraction, and images in scattered light or thermal emission with high spatial resolution are required to determine the true nature of these objects. Existing surveys of young, nearby SFRs tend to have selection biases toward more face-on systems and are dependent on the cloud properties and the science drivers of the survey. In order to determine the true edge-on disk fraction and to confirm or deny that the high disk mass of ESO-Hα 569 is indeed representative of the population of protoplanetary disks, a uniform sample of disk observations at sufficiently high angular resolution will be required.