A Resolved and Asymmetric Ring of PAHs within the Young Circumstellar Disk of IRS 48

In this work, we make use of new data sets acquired using the instrument NaCo (abbreviation for the Nasmyth Adaptive Optics System (NAOS) and the Near-infrared Imager and Spectrograph (CONICA)) commissioned on the Very Large Telescope (VLT) at the Paranal Observatory of the European Southern Observatory (ESO). This data was acquired in different near-infrared bands and at four epochs over a two-year period; these sets are detailed in Table 2. At that time, NaCo was integrated at the Nasmyth platform of the UT4 8.2 m telescope. We used this instrument in Sparse-Aperture Masking (SAM) and full-pupil imaging modes. For both modes, given the large reddening of the target (>16 mag in the visible), we use the NaCo built-in infrared wavefront sensor for the Adaptive Optics (AO) correction.

At epoch 1, full-pupil images were acquired in Ks-, Lp-, and Mp-bands, as well as spectra between 2.87 and 4.11 μm, with ${R}_{\mathrm{average}}\approx 800$. The Lp-band image is displayed in Figure 1 both with and without point-spread function (PSF) subtraction. One can easily see the east–west elongated nature of the circumstellar emission, at ≈100 milli-arcsecond (mas), as well as an asymmetry in the disk semimajor axis; the western extension is ≈1.15 times brighter than the eastern one. However, the ≈100 mas spatial scale is similar to the resolution criterion λ/D = 96 mas at these wavelengths, which enormously complexifies the use of these images.

Zoom In Zoom Out Reset image size

The slit of the spectrograph was orientated along the east–west direction as displayed in the latter figure. The pixel scale of the detector is constant over the whole spectra: 1.6 Å × 2.7 mas. The slit height is 86 mas; given the short scale extents of the disk and the high inclination (≈50°), we expect a large contamination of the disk at the stellar position (assumed to be located at the peak intensity in the spectra). Similarly, given the full width at half maximum of the PSF at ≈50 mas, we expect a large contamination of the star in the disk. To minimize the stellar contamination on the disk, we make use of the spectra at the first null of the star, 1.22λ/D ≈ 115 mas. The PSF-reference image shows a very low stellar flux in the null, ≈3% of the peak flux; this highlights the high strehl ratio of the observation.

The spectra in Figure 1 show the flux of the disk—taken at the first stellar null—relative to the peak flux of the spectra. The extraction aperture is three pixels for the peak, and one pixel interpolated at 1.22λ/D for the null. One can see very strong PAH emission features at 3.3 μm, characteristic of neutral PAH: in this emission line, the disk-to-peak contrast jumps by nearly a factor of two, compared to the disk-to-peak continuum contrast. Note that this emission line is not recorded in the interferometric data or full-pupil images given that it is outside the Lp-band filter bandpass, see Figure 1. Although both east and west spectra are very similar in shape, the Western disk-extension is constantly—at all spectrum wavelengths—1.15 times brighter than its eastern region. The same east–west asymmetry was measured on the full-pupil image.

The most interesting data sets are SAM data in that they offer an angular resolution down to a few tens of mas, equivalent to a few astronomical units given the distance of 120 pc, as shown in Table 3. Except for the longest baselines of ALMA (≳1 km), such resolution cannot be reached with any other instrument nor any other NaCo mode (coronographic, full-pupil, etc). Moreover, new extreme AO systems such as SPHERE usually have their wavefront sensors in the visible and cannot close the AO loop on this extremely extinct target (>16 mag in the visible), which means that NaCo-SAM is the only instrument able to observe IRS-48 at such resolution in visible or in the infrared (IR).

SAM, also called Non-Redundant Aperture Masking (NRM) is an observational technique in which a single mirror is transformed into a multi-pupil Fizeau interferometer. A mask placed in the pupil plane of the telescope transforms the usual Airy PSF into a pattern made of multiple fringes. The technique consists in measuring the amplitude and phase of these interferometric fringes so that information can be retrieved down to a fraction of a resolution element of the longest baseline of the mask (typically $0.5\lambda /{B}_{\max }\approx 50$ mas at 4 μm). However, the interferometric field of view remains very limited, in practice, of the order of $2.5\lambda /{B}_{\min }\approx 250$ mas at 4 μm (though a somewhat lower sensitivity exists until ≈6$\lambda /{B}_{\min }$). Using closure phases (T3)—an estimator robust to atmospheric turbulence and optical aberrations (Jennison 1958; Monnier 2000)—allows high-contrast capability without the inevitable "blind-spot" of other techniques such as coronography.

The observation strategy for the SAM mode was designed to intertwine the science target IRS-48 between two blocks on an interferometric calibrator, here Elia 2–35, 2–37, or 2–15, as detailed in Table 2. The mask used was a seven-hole non-redundant mask with baselines ranging from 1.78 to 6.43 m, allowing the simultaneous measurement of 21 visibility-squared (VIS2) and 35 closure phases (T3, among which 15 are independent). In full-pupil imaging mode, we used the same calibrators to have PSF references.

The data was processed using the SAMP pipeline (Lacour et al. 2011). It includes sky subtraction, bad pixel subtraction and fringe fitting. T3 and VIS2 measurements are calibrated by the two calibrators for each bracket.

Calibration to PSF reference help reduce systematic errors in the data, but hardly corrects all of them. A fitting mechanism to account for that is specifically developed, refer to Section 5. The error bars on the data are obtained from statistical measurements on the data, which account for random errors. T3 errors they average ≈05 in Lp-band and ≈17 in Ks  and M bands.

In addition to the multi-wavelength interferometric data listed above, we compiled SED fluxes for IRS-48 from the literature. All data are listed in Table 4 with the corresponding references. As can be seen, the photometry from the literature has a broad coverage from visible-blue to millimeter wavelengths. No measurements in the ultraviolet exist because of the very high extinction of the source; measurements in centimeter wavelengths are discarded given that they trace cold and large-scale structures outside the region of sensitivity and interest of our observations. For the fitting exercise, the errors on photometry fluxes were assumed to be of 5% if not specified in the literature.

