Improved Constraints on the Disk around MWC 349A from the 23 m LBTI

To estimate the size of the MWC 349A disk, we first fit uniform ellipses to the calibrated kernel phases and squared visibilities. This model is identical to that published in Danchi et al. (2001): a solid ellipse with semimajor axis R_out, position angle θ measured east of north, and axial ratio r. Depending on the disk inclination and geometry, a bright inner disk rim, gas or refractory dust within the sublimation radius, or the central star may be visible. We thus also fit geometric models that include central delta functions accounting for a fraction b of the total flux, beginning with a uniform ellipse plus delta function model. These two models are symmetric and cannot cause non-zero kernel phase measurements. Since the asymmetry in the 1.65 μm Keck image could have resulted from forward scattering from a flared disk, we also consider skewed ellipse models. The skewed ellipse is the uniform ellipse multiplied by a sinusoid in position angle, given by the following (e.g., Schaefer et al. 2010):

$$\begin{eqnarray}I=\left\{\begin{array}{ll}1+{A}_{s}\cos ({\phi }_{s}-\phi ), & {\rm{if}}\ {\left(\tfrac{{x}^{\prime }}{{R}_{\mathrm{out}}}\right)}^{2}+{\left(\tfrac{{y}^{\prime }}{{{rR}}_{\mathrm{out}}}\right)}^{2}\lt 1\\ 0, & {\rm{otherwise}}\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}{x}^{\prime } & = & x\cos (\theta )-y\sin (\theta )\\ {y}^{\prime } & = & x\sin (\theta )+y\cos (\theta )\\ \phi & = & \arctan \displaystyle \frac{y}{x}.\end{eqnarray} \tag{ 3 }$$

Here ${\phi }_{s}$ is the position angle at which the flux is brightest. Given the high temperature and luminosity estimates for MWC 349A, a clearing in the dust disk may be resolved. We thus also fit skewed ring plus delta function models to allow for a compact component (the star plus any gaseous/refractory material within the sublimation radius) and an outer disk. The skewed ring model is the skewed ellipse with an inner hole of radius R_in:

