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

Schuyler G. Wolff, Marshall D. Perrin, Karl Stapelfeldt, Gaspard Duchêne, François Ménard, Deborah Padgett, Christophe Pinte, Laurent Pueyo, William J. Fischer

► To cite this version:

Schuyler G. Wolff, Marshall D. Perrin, Karl Stapelfeldt, Gaspard Duchêne, François Ménard, et al. Hubble Space Telescope Scattered-light Imaging and Modeling of the Edge-on Protoplanetary Disk ESO-Hα 569. The Astrophysical Journal, 2017, 851, 10.3847/1538-4357/aa9981. insu-03692463

HAL Id: insu-03692463

https://insu.hal.science/insu-03692463v1

Submitted on 23 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Schuyler G. Wolff<fnref n="1" />,<fnref n="2" />,*, Marshall D. Perrin<fnref n="3" />,*, Karl Stapelfeldt<fnref n="4" />,<fnref n="5" />, Gaspard Duchêne<fnref n="6" />,<fnref n="7" />,*, Francois Ménard<fnref n="7" />,*, Deborah Padgett<fnref n="5" />,*, Christophe Pinte<fnref n="8" />,<fnref n="9" />,*, Laurent Pueyo<fnref n="3" />, and William J. Fischer<fnref n="3" />,†

Received 2017 June 27; revised 2017 November 7; accepted 2017 November 7; published 2017 December 12

## Abstract

We present new *Hubble Space Telescope (HST)* Advanced Camera for Surveys observations and detailed models for a recently discovered edge-on protoplanetary disk around ESO-Hα 569 (a low-mass T Tauri star in the Cha I star-forming region). Using radiative transfer models, we probe the distribution of the grains and overall shape of the disk (inclination, scale height, dust mass, flaring exponent, and surface/volume density exponent) by model fitting to multiwavelength (F606W and F814W) *HST* observations together with a literature-compiled spectral energy distribution. A new tool set was developed for finding optimal fits of MCFOST radiative transfer models using the MCMC code *emcee* to efficiently explore the high-dimensional parameter space. It is able 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 those derived properties. We confirm that ESO-Hα 569 is an optically thick nearly edge-on protoplanetary disk. The shape of the disk is well-described by a flared disk model with an exponentially tapered outer edge, consistent with models previously advocated on theoretical grounds and supported by millimeter interferometry. The scattered-light images and spectral energy distribution are best fit by an unusually high total disk mass (gas+dust assuming a ratio of 100:1) with a disk-to-star mass ratio of 0.16.

**Key words:** protoplanetary disks – radiative transfer

## 1. Introduction

We seek to understand the initial conditions for planet formation and the physical processes that contribute to the assembly of planets by measuring the properties of young protoplanetary disks. The unique geometry of edge-on circumstellar disks provides a valuable opportunity to study the detailed disk structure, as the bright central star is occulted from view and thus does not pose a contrast problem. The width of the disk’s dark lane (the vertical extent of the $τ = 1$ surface), outer radius, and degree of flaring can be directly measured, and the scale height of the disk can be related to the local disk temperature (Watson et al. 2007). Stapelfeldt (2004) provides a review of the observational advantages of targeting edge-on disks. Previous studies of edge-on disks (EODs) have measured disk inclinations and dust masses from a combination of scattered-light images and millimeter-continuum maps (Wolf et al. 2003; Sauter et al. 2009). Additionally, the change in the dust lane thickness with wavelength allows dust grain properties to be derived (Cotera et al. 2001; Wood et al. 2008; Duchêne et al. 2010; McCabe et al. 2011). However, the sample of edge-on disks with high-resolution observations remains relatively small.

ESO-Hα 569, a young M2.5 star embedded in the Chameleon I star-forming region (SFR), was imaged as part of an *Hubble Space Telescope (HST)* observation program designed to double the sample of edge-on protoplanetary disks

for which high-resolution scattered-light images have been obtained. The sample for the survey was chosen from *WISE* and *Spitzer* surveys of nearby SFRs, which allow the identification of new candidate edge-on disks from their characteristic double-peaked spectral shape. *HST* program 12514 in Cycle 19 obtained high-resolution optical imaging of the top 21 candidates, including the data presented in the current work (Stapelfeldt et al. 2014).

Several of the targets in this sample of edge-on protoplanetary disks, including ESO-Hα 569, are known members of the Chameleon I (Cha I) SFR. Distances to Cha I, one of the nearest SFRs, have been determined in a variety of ways, including through zero-age main-sequence fitting and *Hipparcos* parallaxes of members. Whittet et al. (1997) provide a review of the results and combine measurements to arrive at a distance of $160 \pm 15$ pc. Bertout et al. (1999) confirm this distance after cross-correlating the Herbig & Bell and *Hipparcos* catalogs. Belloche et al. (2011, see their Appendix B1) present a more detailed review of the Cha I distance measurements. Age estimates for Chamaeleon I range from 1 to 2 Myr (Baraffe et al. 1998; Chabrier et al. 2000). The Cha I SFR is characterized by a relatively high extinction with an observed maximum of $A_V \sim 10$ (Cambresy et al. 1997). Such a high extinction would suppress the blue side of the spectral energy distribution (SED) of a young stellar system. The initial mass function for Cha I has a maximum mass of 0.1–0.15 $M_\odot$ (Luhman 2007), while the total mass of Cha I is $\sim 1000 M_\odot$ (Boulanger et al. 1998).

<footnote n="1">¹ Johns Hopkins University Baltimore, MD 21218, USA; wolff@strw.leidenuniv.nl</footnote>

<footnote n="2">² Leiden Observatory, Leiden University 2300 RA Leiden, The Netherlands</footnote>

<footnote n="3">³ Space Telescope Science Institute Baltimore, MD 21218, USA</footnote>

<footnote n="4">⁴ Exoplanet and Stellar Astrophysics Laboratory, NASA Goddard Space Flight Center, USA</footnote>

<footnote n="5">⁵ Jet Propulsion Laboratory, California Institute of Technology, USA</footnote>

<footnote n="6">⁶ University of California, Berkeley, USA</footnote>