The correction of the interstellar extinction was made using Cardelli et al. (1989), improved by O'Donnell's (1994) work in the visible. These extinction curves are parametrized with R_v (reddening slope) and A_v (absorption in V band) such that ${{\rm{A}}}_{v}\equiv {R}_{v}$*Extinction(B_mag–V_mag). The initial values are taken from Brown et al. (2012a), A_v = 11.5 and R_v = 5.5.

We take over the disk structure presented in Bruderer et al. (2014). Given the prominence of PAH in the SED of IRS-48, we merge into this disk model the results of Maaskant et al. (2014) on ionized PAH in several young disks. We obtain a disk structure in three parts: (1) an inner-most dusty disk between 0.4 and 1 au, (2) an outer dusty disk between 55 and 160 au, and (3) a PAH disk in between. For the sake of simplicity, we call this disk model the old model. Its main parameters can be seen in Table 5. This model represents the literature knowledge of IRS-48.

We use MCFOST, a parallelized 3D radiative transfer code based on a Monte-Carlo method (Pinte et al. 2006), able to output SED data and monochromatic images when given an input parameter file. After a few minor adjustments, the fit to the SED literature data is remarkable. However, this model or any small variation of it fails to explain the VIS2 data in either the Lp- or Mp-band. This highlights the degeneracy of SED-fitting exercises.

Indeed, the inner dust disk between 0.4 and 1 au emits very strongly in NIR, while the PAH grains further out are much dimmer. This creates a very high unresolved flux (inner disk and star) compared to the resolved flux from the PAH in the ≈9–55 au region. Hence, the inner regions of the old model disk appears spatially too small, mostly unresolved and the resulting visibilities are close to unity: much larger than the VIS2 data at all baselines. Given the difference of the VIS2 values, the resolved-flux depletion is estimated as ≈2.7 for baselines ≳3 m.

The only solution to make the old model fit the VIS2 is to incorporate a very high flux in the ≈9–35 au region—and in this region only—to make up for the contrast ratio of ≈2 in Lp-band between unresolved and resolved flux.

This implies fundamentally changing the structure, and hence the nature, of the disk around IRS-48.

In our data, VIS2 measurements are not excessively sensitive to asymmetries, essentially because of their relatively large error bars. Hence, they encode a point-symmetrical and smooth disk structure and inform on high-level morphology. We make use of the VIS2 data in Lp-band from all epochs and the Mp-band taken at epoch 1, as well as the SED fluxes from the literature (refer to Table 4). Ks-band from epoch 1 are discarded given that their large uncertainties would not set meaningful constraints on the disk structure. The Lp-band VIS2 data does not show time-variation along the four epochs.

We decided to perform a manual initial exploration of the parameters to identify the most impacting ones. We refine the fit later on with both a genetic algorithm and low-dimensional grids. The genetic algorithm allows for a very fast convergence for cases with a large number of dimensions where basic Markov chain Monte Carlo (MCMC) would fail. Grid calculations only come in a second time and ensure a full exploration of a small subset of the entire parameter space. Verifying smoothness for the outputs of the modeling over these small grid-subsets allowed us to interpolate the results and thus minimize the processing time for these grids. The final convergence was achieved by eye.

For the purpose of measuring visibilities on these MCFOST images, we choose their pixel resolution in "mas per pixel" to display the full interferometric field of view in a constant width of 1001 pixels. Also, in order to account for bandwidth smearing, we generate 15 monochromatic images linearly spaced in wavelength over the bandwidth of each filter (see Table 3) and combine them into a unique polychromatic image.

We perform an SED-fitting "inside-out," or "in increasing wavelengths," so that constraints and parameters obtained from fitting interior structures can remain (mostly) unchanged when fitting structures located further away from the star.

We assume an axial-symmetry perpendicular to the plane of the disk as well as plane symmetry above and below the plane of the disk. We adopt the standard flared disk prescription (Shakura & Sunyaev 1973), described in cylindrical coordinates (r, z) such that the density profile is given by

$$\begin{eqnarray}&&\sum (r,z)=\displaystyle \sum _{0}\,{r}^{p}\exp \left(-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{z}{h}\right)}^{2}\right),\end{eqnarray} \tag{ 1 }$$

where ${\sum }_{0}$ is a density-normalization constant, p is the surface density exponent, and h is the disk scale height; h varies with radius as $h={h}_{0}{(r/{r}_{0})}^{\beta }$, where β is the flaring exponent ($\beta \gt 0$, of order ≈1) and h₀ is the scale height at radius r₀.

The SED shape shows that IRS-48 consists of—at least—two main disk-components: an outer and extended disk beyond ≈55 au radius, and a closer one, that emits most of the PAH light, as found by Maaskant et al. (2014). Our VIS2 data does not constrain the outer disk very well, which mostly lies outside its field of view. Hence, we systematically adopt the simplest axisymmetrical fitting solution consistent with the literature for this outer disk, while we focus our fitting effort on the inner parts.

Similarly to the old model, we adopt a three-disk-component model: two dusty inner- and outer-disks and an additional disk composed of VSP in between. For each of them, the inner radius, outer radius, and total dust mass are free parameters. As an initial step, the parameters "scale height" and "flaring exponent" are shared between all three disks. Other less-impacting parameters, such as the surface density exponent, the grain minimum and maximum sizes, the ratio between carbonaceous (from Li & Greenberg 1997) and silicate (from Draine & Lee 1984) grains, the ratio between ionized and neutral PAH, and the grain-size power-law index or the disk inclination are left-aside for the initial fitting exercise and set to standard values from IRS-48 and disk literature. The flaring exponent = 0.67 (van der Marel et al. 2013), the inclination = 50° (Geers et al. 2007a; van der Marel et al. 2013; Bruderer et al. 2014), the surface density exponent = −1 (Andrews et al. 2013), the inner-rim radius = 55 au (Geers et al. 2007a), the silicate over carbonaceous grains mass ratio = 70/30 (prescriptions span 85/15 in Wood et al. 2002; 45/55 in Kim et al. 1994; 80/20 in Maaskant et al. 2014) and the min/max grain sizes = [0.03, 4000] μm (van der Marel et al. 2013) of the outer disk are set parameters, see Table 5.