$$\begin{eqnarray}I=\left\{\begin{array}{ll}1+{A}_{s}\cos ({\phi }_{s}-\phi ), & {\rm{if}}\ {\left(\tfrac{{x}^{\prime }}{{R}_{\mathrm{out}}}\right)}^{2}+{\left(\tfrac{{y}^{\prime }}{{{rR}}_{\mathrm{out}}}\right)}^{2}\lt 1\\ 0, & {\rm{if}}\ {\left(\tfrac{{x}^{\prime }}{{R}_{\mathrm{in}}}\right)}^{2}+{\left(\tfrac{{y}^{\prime }}{{{rR}}_{\mathrm{in}}}\right)}^{2}\lt 1\\ 0, & {\rm{otherwise}}\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where ${x}^{\prime },{y}^{\prime },\phi ,{\rm{and}}\,{\phi }_{s}$ are identical to those in Equation (2).

We also fit two-dimensional Gaussians to the data to explore models without sharp edges. Like the solid ellipse fits, we first consider simple Gaussians and then add a central delta function and skew. The skewed Gaussian brightness profile is given by the following:

$$\begin{eqnarray}I & = & \,(1+{A}_{s}\cos ({\phi }_{s}-\phi ))\\ & & \times \exp \left[-4\mathrm{ln}2\left({\left(\displaystyle \frac{x^{\prime} }{r\mathrm{HWHM}}\right)}^{2}+{\left(\displaystyle \frac{y^{\prime} }{\mathrm{HWHM}}\right)}^{2}\right)\right],\end{eqnarray} \tag{ 5 }$$

where $x^{\prime}$, $y^{\prime}$, and ${\phi }_{s}$ are defined in the same way as Equation (2). We lastly fit Gaussian ellipses with inner clearings to the data. In order to make the simpler Gaussians a subset of these models, to make the ring model, we start with a simple non-skewed Gaussian. We then subtract a second Gaussian with identical position angle and axis ratio, but with HWHM scaled by f_HWHM. We constrain f_HWHM to be less than 1 to prevent a negative signal in the model images. We lastly apply skew and add a central delta function.

We fit the data using the Markov chain Monte Carlo algorithm emcee (Foreman-Mackey et al. 2013). We apply parallel tempering to ensure that the parameter space is fully explored in the case of multiple likelihood maxima. We calculate the 1σ parameter errors using the 16% and 84% contours from the chains at a temperature of one. To compare the various models, we calculate the Bayesian evidence (e.g., Trotta 2008), the integral of the posterior probability over the parameter priors, or the probability of a model given the data. The evidence ratios, or log evidence differences, between two models give their relative probabilities. Since Bayesian evidence is a noisy statistic with a non-zero false positive probability (e.g., Jenkins & Peacock 2011), we also compute ${\chi }^{2}$ differences to compare the models. For each model, we calculate the difference between its minimum ${\chi }^{2}$ value and that of the most complex model with fewer parameters.

We fit the data once, including the kernel phase and visibility error scalings as nuisance parameters. Since the best fits were nearly identical to the measured scatter, we present results where we fix the error scalings to the observed kernel phase and visibility scatter. We also perform fits to the intra-aperture baselines to understand how the full LBT resolution improves the model parameter constraints.

Table 1 lists the best-fit dual-aperture model parameters and Table 2 lists their corresponding minimum ${\chi }^{2}$ and Bayesian evidence values. The best-fit position angles agree for all models and are also consistent with the best-fit position angle reported in Danchi et al. (2001). For both types of brightness distributions, the Bayesian evidence and ${\chi }^{2}$ difference testing suggest that models including a compact component and skew are significantly better than the simpler models (see Table 2). These models provide a better match to the observations (see Figure 2).

Zoom In Zoom Out Reset image size

Table 1.  Geometric Model Fit Results

Ellipse Models
Model	Rout (mas)	θ (${}^{\circ }$)	r	As	${\phi }_{s}$ (deg)	b	Rin (mas)
Ellipse	$46\pm 2$	$97\pm 3$	0.65 ± 0.03	⋯	⋯	⋯	⋯
Ellipse + δ	$57\pm 2$	$99\pm 3$	0.66 ± 0.03	⋯	⋯	$0.30\pm 0.02$	⋯
Ellipse + δ + Skew	58 ± 2	$98\pm 3$	0.68 ± 0.03	0.17 ± 0.04	$-153{\pm }_{6}^{7}$	0.32 ± 0.01	⋯
Ring + δ + Skew	57 ± 2	98 ± 3	0.68 ± 0.03	0.16 ± 0.04	$-153{\pm }_{6}^{7}$	$0.33\pm 0.02$	$\lt 14$

Gaussian Models
Model	HWHM (mas)	θ (${}^{\circ }$)	r	As	${\phi }_{s}$ (deg)	b	fHWHM
Gaussian	28.2 ± 0.7	97 ± 3	0.64 ± 0.03	⋯	⋯	⋯	⋯
Gaussian + δ	34 ± 1	101 ± 3	0.64 ± 0.03	⋯	⋯	$0.23{\pm }_{0.03}^{0.02}$	⋯
Gaussian + δ + Skew	34 ± 1	101 ± 3	0.66 ± 0.03	0.24 ± 0.02	$-153{\pm }_{5}^{7}$	0.24 ± 0.02	⋯
Gaussian Ring + δ + Skew	$32{\pm }_{3}^{2}$	$101{\pm }_{3}^{4}$	$0.67{\pm }_{0.04}^{0.03}$	0.21 ± 0.06	$-153{\pm }_{8}^{6}$	0.4 ± 0.2	$0.31{\pm }_{0.04}^{0.03}$

Download table as:  ASCIITypeset image

Table 2.  Model Comparison

Model	${\chi }_{\min }^{2}$	dof	${\rm{\Delta }}{\chi }^{2}$ a	Δdofa	Significanceb	$\mathrm{log}Z$
Ellipse	264.5	239	⋯	⋯	⋯	−139 ± 2
Ellipse + δ	234.4	238	30.1	1	$5.5\sigma$	−126 ± 3
Ellipse + δ + Skew	213.9	236	20.5	2	$4.1\sigma$	−120 ± 3
Ring + δ + Skew	213.9	235	0.0	1	⋯	−121 ± 4
Gaussian	255.5	239	⋯	⋯	⋯	−134 ± 2
Gaussian + δ	228.3	238	27.2	1	$5.2\sigma$	−123 ± 3
Gaussian + δ + Skew	208.8	236	19.5	2	$4.0\sigma$	−117 ± 3
Gaussian Ring + δ + Skew	208.5	235	0.3	1	$\lt 1\sigma$	−118 ± 3

Notes.

^aWith respect to the above, simpler model. ^bDerived from the ${\chi }^{2}$ difference test.

Download table as:  ASCIITypeset image

Both the Bayesian evidence and the ${\chi }^{2}$ difference testing suggest that including an inner clearing does not improve the fit significantly. The Ring + δ + Skew model constrains any inner hole to have a radius less than 14 mas, but the best fit is indistinguishable from the Ellipse + δ + Skew model, given the resolution of the observations. While the Gaussian Ring + δ + Skew model has a slightly lower minimum ${\chi }^{2}$ than Gaussian models without an inner clearing, its ${\rm{\Delta }}{\chi }^{2}$ is low enough that it is not preferred at the $1\sigma$ level. It produces nearly identical observables to the Gaussian + δ + Skew best-fit model (see Figure 2). Its evidence value is also comparable to the Gaussian + δ + Skew best-fit model.

Figure 3 shows the posterior distributions for the Ring + δ + Skew model fit using both the intra- and dual-aperture observations. The 23 m LBTI places new and tighter constraints on all of the disk parameters compared to the single-aperture observations. The uniform ellipse model fit to the dual-aperture data results in comparable parameter errors as previous Keck studies (Danchi et al. 2001), but with $\sim 21 \%$ the number of squared visibilities and $\sim 6 \%$ the number of closure phases.

Zoom In Zoom Out Reset image size

We generate radiative transfer models to test whether a disk in radiative equilibrium with the central star can match the observations. We use the open source radiative transfer codes Hyperion (Robitaille 2011) and RADMC-3D (Dullemond 2012) and input the standard density profile for a flared disk:

$$\begin{eqnarray}&&\rho (r,z)={\rho }_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-\alpha }\exp \left(-\displaystyle \frac{1}{2}{\left[\displaystyle \frac{z}{h(r)}\right]}^{2}\right),\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&h(r)={h}_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{\beta }.\end{eqnarray} \tag{ 7 }$$