<footnote n="7">⁷ Université Grenoble-Alpes, CNRS Institut de Planétologie et d'Astrophysique (IPAG) F-38000 Grenoble, France</footnote>

<footnote n="8">⁸ UMI-FCA, CNRS/INSU, (UMI 3386), France</footnote>

<footnote n="9">⁹ Dept. de Astronomía Universidad de Chile Santiago, Chile</footnote>
The remainder of this introduction summarizes prior observations of ESO-Hα 569. In Section 2, we present high-resolution HST scattered-light observations of the ESO-Hα 569 protoplanetary disk and an SED compiled from the literature. In Section 3, radiative transfer modeling efforts to fit these observations to a variety of disk properties are discussed. Both a grid and Monte Carlo Markov Chain (MCMC) approach were used to explore parameter space, and results are given in Sections 3.4 and 3.6. Section 4 discusses these results, including the gravitational stability of the system and places ESO-Hα 569 in context with previous disk observations. Lastly, Section 5 provides a summary and the conclusions.

### 1.1. Prior Studies of 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 ≈ 5000 spectrum from 0.6 to 0.9 µ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 ± 0.08 L_{⊙}$. 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 \ge 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 ~87°1, $L_{bol} = 0.8 ± 0.4 L_{⊙}$, $M_{star} = 0.33 ± 0.03 M_{⊙}$, and $R_{star} = 2.5 ± 0.6 R_{⊙}$, and give an upper limit for the submillimeter disk mass of $0.005 M_{⊙}$. Rodgers-Lee et al. (2014) included Herschel data and found a best-fit inclination of 81°4. More recently, Pascucci et al. (2016) provided 1.3 mm continuum data, which correspond to a disk mass estimate of $0.0046 M_{⊙}$, using an assumed opacity of $\kappa = 2.3 \text{ cm}^2 \text{ g}^{-1}$, a gas-to-dust ratio of 100, and a disk temperature of 20 K.

## 2. Observations

### 2.1. HST Scattered-light Images

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 ~20× 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.

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'' \times 2''$ and was chosen to encompass all disk flux with surface brightness $\ge 3\sigma$ above the background noise. The measured flux density is $0.058 \pm 0.001$ mJy for F606W and $0.21 \pm 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 \pm 0.0''05$ which corresponds to $125 \pm 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.

### 2.2. HST Jet Outflow Images

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''$ (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.0''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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** HST images of the protoplanetary disk ESO-Hα 569. Top: F814W. Bottom: F606W. Both images show the dark dust lane and asymmetries between the top and bottom of the disk, while only F606W establishes the presence of an outflow jet. The 100 au scale bar corresponds to an angular scale of 0.′′625.</figcaption></figure>

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.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** A wider F606W filter image displaying the diffuse nebula extending outward from the disk. The direction of the Hα filament HH 919 is shown by the arrow. The jet lines up well with the reported position of HH 919, consistent with ESO-Hα 569 being the origin of this outflow. The 500 au scale bar corresponds to an angular scale of 3.′′125.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** An image of the jet created by subtracting the F814W image from the F606W image. Contours are drawn from 0.01 to 0.19 mJy arcsec⁻² in intervals of 0.03 mJy arcsec⁻². The 100 au scale bar corresponds to an angular scale of 0.′′625.</figcaption></figure>

## 2.3. Spectral Energy Distribution

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).

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
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Spectral energy distribution for ESO-Hα 569 with upper limits indicated by triangles. The SED exhibits the double-peaked structure typical of an optically thick, edge-on disk. The values were compiled from the literature with more information given in Table 1. The stellar spectrum for an M2.5 star with $T_{\text{eff}} = 3500$ K is overplotted.</figcaption></figure>

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 χ²-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) × 10⁻¹⁶ 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.

### 3. Model Fitting

The scattered-light images and full SEDs together provide a comprehensive data set for ESO-Hα 569 against which properties of the central star and surrounding disk can be tested. The disk geometry can be directly measured from the images, and the distribution of the dust grains within the disk is

traced by the SED and disk morphology. To characterize this system, disk models were constructed to explore parameter space with direct comparisons to the observations. The next section presents the context and challenges for radiative transfer modeling of complex disk structures.