The effective temperature of the central star is set to 9250 K, midway between Follette et al.'s (2015) work (9500 K) and Brown et al.'s (2012a) work (9000 K), who also estimated its spectral type to be A0${}_{-1}^{+4}$ (9000 ± 550 K), see Table 6.

VIS2 data gives stringent constraints on both the flux ratio between the unresolved flux (inner-most disk and star) and resolved flux (further than ≈9 au). To correct for the resolved-flux depletion of about a factor of 2.7 compared to the unresolved flux, one can either (1) decrease the inner-most disk opacity to decrease the unresolved flux or (2) increase the PAH disk opacity further than 9 au to increase the resolved flux. This later solution—increasing PAH emission in Lp band (3.8 μm) and Mp band (4.78 μm)—would immediately lead to a dramatic over-brightness of the main PAH emission features between 4 and 15 μm, with no (or few) degrees of freedom to compensate for it (i.e., structural parameters such as flaring exponent, scale height, etc).

Decreasing the disk brightness within the resolution angle of the instrument would automatically lead to an under-fitting of the NIR flux. This can easily be corrected by increasing the brightness of the star, de facto decreasing the amount of NIR excess arising from dust needed in the fit. Indeed, as discussed by Follette et al. (2015), a hotter and especially larger star questions the existence of a third inner-most disk in the ≈1 au area from the star where this disk was usually needed to account for the 1–3 μm IR excess. We hence decided to release fitting constraints by allowing the stellar luminosity to vary between 10 and 100 L_⊙. This implied adding two additional free parameters to tune the extinction of the target in order to fit the photometric fluxes measured in the visible and NIR (R_v taken in [4.0–10.0] and A_v taken in [10.0–15.0]).

A morphological fit of a simple two-dimensional Gaussian shows that the characteristic location of the resolved flux should lie in a region between 10 and 20 au. The PAH disk is a natural candidate to account for this flux, since previous work showed that the cavity up to ≈55 au was not depleted of such grains (Geers et al. 2007a; Maaskant et al. 2014). However, this means that the bright inner rim of the PAH disk must be fully resolved, i.e., further than ≈7–10 au, in order to account favorably for the contrast ratio between unresolved and resolved flux.

The MIR SED does not show any bump peeking at 9.5 and 18.3 μm, characteristic of silicate grains smaller than ≈2 μm. It does instead show a clear deep at 9.8 μm, de facto eliminating any major silicate grains contribution to the MIR flux; refer to Section 4.5 for a study on CTG in the inner regions. Hence, we decide not to incorporate any small or classical grains in the cavity between the inner  and outer disks.

To model the VSP, we follow the prescription of Natta et al. (1993) and use a mixture of carbonaceous VSG and PAH with an ionized/neutral fraction parameter, from Li & Draine (2001), Weingartner & Draine (2001), and Draine & Li (2007). Following the classification of Tielens (2008; Table 2), we set a_min = 4 Å and a_max = 10 Å for PAH (ionized and neutral), a_min = 10 Å and a_max = 30 Å for VSG (referred to as PAH clusters and VSG in their table). The grain size power-law index for PAH and VSG is set to −4; only the relative abundances of these three species are left as free parameters for fitting. Finally, we assume that these VSP are quantum heated by UV stellar-flux, i.e., outside local thermodynamic equilibrium. The quantum heating routines of MCFOST were benchmarked and validated by Camps et al. (2015).

Figure 4 shows the SED for both models with their respective interstellar extinction. Our new model provides the same fit quality as the old model. It matches very closely all the main features of the SED: the photosphere of the star, the NIR excess, the PAH emission features, and the outer disk. Most importantly, it dramatically improves the fit on the VIS2 data, which it explains entirely: Figure 3 shows the comparison between the VIS2 data in Ks, Lp, and Mp bands, and both models.

Zoom In Zoom Out Reset image size

The solution of a L_* = 42 L_⊙ star with a 3.7 × 10⁻¹⁰ M_⊙ mass VSP ring and no inner-most disk (0.4–1 au) is a complex and precise balance between the fit to the SED and to the VIS2 data. Indeed, the NIR SED requires either a brighter star than the old model, or a similar star with an inner-most disk to account for a large NIR excess. The VIS2 data, however, pinpoints a precise contrast ratio between the resolved VSP ring and the unresolved star and inner-most disk altogether. Adding even the smallest inner-most disk mass creates an excessive NIR emission, which must be compensated with a large increase of the VSP disk-mass to fit the VIS2 data. Such an increase, however, dramatically over-estimates the PAH emission features between 5.5 and 18 μm. This strongly pushes to a solution with no inner-most disk inside the VSP ring.

We follow the isochrones of Siess et al. (2000) young stellar object (YSO) evolutionary models with a luminosity increased from 14 to 42 L_⊙ and an effective temperature of 9250 K to find a larger stellar radius R_* = 2.5 R_⊙, rather than 1.8 or 1.4 R_⊙ as previously reported. The main stellar parameters are listed in Table 6, with a comparison to the previous literature values.

To compensate for the excess stellar luminosity, interstellar extinction was re-evaluated to A_v = 12.9 and R_v = 6.5. IRS-48 is 1.4 V-magnitude more extinct, with a more gray R_v than the extinction inferred in Brown et al. (2012a; A_v = 11.5 and R_v = 5.5). This confirms and reinforces the presence of substantial large dust grains in the line of sight (Kim et al. 1994; Indebetouw et al. 2005), typical of star-forming regions such as Ophiuchus. An equivalent value of A_k = 1.8 is consistent with that of ≈1.6 read on the low spatial-resolution extinction map computed by Lombardi et al. (2008); their estimate for that Ophiuchus-core region suffers, however, from a high discrepancy, most probably due to the unknown depth of their targets into this dense molecular cloud.