Here r and z are the radius and height in a cylindrical coordinate system. The radius value r₀ is where the scale height h is fixed to the constant value h₀ and the midplane density ρ is fixed to the constant value ${\rho }_{0}$. The density constant, ${\rho }_{0}$, can be found by integrating the density over all space with knowledge of the total disk mass. We first consider scale height (β) and density (α) power-law indices (1.25 and 2.25, respectively) consistent with irradiated disks in hydrostatic equilibrium (e.g., D'Alessio et al. 1998; Whitney et al. 2003). We set the disk inner radius at the point where the dust temperature reaches 1500 K to simulate dust sublimation, and use silicate dust with a grain size of $\sim 1\,\mu {\rm{m}}$ (e.g., Bans and Königl 2012). We show results with a disk mass of 0.01 ${M}_{* }$, but also explored ${10}^{-3}\,{M}_{* }$ and $0.1\,{M}_{* }$ disk masses and found that they produce comparable results. We vary the stellar temperature and luminosity within their estimated uncertainties (20,000–35,000 K for temperature and $3\times {10}^{4}\mbox{--}8\times {10}^{5}$ ${\text{}}{L}_{\odot }$ for luminosity.) We set the scale height to outer disk radius ratio at 0.01, and the disk inclination to 75° (Rodriguez & Bastian 1994). We also explore models with higher flaring indices, since MWC 349A may have a centrifugally driven disk wind (e.g., Martín-Pintado et al. 2011).

None of the radiative transfer models for passive irradiated disks match the observations. Figure 4 shows two example disk models for the upper and lower bounds on the temperature and luminosity for MWC 349A. For a low-luminosity MWC 349A, reprocessed light from the inner disk rim can account for the unresolved component in the geometric models. However, in this case, the outer regions of the disk are too cold to produce significant amounts of emission. A high-luminosity MWC 349A is bright on the correct scales along the disk major axis, but due to its inclination it cannot reproduce the visibilities along the minor axis. Asymmetric emission from the vertical wall at the disk inner edge also leads to a large phase signal.

Zoom In Zoom Out Reset image size

We reconstruct images using the BSMEM algorithm (Buscher 1994), assigning uniform closure phase and squared visibility errors of $6\buildrel{\circ}\over{.} 0$ and 0.08, respectively. Degeneracies exist between different reconstructed images from data sets with sparse $(u,v)$ coverage and small amounts of sky rotation. To illustrate this, we reconstruct images from simulated observations of the best-fit model images. We use the same $(u,v)$ coverage and sky rotation and add Gaussian noise at the level measured in the data. We then reconstruct images from both the data and the simulations using multiple priors.

Figure 5 shows images reconstructed from both the data and simulated observations of the best-fit skewed ring plus delta function model. Comparing the two rows of Figure 5 shows that the best-fit geometric model is consistent with images reconstructed using both priors. Comparing the two columns of Figure 5 shows that the reconstructed images depend on the choice of prior image. Additionally, degeneracies exist in the unresolved regions of the reconstructed images. The size and shape of the bright central component in each "Gaussian Prior" image is consistent with the size and shape of the synthesized beam. The fractional flux contained in the central component is roughly the same for the images made using each prior, and is approximately equal to the amount of flux contained within the synthesized main beam in the input model image. Putting a fraction b of the image flux into a central component will create identical closure phases and squared visibilities as long as the central component is unresolved. These degeneracies and the dependence on the prior image make reconstructed images ambiguous and necessitate model fitting in order to understand the source brightness distribution.

Zoom In Zoom Out Reset image size

The compact component in the geometric models accounts for $\lesssim 30 \%$ of the total image flux. Assuming a 3.78 μm flux of ∼100 Jy for MWC 349A (Thompson et al. 1977) yields a 30 Jy flux for the central component. Following Millan-Gabet et al. (2001), we can use the observed MWC 349A V- and L-band fluxes to estimate the amount of compact infrared excess. The emission expected for a star at temperature T with radius ${R}_{\star }$ is the Planck function times the solid angle, ${\rm{\Omega }}=\tfrac{\pi {R}_{\star }^{2}}{{d}^{2}}$, where d is the distance to the star. Using a dereddened V flux of 37.7 Jy and attributing it entirely to the star implies stellar radii of 13–28 R_⊙ depending on the chosen temperature and distance values. With this range of stellar solid angles and temperatures, the amount of unextincted stellar flux expected at 3.78 μm is then 1–3 Jy. Thus at least $\sim 90 \%$ of the compact flux is in excess, and this estimate increases if we include extinction when calculating the stellar flux at 3.78 μm.

Emission from a disk rim can account for the compact infrared excess if the stellar luminosity is low ($3\times {10}^{4}$ ${L}_{\odot }$) and the disk rim is close enough to the star to be unresolved. A higher luminosity ($8\times {10}^{5}$ ${L}_{\odot }$) MWC 349A sets the inner disk radius at ∼40 au in the absence of shielding (Figure 4). This is highly resolved and cannot contribute to a compact infrared excess. If the luminosity is indeed as high as $8\times {10}^{5}$ ${L}_{\odot }$, material such as optically thick gas (e.g., Eisner et al. 2009) or refractory dust (e.g., Benisty et al. 2010) must exist within the theoretical dust sublimation radius to explain the compact infrared excess. Thus the central geometric model component could be caused by a close-in inner disk rim, material within the dust sublimation radius, or some combination of the two.

The inferred stellar radius and compact infrared excess are consistent with both the YSO and B[e] supergiant scenarios for MWC 349A. Comparable stellar radii have been inferred for Herbig Ae/Be stars and B[e] supergiants (e.g., Zickgraf 2006; Fairlamb et al. 2015). Observations of B[e] supergiants suggest compact infrared excesses with comparable fractional flux to that for MWC 349A (e.g., Zickgraf et al. 1986; Kreplin et al. 2012). Large infrared excesses are found in observations of Herbig Ae/Be stars, in which the excess fractional flux can reach $95 \%$ (e.g., Millan-Gabet et al. 2001). Symmetric gaseous emission has been detected within the dust sublimation radius of several Herbig Ae/Be stars (e.g., Kraus et al. 2008; Tannirkulam et al. 2008; Eisner et al. 2009, 2010). This emission is $\lesssim$ au sized and consistent with the size of the compact component in the geometric models, given the distance estimates for MWC 349A. Gaseous emission coming from within the sublimation radius, which may be required if the stellar luminosity is high, would thus support an early age for MWC 349A.

The range of outer radii for the geometric disk models is 44–60 mas, corresponding to $53\mbox{--}102\,\mathrm{au}$ given the MWC 349A distance uncertainties. This is smaller than the gravitational radius for a photoevaporating disk, at which material would no longer be bound and could be lost in an outflow (Hollenbach et al. 1994). The gravitational radius can be written ${r}_{g}={{GM}}_{* }/{c}_{s}^{2}$, where G is the gravitational constant, ${M}_{* }$ is the stellar mass, and c_s is the sound speed. Assuming ${c}_{s}\,=11$ km s⁻¹ (Danchi et al. 2001), ${r}_{g}=219\mbox{--}290\,\mathrm{au}$ depending on the assumed stellar mass. The best-fit outer radii and position angles also agree with radio observations of the bipolar outflow and maser emission. The H30α maser emission spots are separated by 65 mas (Planesas et al. 1992) at a position angle of $107^\circ \pm 7^\circ$. The best-fit model is consistent with the width and orientation of the dark lane seen in VLA data as well (see Figure 6; Martin-Pintado et al. 1993). Thus the geometric model fits are consistent with a disk bound to MWC 349A at the center of the bipolar nebula and with the same orientation as the two maser spots.

Zoom In Zoom Out Reset image size

The best-fit ellipse size in Danchi et al. (2001) increases with wavelength from a major axis of 36 mas at 1.65 μm to 62 mas at 3.08 μm. These ellipse sizes, as well as the best-fit major axis presented here (88–120 mas at $3.78\,\mu {\rm{m}}$) follow a wavelength scaling close to ${\lambda }^{\tfrac{4}{3}}$. This trend is expected for flat, geometrically thin accretion disks as opposed to the ${\lambda }^{2}$ relation expected for flared disks (e.g., Malbet & Bertout 1995). Without complete radiative transfer models to compare to the data, Danchi et al. (2001) interpreted this as evidence for a flat disk around MWC 349A.

The radiative transfer simulations show that passive irradiated disks, which have the majority of their $3.78\,\mu {\rm{m}}$ flux near their inner rim, cannot match the observations given MWC 349A's inclination (75°). For low MWC 349A luminosity ($\sim 3\,\times {10}^{4}$ ${\text{}}{L}_{\odot }$), the extent of the emission is much too small to match the squared visibilities (Figure 4). For a higher stellar luminosity ($\sim 8\times {10}^{5}$ ${\text{}}{L}_{\odot }$), the asymmetric disk rim at larger angular separation causes a phase signal that is too large. A rounded inner disk wall would produce a lower phase signal (e.g., Monnier et al. 2006). This would be consistent with previous interferometric observations of Herbig Ae/Be stars, which could not be fit by models with simple vertical disk rims (e.g., Monnier et al. 2006; Millan-Gabet et al. 2016; Lazareff et al. 2017). However, even a perfectly symmetric ring (see Figure 7) does not match the data, since the large inclination shortens the appearance of the disk on the sky. This results in squared visibilities that fall off too quickly with baseline length. Rounded rim models with an inclination of $\sim 48^\circ$ can match the data; however, this is unlikely given previous constraints on the disk inclination from radio recombination line observations (e.g., Rodriguez & Bastian 1994).

Zoom In Zoom Out Reset image size

In both the high and low-luminosity cases, reproducing the observations requires additional emission, and thus heating, at large radii. The maser emission far from the star supports this scenario, since masers are often caused by shocks, which would heat nearby gas (e.g., Leurini et al. 2016). The presence of an ionized outflow (e.g., Martín-Pintado et al. 2011) is also consistent with heating at large radii. This extended emission may support a young age for MWC 349, since previous observations of Herbig Ae/Be stars suggest the presence of extended envelopes (e.g., Lazareff et al. 2017).

Some studies suggest that MWC 349 may be a hierarchical triple, where A is a close-separation binary surrounded by a circumbinary disk (e.g., Gvaramadze & Menten 2012). Regular brightness variations with a period of nine years (Jorgenson et al. 2000) suggest that MWC 349A may indeed be a close binary system with an orbital separation of ∼13 au (7.7–10.8 mas depending on the distance estimate to MWC 349A). Given their resolution, previous infrared interferometric observations cannot rule out an embedded binary with a separation $\lt 28$ mas (Danchi et al. 2001). Our observations also cannot rule out a close-separation binary morphology for MWC 349A. Model fits that include two point sources within the disk clearing can provide good fits to the data and do not tightly constrain the locations or fluxes of either inner component.