#### 3.1. Radiative Transfer Modeling and Model Fitting of Circumstellar Disks

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.
<table><caption>**Table 1**
Spectral Energy Distribution Photometry and References</caption><thead><tr><th>λ(μm)</th><th>Flux (mJy)</th><th>Source</th><th>Instrument</th><th>Bandwidth (μm)</th><th>Angular Resolution</th><th>Date</th></tr></thead><tbody><tr><td>0.551</td><td>0.030 ± 0.004</td><td>Robberto et al. (2012)</td><td>HST WFPC2</td><td>0.14</td><td>0"0996</td><td>2009 Apr 27</td></tr><tr><td>0.606</td><td>0.058 ± 0.001</td><td>This work</td><td>HST ACS</td><td>0.27</td><td>0"05</td><td>2012 Mar 09</td></tr><tr><td>0.814</td><td>0.21 ± 0.01</td><td>This work</td><td>HST ACS</td><td>0.31</td><td>0"05</td><td>2012 Mar 09</td></tr><tr><td>1.235</td><td>0.66 ± 0.05</td><td>Skrutskie et al. (2006)</td><td>2MASS</td><td>0.16</td><td>~5"</td><td>2000 Jan 25</td></tr><tr><td>1.662</td><td>0.97 ± 0.08</td><td>Skrutskie et al. (2006)</td><td>2MASS</td><td>0.25</td><td>~5"</td><td>2000 Jan 25</td></tr><tr><td>2.15</td><td>0.98 ± 0.09</td><td>Skrutskie et al. (2006)</td><td>2MASS</td><td>0.26</td><td>~5"</td><td>2000 Jan 25</td></tr><tr><td>3.6</td><td>0.58 ± 0.03</td><td>Luhman et al. (2008)</td><td>Spitzer IRAC</td><td>0.75</td><td>~2"</td><td>2004 Jul 04</td></tr><tr><td>4.5</td><td>0.57 ± 0.05</td><td>Luhman et al. (2008)</td><td>Spitzer IRAC</td><td>1.02</td><td>~2"</td><td>2004 Jul 04</td></tr><tr><td>5.8</td><td>0.58 ± 0.05</td><td>Luhman et al. (2008)</td><td>Spitzer IRAC</td><td>1.43</td><td>~2"</td><td>2004 Jul 04</td></tr><tr><td>8.0</td><td>0.67 ± 0.05</td><td>Luhman et al. (2008)</td><td>Spitzer IRAC</td><td>2.91</td><td>~2"</td><td>2004 Jul 04</td></tr><tr><td>3.4</td><td>0.63 ± 0.02</td><td>Cutri et al. (2012)</td><td>WISE</td><td>0.66</td><td>6"1</td><td>2010 Feb 13, 20</td></tr><tr><td>4.6</td><td>0.71 ± 0.02</td><td>Cutri et al. (2012)</td><td>WISE</td><td>1.04</td><td>6"4</td><td>2010 Feb 13, 20</td></tr><tr><td>12</td><td>0.65 ± 0.07</td><td>Cutri et al. (2012)</td><td>WISE</td><td>5.51</td><td>6"5</td><td>2010 Feb 13, 20</td></tr><tr><td>22</td><td>7.5 ± 0.89</td><td>Cutri et al. (2012)</td><td>WISE</td><td>4.10</td><td>12"0</td><td>2010 Feb 13, 20</td></tr><tr><td>24</td><td>8.36 ± 0.77</td><td>Luhman et al. (2008)</td><td>Spitzer MIPS</td><td>5.3</td><td>6"</td><td>2004 Apr 08</td></tr><tr><td>70</td><td>107 ± 10.8</td><td>Luhman et al. (2008)</td><td>Spitzer MIPS</td><td>19</td><td>18"</td><td>2004 Apr 08</td></tr><tr><td>70</td><td>200 ± 100</td><td>Winston et al. (2012)</td><td>Herschel PACS</td><td>25</td><td>5"8</td><td>2011 Jun 23</td></tr><tr><td>160*</td><td>200 ± 200</td><td>Winston et al. (2012)</td><td>Herschel PACS</td><td>85</td><td>12"0</td><td>2011 Jun 23</td></tr><tr><td>250*</td><td>100 ± 100</td><td>Winston et al. (2012)</td><td>Herschel SPIRE</td><td>25</td><td>18"</td><td>2011 Jun 23</td></tr><tr><td>350*</td><td>50 ± 50</td><td>Winston et al. (2012)</td><td>Herschel SPIRE</td><td>25</td><td>25"</td><td>2011 Jun 23</td></tr><tr><td>500*</td><td>50 ± 50</td><td>Winston et al. (2012)</td><td>Herschel SPIRE</td><td>25</td><td>37"</td><td>2011 Jun 23</td></tr><tr><td>870</td><td>72 ± 14</td><td>Belloche et al. (2011)</td><td>APEX/LABOCA</td><td>150</td><td>19"2</td><td>2008 May</td></tr><tr><td>2830</td><td>3.2 ± 0.1</td><td>Dunham et al. (2016)</td><td>ALMA</td><td>55</td><td>~2"</td><td>2013 Nov 29–2014 Mar 08</td></tr></tbody></table>

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

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.
## 3.2. Radiative Transfer Modeling with MCFOST

For this work, the MCFOST radiative transfer code (Pinte et al. 2006, 2009) was used to construct SEDs and 0.8 µm scattered-light images for each of the models. The 0.6 µ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, $\Sigma$, described by a power-law distribution in radius given by $\Sigma = \Sigma_0(R/R_0)^\alpha$, where $\alpha$ is termed the surface density exponent and $R_0$ 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 $\Sigma$ falls off exponentially with some critical radius $R_c$ of material outside of the disk:<fnref n="10" />

$$ \Sigma = \Sigma_c \left(\frac{R}{R_c}\right)^\alpha \exp\left[-\left(-\frac{R}{R_c}\right)^{2+\alpha}\right]. \quad (1) $$

For this work, $R_c = R_{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 $\beta$ is the flaring exponent describing the curvature of the disk and again $R_0 = 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_{sub} = R_{star}(T_{star}/T_{sub})^{2.1} \sim 0.1$ au, where $T_{sub} = 1600$ 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 ($\beta$), surface density exponent ($\alpha$), 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$ µ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 Brugge-man effective mixing rule.

<table><caption>**Table 2**
Modeled Disk Parameters</caption><thead><tr><th>Parameter</th><th>Grid Values</th><th>MCMC Values</th></tr></thead><tbody><tr><td>Distance (pc)</td><td>160 (Fixed)</td><td>160 (Fixed)</td></tr><tr><td>Outer Radius (au)</td><td>125 (Fixed)</td><td>125 (Fixed)</td></tr><tr><td>Minimum Particle Size (µm)</td><td>0.03 (Fixed)</td><td>0.03 (Fixed)</td></tr><tr><td>Inclination (°)</td><td>60–90</td><td>65–90</td></tr><tr><td>Scale Height (H in au at R = 100 au)</td><td>10, 15, 20, 25</td><td>5–25</td></tr><tr><td>Dust Mass (M in M<sub>☉</sub>)</td><td>10<sup>-4</sup>, 3 × 10<sup>-4</sup>, 10<sup>-3</sup></td><td>10<sup>-5</sup>–10<sup>-3</sup></td></tr><tr><td>Surface Density (α)</td><td>-2.0, -1.5, -1.0, -0.5, 0.0</td><td>-2.0 to 0</td></tr><tr><td>Flaring Exponent (β)</td><td>1.1, 1.2, 1.3, 1.4, 1.5</td><td>1.0–1.5</td></tr><tr><td>Maximum Grain Size<sup>a</sup> (µm)</td><td>100, 1000, 3000</td><td>100–3000</td></tr><tr><td>Weight<sup>b</sup></td><td>...</td><td>0.3–0.7</td></tr><tr><td>Grain Porosity</td><td>0.0, 0.25, 0.5</td><td>...</td></tr><tr><td>Structure</td><td>Disk, Tapered-edge Disk</td><td>...</td></tr></tbody></table>

**Notes.**

a Grain size was kept at a constant value of 100 µm for the covariance-based MCMC run.

b During 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.

For comparison with the observed 0.8 µm scattered-light images, each model image was convolved with a Tiny Tim simulated PSF (Krist 1995). The 0.8 µm observations were masked to select only the pixels with flux values ≥3σ 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.

## 3.3. Initial Exploration of Parameter Space via Grid Search

Our initial modeling used a uniform grid sampling, with the explored parameter space shown in Table 2. For

<footnote n="10">¹⁰ Note that some authors give this equation using the notation $\gamma = -\alpha$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Likelihood distributions from the grid search for the disk inclination and scale height computed from the model $\chi^2$ values for the 0.8 µm image (red), the SED (blue) and for the combined data set (gray). The image and SED results favor different regions of parameter space. The sampling of the grid approach is sparse and does not provide an adequate estimate of the uncertainties.</figcaption></figure>

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 µm image and SED for each model along with the combined total $\chi^2_{tot} = \chi^2_{0.8 \mu m} + \chi^2_{SED}$. 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).

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.