The radiative modeling shows that there is no extinction due to the circumstellar disk; its ≈50° inclination is not large enough to extenuate the stellar flux like (nearly) edge-on disks. All extinction is due to interstellar extinction, most of which arises from the local molecular cloud of ρ Oph.

The evolutionary tracks imply that such a star falls onto the ≈4 million year locus, with M_* = 2.5 M_⊙. This is much younger than the 8 or 15 Myr previously inferred, partly solving the evolutionary puzzle on the presence of a disk. A discussion on the age of IRS-48 is to be found in Section 6.1.

This new model consists of a smooth centro-symmetrical and optically thin VSP-only ring between 11 and 26 au. It is responsible for all of the IR-excess until ≈13 μm, where the outer disk starts dominating the dust-emission. In the fitting exercise, the VIS2 data constrained mostly the structure of the VSP ring, while the MIR-SED PAH emission features constrained the relative abundance of the PAH grains and VSG. The VSP density tops at ≈1 × 10⁻¹⁷ kg m⁻³ (see Figure 5).

Zoom In Zoom Out Reset image size

The resulting disk temperature is shown in Figure 6 for neutral PAH grains. Because of the close proximity to the star, the time between two successive absorptions is much lower than the relaxation time of the PAH grains. In this case, MCFOST treats quantum heating as "quasi-equlibrium" and a temperature can be defined (i.e., low dispersion on the temperature probability distribution). It reaches particularly high values for neutral PAH grains, T_max ≈ 1500 K and T_min ≈ 600 K. Ionized PAH grains and carbonaceous VSG follow a very similar temperature distribution, with lower temperatures due to their difference in electrical charge, mass, and structure: T_max ≈ 1100 K and T_min ≈ 400 K for ionized PAH and T_max ≈ 550 K and T_min ≈ 250 K for VSG.

Zoom In Zoom Out Reset image size

As previously highlighted, the interferometric data is mostly blind to the outer disk, and is not sensitive in any case to the southern asymmetry in millimeter wavelengths observed with ALMA. We fitted the simplest centro-symmetric disk model to satisfy the FIR and millimeter-wavelength SED.

This centro-symmetric approach is, however, not void of meaning. Indeed, the outer disk is optically thin in its emission window (for wavelengths ≳13 μm, τ in the mid-plane <10⁻¹, see Figure 7). Hence the precise azimuthal location of grains (azimuthally asymmetric or not) only plays a secondary role in fitting the SED.

Zoom In Zoom Out Reset image size

In order to satisfy the millimeter-wavelength SED, the outer radius of the disk is found at ≈250 au, similarly to Maaskant et al. (2014; 225 au), though SED fitting is not very sentivite to this parameter. For that disk extend, the total mass of the outer disk in the MCFOST model is 9 × 10⁻⁶ M_⊙ (3 M_Earth) for grains between 30 and 4 mm, with a grain-size power law of −3.5. This mass is similar, though lower than Bruderer et al. (2014; 1.6 × 10⁻⁶ M_⊙, 5.3 M_Earth) and 9 M_Earth found by van der Marel et al. (2013) for the same grain-size distribution and bounds. This is due to the fact that the later mass estimate obtained from 685 Ghz ALMA flux was biased as they used a lower temperature, T = 60 K at 60 au, than what is currently observed in our model. Indeed, our star is brighter, hence it heats the outer disk more efficiently: as seen in Figure 6, the inner-rim temperature peaks at T_dust = 148 K at 55 au.

We update their dust mass calculations with the values of our model to find ${\tau }_{685\mathrm{GHz}}=0.14$ for T = 140 K at 60 au, using their expression

$$\begin{eqnarray}&&{F}_{\nu }={\rm{\Omega }}\times {B}_{\nu }({T}_{\mathrm{dust}})\times (1-{\exp }^{-{\tau }_{\nu }}).\end{eqnarray} \tag{ 2 }$$

When we integrate over the whole disk model the optical thickness at ν = 685 GHz (λ = 440 μm), we find ${\tau }_{685\mathrm{GHz}}=0.19$, closely consistent with the previous 0.14 found independently using ALMA flux. Using ${M}_{\mathrm{dust}}=21\times {\tau }_{685\mathrm{GHz}}\sqrt{{a}_{\max }}$ from Draine (2006) and van der Marel et al. (2013) under the same hypothesis, we find M_dust = 3 M_Earth for τ = 0.14 and 4 M_Earth for τ = 0.19, consistent with the outer disk dust mass in our radiative transfer model (3 M_Earth).

Figure 8 shows model images generated in the mid-IR with the corresponding imaging data, at 4.78 μm (Mp-band), 8.6, and 18.7 μm. Although 8.6 and 18.7 μm images were not included in the fit—only SED fluxes were fitted at these wavelengths—our model can reproduce the observed structure accurately. Geers et al. (2007a) noted that the 8.6 μm image is mostly unresolved, with most flux arising from PAH in the inner 30 au, while the 18.7 μm image is highly resolved, with most flux arising from the inner rim of the outer disk located at ≈55 au.

Zoom In Zoom Out Reset image size

Figure 9 compares intensity profiles along the semimajor axis of the observed and modeled disks at 18.7 μm. Several models were computed, where the amount of thermal grains (silicates (70%) and carbonaceous (30%)) with sizes between 0.3 and 300 μm (mass power index of −3.5), are introduced in the VSP ring as settled material (one-fifth of the VSP-ring scale height). Although our axisymmetrical model cannot obviously reproduce the eastern asymmetry, it shows a distinctive dip at the star location corresponding closely to the 18.7 μm observed image. When the amount of classical thermal grains (CTG) increases, the inner-ring becomes brighter due to the efficiency of silicate grains to radiate at 18–19 μm. The central dip progressively disappears until the flux-profile is nearly flat, for a thermal grain mass of ≈10 M_VSP.

Zoom In Zoom Out Reset image size

Note that the VIS2 data fitting was not affected by the adding of CTGs, which emit a negligible amount of flux in the Ks, Lp, and Mp bands. Indeed, the emission of such CTGs is mostly visible through the emission "bumps" of the silicate grains, at ≈9.5 and 18 μm, where the SED is dominated by the outer disk. Until an addition of 5 M_VSP worth of CTGs in the VSP ring, the fit to the SED do not show any variation. From 10 M_VSP, however, the model starts to clearly overestimate the SED fluxes between ≈8 and 15 μm. Refer to Section 6.3 for a discussion.

We fit N point sources to the T3 data in order to highlight the main locations of the asymmetries in the disk. Such simplistic models allow analytical, and hence fast, fitting to the data.

We make use of six analytical models: from one to three point sources in addition to the central star, either with or without the disk modeled with MCFOST. The central star is set as a reference to a flux of 1 and in the middle of the field. Each point source is described with three parameters: (1) angular separation to the central star within [20, 300] mas, (2) position angle (east of north) within [0, 360]°, (3) relative flux within $\left[\tfrac{1}{200},1\right]$ stellar flux.

The disk is scaled and has its flux normalized with respect to that of the central star. It is treated as a "background image" on which T3 are measured at the (U, V) coordinates of the data. Complex visibilities from each unitary model object are added together and T3 values for the global model are obtained from the phases of the complex visibilities at the (U, V) coordinates of each T3 triangle.

The complex visibility ${ \mathcal V }$ at the (U, V) coordinates for N point sources is expressed by

$$\begin{eqnarray}&&\left(1+\displaystyle \sum _{i=1}^{N}{f}_{i}\right)\,{ \mathcal V }=1+\displaystyle \sum _{i=1}^{N}{f}_{i}\,\exp \left(-{\boldsymbol{j}}\displaystyle \frac{2\pi }{\lambda }\,(U\,{\delta }_{i}\,+\,V\,{\alpha }_{i})\right),\end{eqnarray} \tag{ 3 }$$

where λ is the central wavelength of the band. ${\delta }_{i}$, ${\alpha }_{i}$ are the angular separation and position angle of the ith point source with respect to the central star (projected along declination and right ascension axes), and f_i is the flux ratio.

The error estimation of interferometric data is extremely challenging, due to the intrinsic nature of the measurements that go through heavy data reduction. Systematic errors are calibrated using several PSF references (thanks to the observation pattern "Calibrator–Science Target–Calibrator"), while random errors are estimated from the scatter on the data. However, possibly large systematic errors remain even after calibration. The reason for this is a variation of experimental conditions between the measurements on the PSF reference(s) and the science target: either turbulence (i.e., AO correction), airmass, position on the detector, or position angle (i.e., angle-dependent aberrations in the optics).

A thorough treatment of the variance-covariance matrix of T3 requires at least as many independent measurements of each target as independent T3 (Ireland 2013), which is in practice not feasible. Instead, the variance-covariance matrix can be modeled (Kraus et al. 2008), or an additional error term can be added in quadrature (i.e., sum of independent variances) to the error estimation on the data and tuned to normalize the reduced chi-square distribution of the fitting (Hinkley et al. 2011). However, the nonlinear form of the complex visibility expression to evaluate, see Equation (3), makes it impossible to determine the number of degrees of freedom; the reduced chi-square cannot be used to compare models (Andrae et al. 2010). Rather than tuning it using reduced chi-square, we choose to fit the additional variance term ${\sigma }_{\mathrm{systematic}}^{2}$ according to a maximum-likelihood estimation. We define the total variance ${\sigma }_{i,\mathrm{total}}^{2}$ to be

$$\begin{eqnarray}&&{\sigma }_{i,\mathrm{total}}^{2}={\sigma }_{i,\mathrm{random}}^{2}+{\sigma }_{\mathrm{systematic}}^{2}\end{eqnarray} \tag{ 4 }$$

in order to account for error correlation and uncalibrated systematic errors of the T3.

We use an MCMC simulation using the emcee pure-Python implementation of Goodman & Weare's Affine Invariant Ensemble sampler Foreman-Mackey et al. (2013) to maximize the full log-likelihood expression $\mathrm{ln}{ \mathcal L }$ defined by

$$\begin{eqnarray}&&-2\,\mathrm{ln}{ \mathcal L }=\displaystyle \sum _{i}\left[\displaystyle \frac{{({\mathrm{Data}}_{i}-{\mathrm{Model}}_{i})}^{2}}{{\sigma }_{i,\mathrm{total}}^{2}}-\mathrm{ln}\left(\displaystyle \frac{2\pi }{{\sigma }_{i,\mathrm{total}}^{2}}\right)\right],\end{eqnarray} \tag{ 5 }$$

where the index i iterates over each (U, V) coordinate of a data set. While the first term under the sum favors high values of ${\sigma }_{i,\mathrm{total}}$ (hence higher values of the fitted additional variance ${\sigma }_{\mathrm{systematic}}$), the second term penalizes the increase of ${\sigma }_{\mathrm{systematic}}$. The balance between these terms ensures the convergence of the estimation.

In order to compare models, we calculate the Bayesian information criteria (BIC) at the maximum likelihood of each model. The BIC is defined by Schwarz (1978), such as

$$\begin{eqnarray}&&{\mathrm{BIC}}_{\mathrm{Model}}=2\,(k\mathrm{ln}(n)-\mathrm{ln}{{ \mathcal L }}_{\mathrm{Model}}{| }_{\max }),\end{eqnarray} \tag{ 6 }$$

where k is the number of free parameters in the model and n is the sample size. The presence of k weighted with the sample size ensures that the more simple models will be preferred (to a certain extent) over more complex models. We estimate the relative probability of all models with respect to the null hypothesis of a unique central star without any disk using likelihood ratio tests (LR-Test) in its most general formalism:

$$\begin{eqnarray}{\mathrm{LR}}_{\mathrm{Model}/\mathrm{null}} & = & \exp \left(\displaystyle \frac{{\mathrm{BIC}}_{\mathrm{null}}-{\mathrm{BIC}}_{\mathrm{Model}}}{2}\right)\\ & = & \displaystyle \frac{1}{{n}^{3N}}\displaystyle \frac{{{ \mathcal L }}_{\mathrm{Model}}{| }_{\max }}{{{ \mathcal L }}_{\mathrm{null}}},\end{eqnarray} \tag{ 7 }$$