### 3.4. Results and Conclusions from Grid Search

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 75°, and a scale height of 20 au at a reference radius of 100 au. The preferred maximum particle size is 3000 µm, the dust mass is $3 \times 10^{-4} M_\odot$, the flaring exponent $\beta$ is 1.3, and the surface density exponent $\alpha$ is -0.5. The separate SED and image fits for the $\alpha$ and $\beta$ exponents favor opposing extremes of parameter space, but the combined $\chi^2_{tot}$ 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.

#### 3.4.1. Porosity

Porous grains were initially included in the modeling parameters to provide a better fit to the flux ratio between the top and bottom disk nebulae. Porous grains are generally more forward-scattering, which would increase the flux ratio without needing to increase the line-of-sight inclination. However, the SED fitting strongly favored non-porous grains. A porosity of $\ge 0.5$ produced a strong dip in the SED around the 10–20 µm silicate feature that was not observed for this target. The overall SED+image fits also favor the non-porous grains, though not as strongly as the SEDs alone. The flux ratio issue was subsequently solved by invoking a tapered-edge surface density model for the disk structure. For subsequent modeling, only non-porous grains were used.

#### 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.

### 3.5. Model Optimization via MCMC

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

$$P(\Theta|D) \propto P(D|\Theta)P(\Theta), \quad (2)$$

where $D$ represents the observations and $\Theta$ the free parameters in the model. Here, $P(D|\Theta)$ is the likelihood of the data given the model and $P(\Theta)$ 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
<figure><img src="image_8.png" /><figcaption>Figure 6. Surface brightness profiles for two vertical image cuts through the data (left) and through the sharp-edged (top middle panel) and tapered-edged (bottom middle panel) disk models. The residuals for the two models are plotted on the same scale (smaller panels at right). The tapered-edge model does a much better job of fitting the shape of the disk, especially the depth of the disk midplane and the diffuse outer regions.</figcaption></figure>

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,<fnref n="11" /> 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:

$$ \ln[P(D|\Theta)] = -0.5N \ln 2\pi + \sum_{i=1}^{N} (-\ln|\sigma_i|) - \frac{1}{2} \chi_i^2. \quad (3) $$

Here, $N$ is the number of data points and $\sigma$ is our uncertainty. The MCMC approach inherently requires a single goodness-of-fit metric, so it is essential to combine the SED and image

<footnote n="11">¹¹ https://github.com/swolff9/mcfost-python</footnote>
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^2_{SED} = 1.3$, $\chi^2_{0.8 \ \mu m} = 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):