where N is the number of point sources in the model.

We compute the relative probability for all models on each of the six data sets, which we tabulate as log of Equation (7) in Table 7. All six data sets except Lp-band epoch 3 favor the 2 point-source model rather than the 0, 1, or 3 point-source models, and all models favor the presence of a disk over a disk-less star (with the exception of Ks-band epoch 1, where the relative probabilities of a star with and without a disk are nearly identical). Since the fits on all T3 data sets were performed independently, there is a very clear statistical significance of a model that comprises a disk and two additional point sources. Note that given the relatively low flux from the point sources compared to that of the central star and the disk altogether, the additional two point sources do not significantly impact the fit of the radiative transfer disk model on the VIS2 data.

Finding the most probable model among several does not, however, ensure that the fit, hence the model, is satisfactory. We define residuals to be the normalized error of the model such that

$$\begin{eqnarray}&&{\mathrm{Residuals}}_{i}=\displaystyle \frac{{\mathrm{Data}}_{i}-{\mathrm{Model}}_{i}}{{\sigma }_{i,\mathrm{total}}},\end{eqnarray} \tag{ 8 }$$

where the index i iterates over each (U, V) coordinate of a data set. Table 8 shows that the additional systematic error ${\sigma }_{\mathrm{systematic}}$ in Ks and Lp bands is ≈08, and reaches 4° in the Mp band. These uncalibrated systematic errors are consistent with previous performance and slightly worse SPHERE-SAM mode: ≈05 (SPHERE SAM commissioning team, private communication). This table also shows that on all data sets, the rms of residuals calculated using ${\sigma }_{i,\mathrm{total}}$ is very close to unity, meaning this model does not under- or over-fit the T3 data. Note that the fact that rms values of residuals is unity is a pure consequence of the likelihood maximization of the additional variance parameter, not an a posteriori normalization. The average of residuals is only a fraction of unity, which denotes only a slight deviation from normal residuals (i.e., unbias).

Table 8.  Top Table: Overview of the Quality of the T3 Data Fit with the Favored Model (2 Point Sources and Disk), for all Data Sets

Data Set	$\langle {\sigma }_{i,\mathrm{random}}\rangle$ (°)	${\sigma }_{\mathrm{systematic}}$ (°)	$\langle {\sigma }_{i,\mathrm{total}}\rangle$ (°)	$\mathrm{rms}({\mathrm{Residuals}}_{i})$	$\langle {\mathrm{Residuals}}_{i}\rangle$
Ks epoch 1	1.6	0.86 ± 0.17	1.82	0.97	0.35
Lp epoch 1	0.27	0.68 ± 0.17	0.73	1.01	0.20
Mp epoch 1	1.9	4.0 ± 0.6	4.5	1.00	0.10
Lp epoch 2	0.44	0.63 ± 0.23	0.77	0.97	−0.21
Lp epoch 3	0.25	0.86 ± 0.11	0.89	1.00	0.18
Lp epoch 4	0.26	0.97 ± 0.11	1.0	1.00	0.26

	First Point-source			Second Point-source
Data Set	Sep. (mas)	PA (°)	Δmag	Sep. (mas)	PA (°)	Δmag
Ks epoch 1	101 ± 4	281 ± 4	4.9 ± 0.3	69 ± 5	16 ± 8	4.8 ± 0.3
Lp epoch 1	105 ± 2.5	270.3 ± 1.1	3.35 ± 0.12	${41}_{-3}^{+9}$	41 ± 4	2.7 ± 0.5
Mp epoch 1	89 ± 11	269 ± 4	2.9 ± 0.3	39 ± 9	52 ± 12	${1.8}_{-0.6}^{+1.0}$
Lp epoch 2	111 ± 2.5	275.1 ± 1.3	3.20 ± 0.20	${34}_{-2}^{+15}$	47 ± 4	${1.8}_{-1.0}^{+1.7}$
Lp epoch 3	109 ± 3	276.5 ± 1.3	3.40 ± 0.05	${39}_{-2}^{+10}$	44 ± 4	${2.3}_{-0.5}^{+0.7}$
Lp epoch 4	119 ± 4	284.6 ± 1.1	3.25 ± 0.06	50 ± 4	43.3 ± 1.1	2.72 ± 0.22

Note. The total resulting error ($\langle {\sigma }_{i,\mathrm{total}}\rangle$) defined by Equation (4), composed of the maximum-likelihood estimation of the additional error term (${\sigma }_{\mathrm{systematic}}$), and the average random errors of the data ($\langle {\sigma }_{i,\mathrm{random}}\rangle$). The rms of residuals remains close to unity, while its average remains close to zero, meaning no or little bias. Bottom table: maximum-likelihood estimation of the point-source parameters for the most probable model, consisting of two point sources in addition to a central star with the disk found in the radiative transfer VIS2 data fitting, for all data sets.

Download table as:  ASCIITypeset image

The results for the maximum-likelihood parameter estimation for the favored model with a disk and 2 point source are tabulated in Table 8. Although point sources were free to take any positions and fluxes, the parameter estimation locates both point sources at a consistent spot at all epochs and wavelengths. The first point source is ≈105 mas on the west of the star with a contrast ratio of ≈3.3 mag in the L band; the second point-source location is somewhat noisier due to its proximity to the star, though still consistently found at ≈45 mas on the northeast with a contrast ratio of ≈2.5 mag in the L band. These locations correspond closely to the inner rim of the VSP disk.

Fitting point sources to T3 data is equivalent to finding the photo-center of asymmetrically emitting regions; it does not constrain their size nor shape. Additional fits were carried using uniform disks (3 point-source parameters and an additional diameter), Gaussian blobs (3 point-source parameters and a semimajor axis) and two-dimensional Gaussian blobs (3 point-source parameters, 2 semimajor axes, and a position angle) instead of point sources. However, the data was not able to constrain well the additional parameters well: the point sources are mostly unresolved (i.e., angular extension ≲3 au).