$$ \ln[P(D|\Theta)] = -\frac{1}{2}(R^T C^{-1} R + \ln[\det(C)] + N \ln[2\pi]), \quad (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_{pix} \times N_{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} \delta_{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$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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Mean of the autocorrelation of the residuals from subtractions between our 0.8 µm observed image and a randomly selected subset of 1000 model images spanning the range of the priors. Residuals are most strongly correlated between pixels that are horizontally adjacent, as expected for an edge-on disk with its major axis oriented horizontally. The slight anticorrelation in the vertical direction is likely due to dark lane structure between the two bright lobes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Slices through the mean autocorrelation shown in Figure 7. Both the vertical (blue) and horizontal (red) slices along with several Matérn kernels are presented for comparison. We conservatively adopt the wider correlation scale from the horizontal axis to generate the global covariance matrix. It is not unsurprising that the autocorrelation image is more broadly extended in the horizontal direction, where the disk is elongated, than in the vertical, where the gradients in the disk are much sharper.</figcaption></figure>

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.

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

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Top: covariance matrix ($C_{i,j}$) used to compute the log likelihood of the model images given the observations. The matrix combines information about the noise in the observations, the covariances between adjacent pixels, and the pixel mask. Bottom: a zoomed-in region illustrating the contribution of the autocorrelation function between adjacent pixels. To generate this, the 2D 50 × 50 pixel image is first unwrapped into a 1D 2500 pixel array by stacking each row horizontally. The diagonal of the covariance matrix gives the uncertainties associated with each pixel (where $i = j$). The other elements of the covariance matrix dictate the covariances between the various pixel pairs ($i, j$), given by the autocorrelation shown in Figure 8 and depending on the distance between pixels $i$ and $j$ in the 2D detector frame (not in the 1D unwrapped image).</figcaption></figure>

$\chi^2 \sim 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 allowed parameter ranges were adjusted slightly for the MCMC modeling compared to the grid fit. The computation time for the grid modeling depended both on the number of free parameters and on the size of the allowed parameter ranges, while the MCMC modeling time depended only on the number of free parameters. Therefore, we were able to widen the prior distributions for the MCMC modeling, being careful to widen allowable ranges for those parameters that were best fit at the edges of the grid distribution, such as the scale height and disk mass. Parameter ranges are shown in Table 2, column 3. During the IDL grid search modeling phase, we found the image and SED fits both prefer a large maximum grain size. In order to limit the computation time in the MCMC fits, the maximum grain size was fixed to be 3000 µm, and as noted above, the porosity was set at zero.

One downside to the MCMC over the grid search approach is that the chain does not work well with discrete parameter distributions. For example, the abrupt distinction between the tapered- and sharp-edged disk models could not have been tested using MCMC. Given the strong support for the tapered-edge disk model as described in Section 3.4.2, we selected an exponentially tapered outer edge for the MCMC run.

An MCMC run was conducted using the covariance-based log-likelihood goodness-of-fit metric with two temperatures with 50 walkers. Uniform prior distributions were used for all of the parameters (with the dust mass uniformly distributed in log space). The chain was run for $N_{\text{steps}} = 10,000$, with an initial burn-in stage of $N_{\text{burn}} = 0.2 N_{\text{steps}}$. This resulted in a total of 21,000 models requiring ~2 weeks of computation time parallelized over only 10 cores. This was a significant improvement over the grid search approach, which necessitated generating ~200,000 models. As a test of convergence, we compute the integrated autocorrelation times ($\tau_x$) for each of the parameters and use these to estimate the effective sample size, $ESS = N_{\text{samples}}/(2\tau_x)$ (a measure of the effective number of independent samples in the correlated chain). The ESS varied from 761 to 12,075, with the surface density exponent being the least well-constrained parameter. The Monte Carlo standard error for each parameter decreases with increasing effective sample size as $\sigma_i/\sqrt{ESS}$, where $\sigma_i$ is the standard deviation for the posterior distribution (see the discussion in Sharma 2017). For example, to measure the 0.025 quantile to within ±0.01 with a probability of 0.95 requires 936 uncorrelated samples (this corresponds to roughly 10% errors in the best-fit parameter values assuming the tail of the posterior is well-described by a normal distribution), which is achieved for all parameters except the surface density distribution where the 0.025 quantile was only confined to within roughly ±0.0125 (Raftery & Lewis 1992).

### 3.6. Results and Conclusions from 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^{+2.6}_{-4.8}$ degrees, a scale height of $16.2^{+1.7}_{-2.0}$ au at a reference radius of 100 au, a disk dust mass of $0.00057^{+0.00017}_{-0.00022} M_{\odot}$ (assuming a gas-to-dust ratio of 100), a surface density exponent ($\alpha$) of $-1.77^{+0.94}_{-0.14}$, and a flaring exponent ($\beta$) of

<table><caption>**Table 3**
MCMC Best-fit Parameters</caption><thead><tr><th>Parameters</th><th>Best-fit Values</th></tr></thead><tbody><tr><td>Inclination (°)</td><td>83.0<sup>+2.6</sup><sub>-4.8</sub></td></tr><tr><td>Scale Height (au)</td><td>16.2<sup>+1.7</sup><sub>-2.0</sub></td></tr><tr><td>Dust Mass ($M_{\odot}$)</td><td>0.00057<sup>+0.00017</sup><sub>-0.00022</sub></td></tr><tr><td>Surface Density $\alpha$</td><td>-1.77<sup>+0.94</sup><sub>-0.14</sub></td></tr><tr><td>Flaring $\beta$</td><td>1.19<sup>+0.09</sup><sub>-0.08</sub></td></tr></tbody></table>

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

$1.19^{+0.09}_{-0.08}$. 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.

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.

### 3.7. Dust Mass

The best-fit disk dust mass is $0.00057^{+0.00017}_{-0.00022} M_{\odot}$, which corresponds to a disk mass of $0.057 M_{\odot}$ (assuming the standard ISM gas-to-dust ratio of 100). This is 16% of the stellar mass for a $0.35 M_{\odot}$ 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
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Results from the covariance-based MCMC. Top: the model image (middle) corresponding to the best-fit parameters given in Table 3 compared to the 0.8 µm observed image (left). The right panel shows a contour highlighting the shape of the best-fit model disk in red, with contours scaled to the observed 0.8 µm image shown in blue. One hundred randomly chosen models drawn from the MCMC chain are depicted in gray. Bottom: the SED for the same model as above is shown in red and compared to the literature values in blue. The gray curves present the same 100 randomly selected models drawn from the chain. While the MCMC results provide a reasonably good fit to both the image and SED, the compromise between the two data sets, inherent in the covariance framework, lead to imperfect solutions. For example, the best-fit model underpredicts the flux in the 20–100 µm region. The green and purple lines shown in both the SED and image contours highlight two of the models that are poor fits to the observations. The purple model overpredicts both the flux in the SED and the surface brightness ratio between the top and bottom nebulae in the image.</figcaption></figure>

observations, from which they found a total disk mass (gas +dust assuming a gas-to-dust ratio of 100) of 0.057 ± 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(\text{gas} + \text{dust}) = \frac{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$, $\text{cm}^2\text{g}^{-1}$; 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 \text{cm}^2\text{g}^{-1}$ at 870 µm) and use $T = 20$ K, then that yields an estimated disk mass $M(\text{gas} + \text{dust}) = \frac{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
<figure><img src="image_13.png" /><figcaption>Figure 11. MCMC results using the covariance log-likelihood estimation. The blue crosshairs indicate the best-fit value for each parameter. Shading indicates the density of the parameter space sampling, while the red contours are drawn at the 1σ–4σ levels. All parameters are well-constrained, except for the surface density exponent (α). Dashed vertical lines represent the 16th, 50th, and 84th percentiles of the samples in the marginalized distributions.</figcaption></figure>

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.

### 3.8. Scale Height

The best-fit scale height of 16.2$^{+1.7}_{-2.0}$ 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):

$$H(r) = \sqrt{\frac{k_B T(r) r^3}{G M_{\text{star}} \mu m_p}}, \qquad (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 ($\lesssim 10$ µm) that dominate the opacity at visible wavelengths.

### 3.9. Flaring Exponent

The best-fit flaring exponent was $\beta = 1.19^{+0.09}_{-0.08}$. 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 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)(R/R_0)^{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 $\sim 1\sigma$ 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-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.

### 3.10. Surface Density Exponent

The surface density exponent is best fit by $\alpha = -1.77^{+0.94}_{-0.14}$, 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.

## 4. Discussion

### 4.1. Mass and Stability of the Disk

The best-fit dust mass ($0.00057 M_{\odot}$ or $190 M_{\oplus}$) and the associated total (gas + dust) disk mass ($60 M_{\text{Jup}}$, assuming a gas-to-dust ratio of 100) imply a disk mass to star mass ratio ($M_{\text{D}}/M_{\text{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 $\sim 50 M_{\text{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_{\text{Jup}}$ (implying a dust mass of $\sim 16 M_{\oplus}$).

This trend of low $M_{\text{D}}/M_{\text{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_{\odot}$), 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 ($\sim 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 $\ge 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_{\text{dust}} \ge 10^{-3} M_{\odot}$, which gives
$M_D/M_{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_{disk} \simeq 2 \times 10^{-2} M_{\odot}$, which gives $M_D/M_{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_{dust} = 1.0-1.75 \times 10^{-4} M_{\odot}$, which also gives $M_D/M_{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 $^{13}$CO line emission and find only $\sim 1.3 M_{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 $\sim 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 $\sim 16\%$ of the mass of the central star. We investigate the stability of the disk via the Toomre $Q$ parameter,

$$Q = \frac{c_s \kappa}{\pi G \Sigma}, \qquad (6)$$

where $c_s$ is the sound speed in the disk, $\kappa$ is the epicyclic frequency, and $\Sigma$ is the surface density profile of the disk. For a vertically isothermal disk with a Keplerian velocity, $\kappa = \Omega = \sqrt{GM_{star}/R^3}$ and a sound speed $c_s = \sqrt{k_B T/\mu m_p}$, where a reduced mass ($\mu$) 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.

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.

<figure><img src="image_14.png" /><figcaption>Figure 12. Radial profile of the Toomre $Q$ parameter for our best-fit disk. The disk appears to be stable at all radii.</figcaption></figure>

## 4.2. ESO-Hα 569 Compared to Other Cha I Disks

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_{\odot}$, which corresponds to 0.5% of the stellar mass, while the tail of the distribution stretches to 0.1 $M_{\odot}$ 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 under-luminous 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.

## 4.3. A Deficit of Edge-on 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 \pm 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
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Evolution of the shape of the image and SED for different dust masses. Our best-fit model is shown in blue. The other models use the same parameter values except for the mass, which is some fraction of the best-fit dust mass as indicated in the legend. For a fixed inclination, decreasing the dust mass moves photons from the thermal peak in the SED to the scattered-light peak. Decreasing the mass by a factor of 10 generates a flat SED without the double-peaked structure. Likewise, if the dust mass is one-tenth the best-fit value, the double-nebula shape begins to disappear in the scattered-light image and is not seen at all in the $1.9 \times 10^{-5} M_{\odot}$ model.</figcaption></figure>

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.

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.

## 5. Summary and Conclusions

We resolved the disk around ESO-Hα 569 in scattered light with HST/ACS and unambiguously confirmed that it is an optically thick protoplanetary disk viewed nearly edge-on. We performed radiative transfer modeling using a variety of fitting

techniques to constrain the geometry and grain properties of the disk. We successfully combine a covariance-based log-likelihood estimation with an MCMC framework to simultaneously fit the scattered-light image and literature-compiled SED for ESO-Hα 569. Our main results are as follows:

1. We find that a tapered-edge disk structure, with an exponential falloff of material outside of the apparent outer radius, is necessary to generate the diffuse scattered-light emission above the disk midplane, the flux ratio between the top/bottom nebulae of the disk, and the width of the dark lane simultaneously.

2. The best-fit disk mass of $0.057 M_{\odot}$ is abnormally large, especially considering the small central object, though multiple millimeter-continuum observations support this estimate. Assuming a gas-to-dust ratio of 100, the disk mass is 16% of the mass of the central star, establishing the ESO-Hα 569 disk as a clear outlier in the Cha I SFR. Despite the high mass, the disk appears gravitationally stable at all radii.

3. The vertical structure of the disk as defined by the scale height and the power-law flaring exponent is well-constrained. The best-fit model has a mass-averaged disk temperature of ~23 K, similar to other disk observations. The scale height is self-consistent with the modeled temperature profile, supporting a flared disk model in which the gas and dust are well-coupled.

A large effort was put into simultaneous and consistent fitting of the images and the SEDs, resulting in a disk model that is a good compromise between the two. But naturally, a
separate fit to each individual observable is capable of yielding a better fit to that one, at the cost of an inferior fit to the other. This is likely due to (1) limitations in the parameterization of the ongoing complex physical processes in protoplanetary disks. In this work, a fairly simple analytic disk structure formalism with a single grain population was used. (2) The inability of our data set to constrain some aspects of relevant physics and processes (e.g., neither the SED nor the scattered-light image provides much information on the innermost regions of the disk).

Using a combination of different observables (spectral data, images in scattered light, and thermal emission, and polariometry data to constrain grain properties) helps to break degeneracies between various model parameters. However, care must be taken to determine the correct approach for the relative weighting of observables with different noise properties and model sensitivities. Now that high-contrast imaging systems designed to study these circumstellar environments in greater detail are coming on line, there is a plethora of great observations for disks that formed under a range of initial conditions in a wide range of evolutionary stages. We may be entering an era where we have statistically significant numbers of circumstellar disk observations to employ population synthesis techniques. This is an important step if we hope to understand the inherent physics in the disk and planet formation processes. The tools we have been developing take us a step closer to being able to consistently make fits and measurements to e.g., the entire known sample of edge-on disks.

To better constrain the ESO-Hα 569 disk and stellar parameters, we would need to incorporate resolved images at multiple wavelengths. In this work, we chose not to model the 0.6 µm image because of the contamination from the jet. However, we recently obtained resolved images in the HST F475W filter and will use this to probe the diffuse scattering material high up above the disk in a forthcoming paper. Additionally, an ALMA Cycle 4 program (PI: F. Ménard) was awarded to map the thermal emission from 15 confirmed edge-on disks from our HST sample at 870 µm and 2 mm to probe dust settling, migration, and grain growth. Spatially resolved millimeter observations should go a long way toward disentangling many of the outstanding uncertainties regarding this disk’s structure. Looking forward, with the launch of JWST, the MIRI MRS integral field spectrograph will provide spatially and spectrally resolved data across the entire 5–30 µm range for many disks. This would not only help to constrain the structure of the disk in a regime where the current SED fit particularly struggles, but also provide valuable and detailed information about the dust species within the disk.

The authors thank K. Ward-Duong for helpful discussion. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under grant No. DGE-1232825. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program #12514. C.P. acknowledges funding from the European Commission’s 7th Framework Program (contract PERG06-GA-2009-256513) and from Agence

<table><caption>**Table 4**
Parameter values for the Synthetic Data Set</caption><thead><tr><td>Parameters</td><td>Values</td><td>Notes</td></tr></thead><tbody><tr><td>Inclination</td><td>71°6</td><td>Allowed to vary</td></tr><tr><td>Scale Height (R = 100 au)</td><td>25.6 au</td><td>Allowed to vary</td></tr><tr><td>Dust Mass</td><td>4.94 × 10<sup>-4</sup> M<sub>☉</sub></td><td>Held constant<sup>a</sup></td></tr><tr><td>Surface Density α</td><td>-1.76</td><td>Held constant</td></tr><tr><td>Flaring β</td><td>1.54</td><td>Held constant</td></tr></tbody></table>

**Note.**

a This value was held constant for all runs; however, a value of 4.94 × 10⁻⁵ M☉ (0.1 times the actual value) was used to test the robustness of the fitting techniques to systematic model errors.

Nationale pour la Recherche (ANR) of France under contract ANR-2010-JCJC-0504-01. F.M. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013.

*Software: MCFOST* (Pinte et al. 2006), Tiny Tim (Krist 1995), emcee (Foreman-Mackey et al. 2013), mcfost-python (https://github.com/swolff9/mcfost-python).

## Appendix

Here we provide a simplified disk model-fitting effort designed to illustrate the effect of the two “goodness-of-fit” metrics used in the MCMC explorations of parameter space: χ² and covariance log-likelihood-based estimation described in Sections 3.4.1 and 3.4.2, respectively. Although this demonstrates the power of the two tools, we recognize that it is not a comprehensive test of performance. A full benchmarking effort of the mcfost-python package is beyond the scope of this paper.

To test the ability of both fitting metrics, we generate an MCFOST model with known parameter values, add randomly generated 1σ noise to both the MCFOST-produced image and SED, and attempt to retrieve the parameters. The model was randomly drawn from the ESO-Hα 569 MCMC chain described above. We perform a fit to this synthetic data set using both the χ² and covariance log-likelihood-based estimation. For simplicity, we choose only to fit the scale height and inclination of our modeled disk. We expect both methods to recover the known parameter values within the uncertainties. Parameter values used for the synthetic data set are shown in Table 4.

To illustrate the power of the covariance framework over the χ²-fitting technique, we perform the same test, but purposefully input a disk dust mass too low by a factor of 10 into the MCFOST parameter file. This will test how robust the covariance framework in the presence of clear limitations in the model’s ability to fit the data.

In an effort to conduct these tests as close to the MCMC results reported above, we use the same Parallel Tempered ensemble sampler with two temperatures and 50 walkers. With only two free parameters, the chains converged more quickly, requiring only $N_{steps} = 10,000$ with $N_{burn} = 0.2 N_{steps}$. The allowable parameter ranges for the inclination and scale height were the same as reported in Table 2.

Figures 14 and 15 present the results for the covariance and χ²-fitting techniques, respectively. Both methods retrieve the input inclination and scale height within the uncertainties
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Left: MCMC results of the covariance log-likelihood estimation fit using a synthetic data set. We fit only the scale height and inclination of the modeled disk. The blue lines correspond to the known values for each parameter. The correct parameter values were retrieved, and the distributions are sharply peaked. Right: same as the left panel, but the MCMC run was conducted using an incorrect disk dust mass in the MCFOST parameter files. Even assuming a depleted dust mass, the scale height of the disk is still recovered, while the best-fit inclination is ~2σ discrepant. The covariance framework is less sensitive to any global limitations of the disk model to fit the given data set.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Same as Figure 15 but the MCMC was run using the $\chi^2$ log-likelihood-based estimation rather than the covariance framework. Unlike the covariance case, the $\chi^2$ fitting metric has a difficult time retrieving any of the correct parameter values when the incorrect disk dust mass was used to generate each MCFOST model. The disk scale height is not well-constrained at all, and the best-fit inclination is ~14σ discrepant.</figcaption></figure>

when using the correct dust mass. However, when the dust mass is set to one-tenth the actual value, both fitting methods struggle to retrieve the correct parameter values. The covariance run successfully recovered the disk scale height, though the uncertainties are larger than the correct dust mass case. The inclination was found to be $77.7_{-2.6}^{+5.1}$ degrees, which is only ~2σ discrepant from the true value. With the incorrect disk mass, the $\chi^2$ run was unable to recover either parameter. The scale height of the disk is not well-constrained at all, while the likelihood distribution for the inclination is sharply

peaked at $79.8_{-0.6}^{+1.1}$ degrees, which is ~14σ discrepant from the true value. It is unsurprising that the covariance framework is much more robust to global limitations of the models to fit the data set.

## ORCID iDs

Schuyler G. Wolff https://orcid.org/0000-0002-9977-8255
Marshall D. Perrin https://orcid.org/0000-0002-3191-8151
Gaspard Duchêne https://orcid.org/0000-0002-5092-6464
Francois Ménard https://orcid.org/0000-0002-1637-7393
Deborah Padgett https://orcid.org/0000-0001-5334-5107

Christophe Pinte https://orcid.org/0000-0001-5907-5179

William J. Fischer https://orcid.org/0000-0002-3747-2496

References

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Andrews, S. M., & Williams, J. P. 2007, ApJ, 659, 705

Bally, J., Walawender, J., Luhman, K. L., & Fazio, G. 2006, AJ, 132, 1923

Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Belloche, A., Schuller, F., Parise, B., et al. 2011, A&A, 527, A145

Bertout, C., Robichon, N., & Arenou, F. 1999, A&A, 352, 574

Boulanger, F., Bronfman, L., Dame, T. M., & Thaddeus, P. 1998, A&A, 332, 273

Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437

Cambresy, L., Epchtein, N., Copet, E., et al. 1997, A&A, 324, L5

Carmona, A., Pinte, C., Thi, W. F., et al. 2014, A&A, 567, A51

Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Chiang, E. I., Joung, M. K., Creech-Eakman, M. J., et al. 2001, ApJ, 547, 1077

Chiang, H.-F., Looney, L. W., & Tobin, J. J. 2012, ApJ, 756, 168

Cleeves, L. I., Oberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110

Comerón, F., Reipurth, B., Henry, A., & Fernández, M. 2004, A&A, 417, 583

Cotera, A. S., Whitney, B. A., Young, E., et al. 2001, ApJ, 556, 958

Cutri, R. M., et al. 2014, VizieR On-line Data Catalog, 2311, 0

Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2015, ApJ, 812, 128

David, T. J., Hillenbrand, L. A., Cody, A. M., Carpenter, J. M., & Howard, A. W. 2016, ApJ, 816, 21

Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503

Duchêne, G., McCabe, C., Pinte, C., et al. 2010, ApJ, 712, 112

Duchêne, G., Ménard, F., Stapelfeldt, K., & Duvert, G. 2003, A&A, 400, 559

Dunham, M. M., Offner, S. S. R., Pineda, J. E., et al. 2016, ApJ, 823, 160

Espaillat, C., Furlan, E., D'Alessio, P., et al. 2011, ApJ, 728, 49

Flaherty, K. M., Muzerolle, J., Rieke, G., et al. 2012, ApJ, 748, 71

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Glauser, A. M., Ménard, F., Pinte, C., et al. 2008, A&A, 485, 531

Goodman, J., & Weare, J. 2010, Comm. App. Math. Comp. Sci., 5, 65

Guidi, G., Tazzari, M., Testi, L., et al. 2016, A&A, 588, A112

Guilloteau, S., Dutrey, A., Piétu, V., & Boehler, Y. 2011, A&A, 529, A105

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Huélamo, N., Bouy, H., Pinte, C., et al. 2010, A&A, 523, A42

Hughes, A. M., Wilner, D. J., Qi, C., & Hogerheijde, M. R. 2008, ApJ, 678, 1119

Isella, A., Carpenter, J. M., & Sargent, A. I. 2010, ApJ, 714, 1746

Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714

Krist, J. 1995, in ASP Conf. Ser. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (San Francisco, CA: ASP), 349

Lebreton, J., Augereau, J.-C., Thi, W.-F., et al. 2012, A&A, 539, A17

Long, F., Herczeg, G. J., Pascucci, I., et al. 2017, arXiv:1706.03320

Luhman, K. L. 2007, ApJS, 173, 104

Luhman, K. L., Allen, L. E., Allen, P. R., et al. 2008, ApJ, 675, 1375

Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, arXiv:1704.02842

McCabe, C., Duchêne, G., Pinte, C., et al. 2011, ApJ, 727, 90

Millar-Blanchaer, M. A., Graham, J. R., Pueyo, L., et al. 2015, ApJ, 811, 18

Milli, J., Mawet, D., Pinte, C., et al. 2015, A&A, 577, A57

Muzerolle, J., Flaherty, K., Balog, Z., et al. 2009, ApJL, 704, L15

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Pinte, C., Padgett, D. L., Ménard, F., et al. 2008, A&A, 489, 633

Pohl, A., Sissa, E., Langlois, M., et al. 2017, A&A, 605, A34

Raftery, A., & Lewis, S. 1992, Bayesian Statistics, Vol. 4 (Oxford: Oxford Univ. Press), 763

Ribas, Á., Bouy, H., Merín, B., et al. 2016, MNRAS, 458, 1029

Ricci, L., Carpenter, J. M., Fu, B., et al. 2015, ApJ, 798, 124

Robberto, M., Spina, L., Da Rio, N., et al. 2012, AJ, 144, 83

Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256

Rodgers-Lee, D., Scholz, A., Natta, A., & Ray, T. 2014, MNRAS, 443, 1587

Sauter, J., Wolf, S., Launhardt, R., et al. 2009, A&A, 505, 1167

Sharma, S. 2017, arXiv:1706.01629

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Stapelfeldt, K. 2004, in IAU Symp. 202, Planetary Systems in the Universe, ed. A. Penny (Cambridge: Cambridge Univ. Press), 291

Stapelfeldt, K. R., Duchêne, G., Perrin, M., et al. 2014, in IAU Symp. 299, Exploring the Formation and Evolution of Planetary Systems, ed. M. Booth, B. C. Matthews, & J. R. Graham (Cambridge: Cambridge Univ. Press), 99

Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63

van der Plas, G., Ménard, F., Ward-Duong, K., et al. 2016, ApJ, 819, 102

Watson, A. M., Stapelfeldt, K. R., Wood, K., & Ménard, F. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 523

Weidenschilling, S. J. 1977, Ap&SS, 51, 153

Whitney, B. A., Indebetouw, R., Bjorkman, J. E., & Wood, K. 2004, ApJ, 617, 1177

Whitney, B. A., Wood, K., Bjorkman, J. E., & Cohen, M. 2003a, ApJ, 598, 1079

Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003b, ApJ, 591, 1049

Whittet, D. C. B., Prusti, T., Franco, G. A. P., et al. 1997, A&A, 327, 1194

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Winston, E., Cox, N. L. J., Prusti, T., et al. 2012, A&A, 545, A145

Woitke, P. 2015, European Physical Journal Web of Conferences, 102, 00007

Woitke, P., Pinte, C., Tilling, I., et al. 2010, MNRAS, 405, L26

Wolf, S., Padgett, D. L., & Stapelfeldt, K. R. 2003, ApJ, 588, 373

Wolff, S. G., Perrin, M., Ren, B., & Pinte, C. 2017, mcfost-python, v.1.0, Zenodo, doi:10.5281/zenodo.839863

Wood, K., Whitney, B. A., Robitaille, T., & Draine, B. T. 2008, ApJ, 688, 1118