We determine the color–temperature of both point sources at epoch 1, using their contrast ratios in Ks, Lp, and Mp bands. We assume a blackbody behavior constrained by its effective temperature and its radiating surface, assumed circular and defined by its radius. We use emcee to maximize the likelihood of the temperature parameter T taken in [100, 10,000] K and the logarithm of the radius parameter R taken in [0.1, 40] R_⊙.

The maximum-likelihood color temperature found for the first point source is ≈1350 K, see Figure 10. However, 1σ error bars yield several hundreds of Kelvin because of large error bars on the point-source contrast ratios obtained with T3 data. Still, this temperature corresponds closely to PAH-temperatures found with radiative transfer modeling at the unprojected point-source location, ≈12.5 au; refer to Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second point source shows a somewhat lower temperature, most likely of ≈1000 K. Here again, 1σ error bars yield several hundreds of Kelvin large: from 600 to 1300 K, refer to Figure 10. Due to its proximity to the star, at or below $0.5\lambda /D$, the angular separation of the second point source is particularly degenerate with its contrast ratio. However, the position angle remains mostly unaffected and shows no significant angular movement with epochs.

Discrepancy on the epoch 1 locations of the second point source between the different bands may be due to several factors, given that fitting T3 data only gives the location of the photo-center of an asymmetrically bright region. First, this location may change with wavelength due to the intrisic cause of the asymmetry: a clump of segregated grains, an arm of material showing decreasing temperatures as it spirals out from the star, or visual effects of optical depth, opacity, and shadows. Second, the data—especially the sample size and the error bars—allowed to fit only two point sources. This second point source may, however, well be the (poorly defined) photo-center of a more complex system of asymmetries, unresolved by the data.

The first point source in the Lp band shows a clear rotational movement through different epochs. It is unlikely related to the (U, V) field rotation, since it shows much smaller variations over the four epochs, and not monotonic. This movement is compatible with a circular orbit in the VSP disk and in the same direction as the gas kinematic found by Brown et al. (2012a) and Bruderer et al. (2014; with south as the near side). If we assume a circular orbit, the origin of the over-brightness is located at a radius of ≈14.2 au, and a height of ≈3.5 au (see below). This over-brightness location corresponds to the densest region of the VSP ring in our radiative transfer modeling, as Figure 5 tells. However, the nature of the over-brightness is not clearly constrained. It may be (1) an orbiting over-density of VSP material, or (2) the central star orbiting a center of mass, shifted from the center of the disk due to a massive companion.

Figure 11 shows the Lp-band point-source location in the disk as well as two circular-orbit solutions. The first solution is an orbit with a semimajor axis of a = 14.2 au. It lies flat 3.5 au above the mid-plane (hence inclination ≈0°). This suggests that either the clump of material producing the over-brightness is mostly pressure-supported or that only the top layer of it is detected. This latter option, however, implies that the center of mass and the photo-center of the clump are not co-located, inconsistent with the high transparency of this VSP ring (see Figure 7). The second solution is an inclined orbit, i = 14°, with a similar semimajor axis, a = 14.6 au; this implies a gravity-dominated clump of material, orbiting the star.

In both orbital solutions, the over-brightness moves at an unprojected Keplerian velocity of 12.7 ± 1.9 km s⁻¹ in average over the four epochs (similar values are obtained for other intervals, between epochs 1–3 and 3–4); the total orbital period at this pace is ≈34 years. The theoretical Keplerian velocity at such a distance and given the stellar mass is $\nu =\sqrt{G\times {m}_{* }/r}\,=12.5$ km s⁻¹.

Finding a brighter star, hence younger (4 Myr), partly solves an evolutionary puzzle on this target.

First, the discrepancy between the molecular Ophiuchus cloud age of ≈1 Myr (Luhman & Rieke 1999) and the stellar age is widely decreased, even though a fair difference remains that cannot be explained to date because this is the case for all class II disks in Ophiuchus. IRS-48 might have been part of an earlier star formation episode in Ophiuchus, as postulated by Brown et al. (2012a).

Second, the disk of our model is entirely dust-depleted up to 55 au, except for a thin VSP disk between 11 and 26 au. This implies that the evolutionary stage of the disk is probably more advanced than what was previously thought. Accretion was estimated at ≈10^−8.5 M_⊙ yr⁻¹ by Salyk et al. (2013) and Follette et al. (2015) using Paβ and Brγ lines. This non-negligible accretion rate raises the question of the origin of such material and does not seem consistent with either the old model (where one finds 8.0 × 10⁻¹² M_⊙ dust in the inner-most disk, and 8.0 × 10⁻¹⁰ M_⊙ PAH until 50 au), or our new model (3.7 × 10⁻¹⁰ M_⊙ VSP at 11–26 au). Moreover, the disk was found to be gas-depleted until 20 au, which puts a very hard constraint on how much material is to be found in the first tens of astronomical units, available for accretion. This means that some active replenishment processes take place to sustain accretion, or that IRS-48 is on the verge of becoming a debris disk.

Finally, a younger age for IRS-48 together with a later disk-evolution stage is consistent with the fact that massive stars erode their disk at a larger pace than solar-type stars. Indeed, earlier studies on Upper Scorpius OB association (Carpenter et al. 2006; Williams & Cieza 2011) showed that ≈20% of 127 K- and M-stars are surrounded by optically thick disks, while none of the 30 F- and G-stars showed such disks for a wavelength ≈16 μm, implying that primordial circumstellar disks are consumed faster by a massive star. This result is confirmed by Ribas et al. (2015) who show that only 2% of ${\rm{M}}\geqslant 2$ M_⊙ stars aged $\geqslant 3$ Myr show a protoplanetary disk.

6.2.1. 30 au CO Gas Ring

6.2.2.  VSP Ring Inner Truncation, Binarity

6.2.3. 26–55 au Cavity

Both the presence of the outer-rim of the VSP ring and the inner-rim of the outer disk require the existence of an explanatory phenomenon. Based on the millimeter asymmetry observed in the southern part of the disk, a 2 M_⊙ star and Hill radius estimates, van der Marel et al. (2013) suggested the presence of a planet at 17–20 au with M ≳ 10 M_Jup.

If we carry similar first-order Hill radius estimates with the new model disk structure and a higher-mass star, we speculate a M_p ≈ 3.5 M_Jup planet at a radius of a_p ≈ 40 au. Indeed, Dodson-Robinson & Salyk (2011) find that the outer edge of the gas cavity created by a planet is expected at ≈5 Hill radii r_H, defined by $3\times {({r}_{{\rm{H}}}/{a}_{p})}^{3}={M}_{p}/{M}_{* }$. Pinilla et al. (2012) find that for planets with masses between 1 and 3 M_Jup, dust dynamically accumulates at 7 r_H from their location. With the later planetary parameters, we are able to obtain the dust bump at 63 au (van der Marel et al. 2013), together with a VSP-ring bump at 17 au, and gas densities peaking at 26 and 56 au (all radii given as distances from the star). These radii correspond closely to the structure of our model; a schematic picture is shown on Figure 12.

Zoom In Zoom Out Reset image size

Dullemond et al. (2007) suggested that PAH and VSG do not settle well compared to CTG, and found that, under conditions of low turbulence, strong PAH emission features can "hide" a substantial amount of settled thermal silicate grains.

We were indeed able to add up to five times the mass of the VSP ring (5 M_VSP) worth of settled CTG in the VSP ring, without affecting the VIS2 data nor the SED fits. However, adding such grains in the close vicinity of the star affects the MIR emission in such a way that the 18.7 μm model image does not display the central dip with two bright spots at ±55 au from the star like the data image does (Figure 8). Instead, a contiguously bright region lies along the whole semimajor axis of the disk. This result shows that the VSP ring most probably does not contain more thermal grains than 5 × M_VSP ≈ 2 × 10⁻⁹ M_⊙.

The presence of PAH in many circumstellar disks reveals that in these objects the dust coagulation process was apparently not effective enough to remove the PAH from these disks, or that some other process continuously replenishes PAH grains (Dullemond et al. 2007). Out of theoretical and laboratory works, Jochims et al. (1994) showed that the largest and most regular PAH species are extremely stable to destructive radiations, i.e., photo-evaporation. Indeed, large PAHs of 20–30 carbon atoms will preferably relax through PAH emission in bands rather than photo-fragmentation. As a comparison, one of the most stable PAH species, circumcoronene C54H18, has a typical size of 4.87 Å. PAH used in our model span 4–10 Å (20–100 carbon atoms); VSG span 10–20 Å (100–1000 carbon atoms), see Tielens (2008).

Geers et al. (2006) find a PAH-to-dust mass fraction of 6% based on the abundance of 5 × 10⁻⁵ carbon atoms locked in PAH molecules per H nuclei (C_VSP/H). If we take updated values for PAH and VSG abundances from Tielens (2008; their Table 2), we find C_VSP/H = 3.5%. Our model showing a VSP-to-dust mass fraction >20% in the VSP ring between 11 and 26 au clearly highlights a depletion of a factor of ≈5–6 of classical dust grains up to 0.3 mm compared to VSP.

We are left with the problem of where one shall study the disappearance rate of VSP through grain growth, radiation pressure, or inward drift, versus the rate of some replenishment process, e.g., the collisional destruction of larger grains or inward drift from an outer reservoir.

There is presently little experimental data on the size distribution resulting from sub-micrometer grain–grain collisions. Tielens (2008) report that theoretical studies on grain–grain collisions result in a mass distribution with a power index of ≈−3.3. Although the minimum size is not constrained in this later study, given the layered structure of graphitic materials, it is expected that the smallest fragments are likely small, two-dimensional structures (e.g., PAH-like molecular species). Further processing by atomic reactions, as well as UV processing, may then quickly transform these species into compact PAH (Tielens 2008). In the frame of interstellar medium shocks, Jones et al. (1996) calculate that the catastrophic destruction (i.e., complete destruction of the grains, rather than cratering) of a 1000 Å grain by a 50 Å grain needs a critical velocity of approximately 75 km s⁻¹; larger grains need larger collisional energy. They also conclude that interstellar grain-growth processes shall be much more efficient than catastrophic destruction in order to explain the observation that most of the interstellar dust comes in grains larger than ≈1000 Å.

The orbiting velocity of the first point source in the VSP ring around IRS-48 is ≈12.7 km s⁻¹, similar to the Keplerian velocity at that distance for the 2.5 M_⊙ star. Moreover, A0 stars do not generally show strong winds or emit high-energy radiation. It seems very unlikely that IRS-48 possesses a high catastrophic collision rate within its VSP ring at 11–26 au. It is therefore expected that PAH and VSG cannot be replenished easily through grain–grain collisions in the VSP ring, in spite of a possible co-located reservoir of CTG up to five times the mass of VSP.

The radiation pressure estimate of the star can be achieved by calculating the stellar constant W ≈ 45,000 W m⁻² at 1 au. The radial acceleration from radiation pressure produced on a particle at radius r is obtained with a_p = W S_eff c⁻¹ r⁻², where S_eff is the effective cross-section of the particle and c is the speed of light in vacuum. It is found to be 39 (resp. 42, 17) times larger than the gravitational potential (G.M_*.r⁻²) on the neutral PAH particles (respectively, ionized PAH and VSG). This probably shows that the inward drift of such small particles is very low. The previous IRS-48 stellar model with radius 1.8 R_⊙ shows a somewhat lower but still overwhelmingly strong radiation pressure compared to gravitational potential for these species: 16, 17, and 7 times larger on average for neutral, ionized PAH and VSG.

The fact that the VSP ring was observed to be consistently similar at a radius between ≈11 and 26 au through the four epochs of observation spanning two years time, highlights that the VSP ring is governed by a much more complex balance of physical processes. In order to remain observationally in place, it must either be (1) constantly replenished by the outer-reservoir outside 55 au—possibly channeling triggered by the 40 au companion—or (2) kept in place as-is through (2a) gravitational shepherding from additional companions, or (2b) coating of large grains—unaffected by radiation presure—by PAH molecules.