© ESO, 2012

1. Introduction

The class of objects known as massive young stellar objects (MYSOs) are thought to represent massive stars at the earliest phases of evolution. The inherently large distances (several kpc) and high extinctions (A_V ≳ 20 mag) to typical MYSOs complicate attempts to characterize these objects observationally, and many studies (e.g. Grave & Kumar 2009) rely on fitting the spectral energy distribution (SED) to untangle the underlying physics involved in the formation of massive stars.

However, it has been shown that pure SED fitting, i.e. without making use of any spatial information, is a highly degenerate approach, at best (e.g. Thamm et al. 1994; Men’shchikov & Henning 1997). Different physical scales and temperature regimes of the circumstellar material can be resolved using a variety of techniques at different wavelengths, ranging from weakly-resolved single-telescope measurements, spectroscopic (spectroastrometric) measurements, and interferometric measurements, both at near-/mid-infrared and (sub-)millimeter wavelengths. The temperatures and regions thus probed range from  ~10³−10⁴ K gas within the dust sublimation radius using spectroastrometry of molecular gas lines (e.g. Wheelwright et al. 2010),  ~10³ K dust at the inner dust disk radius using near-infrared interferometry (e.g. Kraus et al. 2010),  ~10²−10³ K dust farther out in the disk (and possibly into outflow cones) using mid-infrared interferometry (e.g. Linz et al. 2009), and the broad distribution of cold  ~10¹−10² K dust using millimeter-wave interferometry (e.g. Beuther et al. 2007) and single-dish measurements (e.g. Beltrán et al. 2006).

In this paper, we present observations and the results of multi-wavelength modeling of one particular MYSO, AFGL 4176 (also known as IRAS 13395-6153, G308.9+0.1). This object was identified as a candidate MYSO based on the similarity of the far-infrared spectrum to that of the Becklin-Neugebauer (BN) object in Orion (Henning et al. 1984, 1990). The class of so-called BN-type objects are luminous (≳10³ L_⊙) sources, extremely red (usually with no optical counterpart), with a broad peak in emission (flat spectrum in νF_ν) between 50 and 200 μm and deep absorption due to silicates at 10 μm. They are generally interpreted as young (possibly still forming) massive stars, surrounded by a thick dust envelope, despite the unclear nature of the prototype object Orion BN (see the work by Tan 2004, for example).

At radio wavelengths, using the Molonglo Observatory Synthesis Telescope (MOST), Caswell et al. (1992) reported a small-diameter source (≤ 02) roughly coincident with the IR source AFGL 4176. The integrated flux densities are 0.31 Jy at 843 MHz and 0.40 Jy at 1415 MHz, confirming the emission is thermal bremsstrahlung radiation arising from an H II region. Using the Australian Telescope Compact Array (ATCA), Phillips et al. (1998) reported compact continuum emission at 8.6 GHz spanning about 10′′, with a collection of 6.7 GHz methanol masers on the northern edge. Later, Ellingsen et al. (2005) showed the 8.6 GHz continuum emission extends to about 1′. Finally, the 1.2 mm dust continuum was mapped by Beltrán et al. (2006), and comprises a roughly circular,  ~40′′ clump at the same location with an estimated mass of 1120 M_⊙ at 5.3 kpc, and a fainter clump (estimated mass 177 M_⊙ at 5.3 kpc) about 15 to the east.

In the near infrared, AFGL 4176 is a point source for seeing-limited observations. Using speckle interferometry, however, Leinert et al. (2001) found a faint extended halo in the K band, with a size of . The positions from the Two Micron All Sky Survey (2MASS) and the methanol masers reported by Phillips et al. (1998) are within 07. De Buizer (2003) detected a large number of shocked H₂ emission blobs around AFGL 4176 by means of narrow-band NIR imaging. However, their spatial distribution does not reveal a preferred direction for a potential jet, and in general, the association of AFGL 4176 with a jet or outflow remains unclear. The source was also observed in the mid- through far- infrared by the space missions IRAS, MSX, ISO, Spitzer and Herschel. We note that the coordinates listed in the IRAS Point Source Catalog appear to be erroneous when compared with the MSX 12 μm and Spitzer 24 μm images, and adopt the 2MASS position of (J2000) as the location of the infrared source. The ISO spectrum shows deep silicate absorption features at 9.7 and 18 μm, as well as absorption due to water ice at 3 μm, CO₂ ice at 4.3 μm, and emission lines of [Ne II] at 12.8 μm; [S III] at 12.0, 18.7 and 33.5 μm; [O I] at 145.5 μm and possibly 63.2 μm; and [C II] at 157.7 μm.

The distance to AFGL 4176 remains poorly constrained. Using both the velocity curve of Brand & Blitz (1993) with updated parameters from Levine et al. (2008), and the velocity curve of Reid et al. (2009), we find that the combination of the source’s LSR velocity of  ~−51 km s^-1 (e.g. Fontani et al. 2005; De Buizer et al. 2009) and galactic coordinates are forbidden¹. Consequently, a kinematic distance cannot be determined for AFGL 4176 using these rotation curves. Despite this, several authors (e.g. Beltrán et al. 2006; Fontani et al. 2005; Grave & Kumar 2009) report kinematic distances of around 5 kpc. Saito et al. (2001) showed that the spread in velocities of C¹⁸O molecular clumps in the vicinity is high (~30 km s^-1), meaning that kinematic estimates of the distance to this object should be highly uncertain, even after improvement of the galactic structure parameters. To date, other means of the distance measurements (e.g. maser parallax measurements) have not been applied to AFGL 4176. Because of this uncertainty, we consider the observational data in terms of both “near” (3.5 kpc) and “far” (5.3 kpc) distances, corresponding to the edges of the large C¹⁸O complex as determined by Saito et al. (2001).

The present paper is similar in scope to previous infrared interferometric studies of MYSOs by Linz et al. (2009); Follert et al. (2010); de Wit et al. (2010); Kraus et al. (2010). However, in this work on AFGL 4176, we go beyond traditional SED modeling techniques and use a far more comprehensive approach, including not only flux values, but also spatial information at mid-infrared, far-infrared, and sub-millimeter wavelengths directly into the fitting process. Such an approach helps to break the degeneracies inherent in modeling the SEDs of deeply-embedded sources, and provides structural information on scales ranging from tens of AU up to a parsec. We begin by presenting and discussing new observations of AFGL 4176, then go to one-dimensional radiative transfer modeling of the entire data set. Next, we use a geometric component analysis of the mid-infrared interferometric data to explore deviations from spherical symmetry, which we interpret as strong evidence for a circumstellar disk around this object. Finally, we compare the results of our study with those of previous studies of this object, as well as studies of other MYSOs.

2. Observations and data reduction

2.1. Mid-infrared interferometry

AFGL 4176 was observed with the two-telescope mid-infrared interferometer MIDI on the Very Large Telescope Interferometer (VLTI) of the European Southern Observatory on Cerro Paranal in Chile. Observations were conducted in 2005–2007 and 2011–2012 with both the 8.2 m Unit Telescopes (UTs) and the 1.8 m Auxiliary Telescopes (ATs) as part of guaranteed-time observations (proposal IDs 074.C-0389, 075.C-0755, 076.C-0757, 077.C-0440, 078.C-0712 and 084.C-1072). In total, 37 points in the uv plane were obtained for projected baselines ranging from 5 to 62 m and a variety of position angles, with six of the measurements obtained with the UT telescopes, while the remaining 31 were obtained with the AT telescopes. This represents, by far, the largest collection of long-baseline mid-infrared interferometric observations presented for a MYSO to date. A summary of the observations is shown in Table 1, and the uv coverage is shown in Fig. 1. For clarity, the individual observations have been color-coded by projected baseline length: purple indicates measurements with the longest baselines (60–64 m); red (36–49 m) and green (25–32 m) show those with intermediate baselines; blue (13–16 m) and orange (5–11 m) mark the shortest baseline measurements.

We used the prism as the dispersive element, resulting in spectrally-resolved visibility amplitudes and differential phases with a resolution of λ/Δλ = 35 in the wavelength range 8−13 μm. All observations were carried out in the HIGHSENS mode, meaning the correlated flux and the photometric measurements are done separately. For a detailed description of this observing procedure, see the work by Leinert et al. (2004). In general, every science measurement was preceded and/or followed by a calibration measurement of a bright single star of known diameter.

The raw data were reduced using the MIA+EWS package version 1.7.1. Masks were generated for each interferometric measurement from the corresponding photometric data. From the calibrator measurements taken throughout the night, we computed the mean sensitivity function, relating counts/s to Jy of correlated flux. The uncertainty in this sensitivity function was taken to be the standard deviation over all the calibrator measurements in the night, and is typically on the order of 10%.

Fig. 1 uv coverage of MIDI observations. North is up, east is left. The positions have been color–coded by projected baseline length, and the dashed lines indicate the mean projected baselines for each color (orange – 8.1 m; blue – 15.2 m; green – 28.5 m; red – 42.5 m; purple – 61.7 m).

Because the quality of the photometric data from the ATs is generally poor, we assume AFGL 4176 to be non-variable and use an average total spectrum from the UTs for computing the visibility amplitude from the correlated flux. Uncertainties in the final visibilities were calculated using standard error propagation formulae.

2.2. Far-infrared imaging

The AFGL 4176 region was covered by far-infrared observations taken by the Herschel Space Observatory (Pilbratt et al. 2010), within the framework of the Hi-GAL survey (Molinari et al. 2010). In this paper, we present bolometer data at 70 and 160 μm, obtained with the PACS instrument (Poglitsch et al. 2010) on August 16, 2010. These data are scan maps with a scan rate of 60′′/s, which slightly degrades the spatial resolution compared to the theoretical diffraction limit. After onboard preprocessing, the data are transmitted with a downlink rate of 5 frames/s for the 70 μm filter, and 10 frames/s for the 160 μm filter. The final resulting spatial resolution for the 70 μm data is about 92, and 126 for the 160 μm data, although the true point spread function (PSF) is asymmetric. For a more detailed description of the observing procedure used in the Hi-GAL program, we refer to the description of the observing strategy by Molinari et al. (2010).

We extracted the level-0 data from the Herschel Science Archive of the relevant 2° × 2° strip covering AFGL 4176 (OBSIDs 1342203085 and 1342203086). Despite the occurrence of very bright far-infrared sources in the Hi-GAL fields, these observations were performed with the normal high-gain settings for the PACS bolometers. The level-0 data were processed with the HIPE software (Ott 2010) version 8.0.3178 up to the so-called level-1, where internal data units have been translated to Janskys, and outliers (i.e. glitches) in the data and saturated pixels have been removed or masked.

In addition to these standard procedures, we applied a task for eliminating the electronic cross-talk within the bolometers, which can cause the appearance of ghosts close to strong compact sources. Furthermore, we use a module recently developed by the instrument team which applies a non-linearity correction to the data. This is especially important for our purpose, since in the high-gain mode used for the Hi-GAL program, very bright sources will trigger a non-linear response in the bolometer behavior, resulting in an altered PSF and decreasing measured fluxes by 10−15%.

The removal of 1/f noise and effects resulting from temperature drifts within the bolometers, as well as the final mapping, was done using the version 7 of the Scanamorphos program (Roussel 2012), for which we used the “galactic”, “noglitch” and “parallel” options. Final maps were produced with a pixel scale of 32 for both filters, which is the native plate scale for the 70 μm bolometers, and half the native plate scale for the 160 μm bolometers.

2.3. 870 μm imaging

AFGL 4176 was also imaged at 870 μm as part of the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL). We refer to the publication by Schuller et al. (2009) regarding the general outline of this sub-millimeter continuum survey, and the technical details of the observations and the data reduction. In the present work, we use a map from the re-reduction of the data in 2011, done by the Schuller, using refined methods for thermal drift removal within the bolometer reduction software BOA (Schuller et al. 2009). Furthermore, the treatment of strong sources within the iterative process of the data reduction has been improved. Finally, the data have been mapped onto a 6-arcsec grid, using version 2.18 of the mapping software SWarp (Bertin 2010), which prevents the formation of spurious Moiré patterns still present in earlier ATLASGAL maps.

3. Results

Fig. 2 Above: acquisition image of AFGL 4176 from MIDI obtained on the 8.2 m UT3 telescope using the SiC filter, which has a central wavelength of 11.8 μm and a width of 2.32 μm. Below: comparison of a cut through the acquisition image of AFGL 4176, and an analogous image of the unresolved calibrator star HD 102461.

We show an AO-corrected acquisition image of AFGL 4176 obtained with the 8.2 m telescope UT3 on Mar. 2, 2005 in the upper panel of Fig. 2. The image was obtained using the SiC filter, which has a central wavelength of 11.8 μm and a width of 2.32 μm. In the lower panel, we compare a cut through the acquisition image of AFGL 4176 with an analogous cut through the acquisition image of the unresolved calibrator star HD 102461, which was observed immediately after the science source. From the similarity of the spatial cut through the image of AFGL 4176 to that of the calibrator, we conclude that AFGL 4176 is at most marginally resolved in diffraction-limited 8 m-class telescope images at N-band wavelengths.

In Fig. 3, we show the average N-band spectrum of AFGL 4176 from the 8.2 m UT telescopes, used to derive the visibility amplitude from the correlated flux measurements. The visibility amplitudes as a function of wavelength are shown in Fig. 4, where the measurements have been ordered by projected baseline length, and are presented using the same color scheme as in Fig. 1. Finally, in Fig. 5, we show the differential phases measured with MIDI for the measurements with B > 36 m (the differential phases at shorter baselines are not shown, as they are essentially zero).

We note the clear presence of a zero crossing in both the visibility amplitude and differential phase for the two measurements with B = 61.6 and 61.8 m, indicated with a vertical dashed line in Figs. 4 and 5. We previously interpreted this as evidence for circular symmetry at the smallest scales probed (Boley et al. 2011). To test this hypothesis, we targeted this spatial frequency for follow-up observations with the ATs at other position angles on Feb. 24, 2012 (the remaining “purple” measurements in Figs. 1, 4 and 5). However, while a non-zero differential phase is measured, a zero crossing is clearly not present in these follow-up observations, implying deviations from spherical symmetry at the smallest angular scales probed by MIDI.

In Fig. 6, we show the images at 70 and 160 μm (Hi-GAL) and 870 μm (ATLASGAL). The approximate beam size for each image is 92, 126 and 192, respectively, and is shown as a circle in the lower left-hand corner of each image. The position of the infrared source is indicated as a cross in each image. The position of the emission peak in the Hi-GAL data corresponds to that of the 2MASS source to within the pointing accuracy of the observations. In the ATLASGAL data, the emission peak lies  ~4′′south of the IR position, although this is also within the rms pointing accuracy of the ATLASGAL observations.

Fig. 3 Average N-band spectrum measured with MIDI on the UTs.

At the Hi-GAL wavelengths, within a 36′′ aperture, we measure flux densities of 2382 Jy and 1383 Jy at 70 μm and 160 μm, respectively. While the absolute flux calibration on these data is accurate to within 5% (Herschel PACS Report PICC-ME-TN-037, Apr. 12, 2011), there is significant, non-uniform background emission in the region which is difficult to account for. Therefore, these values should be treated as upper limits, as it is not clear where the background begins to dominate and the extended emission associated with the source stops. The 870 μm emission is also quite extended; within a 19′′ aperture, which covers just the central emission peak, we measure an integrated flux density of 2.86 Jy, while in a larger 72′′ aperture, covering the extended emission, we measure an integrated flux density of 10.6 Jy.

Fig. 4 Complete set of visibility amplitudes measured with MIDI. The plots are sorted by projected baseline length, indicated as B in the inset text, together with the position angle θ. The color coding is the same as for Fig. 1.

Fig. 4 continued.

Fig. 5 Differential phases measured with MIDI for uv points with projected baselines longer than 36 m. The vertical dashed line for the two measurements with B = 61.6 and 61.8 m shows the location of the phase flip, indicative of a zero crossing in the visibility amplitude.

Fig. 6 The 70 and 160 μm images obtained from the Hi-GAL survey with the Herschel Space Observatory, and the 870 μm image from the ATLASGAL survey with the APEX telescope. The circle in the lower left corner of each image indicates the beam size. The two Hi-GAL images are shown with logarithmic scaling of the intensity, while the ATLASGAL image is shown with linear scaling.

In Table 2 we list the photometric data and effective aperture sizes from previous works and surveys which we consider in modeling our source. Unless noted otherwise, we compare our model data (see Sect. 4.1) to the observations by considering the flux contribution within a Gaussian aperture with the same FWHM as the observational data.

In Fig. 7, we show J-, H- and K_s-band images of AFGL 4176 and its immediate surroundings from the first data release of the ESO public survey VISTA Variables in the Vía Láctea (VVV; Saito et al. 2012). These images were obtained on the 4 m VISTA telescope on Cerro Paranal, and surpass 2MASS both in spatial resolution and sensitivity. The seeing for the J, H and K_s images presented here is 082, 086 and 094, respectively. AFGL 4176 is present as a bright, compact source, and saturated in both the H and K_s images.

Finally, the region around AFGL 4176 was mapped in several molecular lines as part of the Millimeter Astronomy Legacy Team 90 GHz (MALT90) survey (Foster et al. 2011) with the Mopra 22 m telescope. Of the 16 spectral lines which were mapped as part of the survey, HCO+ was the brightest in AFGL 4176. The HCO+ line data, shown in Fig. 8, have a spectral resolution of  ~0.11 km s^-1. This line is a tracer of dense molecular gas, and has been used as an outflow indicator (e.g. Hofner et al. 2001). However, in the case of AFGL 4176, the spectrum is very well approximated by a Gaussian (FWHM 3.5 km s^-1), and shows no signs indicating the presence of an outflow.

4. Analysis

4.1. One-dimensional radiative transfer modeling

As shown in Sect. 3, AFGL 4176 shows evidence for deviations from spherical symmetry only at the scales probed by our mid-infrared interferometric observations, which correspond to scales of roughly 20 to 200 mas. Furthermore, no outflow component has been revealed in searches by De Buizer (2003) and De Buizer et al. (2009), nor in the HCO+ data presented in this work (Fig. 8). Since there is no evidence for asymmetry or preferred geometry on scales larger than a few hundred AU, we begin with one-dimensional models to explore the overall structure of this source.

We model AFGL 4176 using the one-dimensional radiative transfer code MoDust, which was developed by de Koter and Bouwman. MoDust uses the Feautrier method to solve the radiative transfer equation (Feautrier 1964), and has been used previously by e.g. Malfait et al. (1999); Kemper et al. (2001); Bouwman et al. (2001). Using MoDust, the dust temperature structure is iteratively solved for a given density structure and dust properties. Different opacities can be given for different temperature regimes, thus making it possible to dynamically determine the distributions of different types of dust throughout the model.

Using this code, we produce radial intensity profiles of the source image on the sky at all wavelengths for which we have observations, i.e. from 1.25 μm to 1.2 mm. From these radial intensity profiles, we calculate aperture-corrected flux values for the SED, PSF-convolved images for comparison with the Hi-GAL and ATLASGAL data², and interferometric visibilities for comparison with the MIDI data. We do not attempt to match the near-IR VVV images, as AFGL 4176 is saturated in them. The relative goodness of fit of models was assessed by constructing a χ² estimate which includes roughly equal contribution from these three measures.

4.1.1. Constraints and uncertainties on fundamental parameters

As noted in the introduction, the distance to AFGL 4176 cannot be reliably determined at present. We therefore chose to investigate both near and far models, at distances of 3.5 and 5.3 kpc, respectively.

Under the assumption that the H II region seen in radio observations is powered by the AFGL 4176 source, the number of ionizing photons, and consequently a minimum inferred spectral type, can be determined. The radio flux measurements at 843 and 1415 MHz of Caswell et al. (1992) appear to be consistent with thermal bremsstrahlung radiation from an H II region. For a distance of 5.3 kpc, using the main-sequence luminosities for O stars from Martins et al. (2005) and the stellar atmosphere models of Kurucz (1979), we find a minimum spectral type of the central source of about O8. This corresponds to an effective temperature of about 33 000 K, which we adopt for the temperature of the central blackbody source in our models.

For the amount of foreground extinction, we adopt a simple prescription of 1 mag of visual extinction per kiloparsec, giving us A_V = 3.5 and 5.3 mag for the near and far models, respectively. This small amount of additional (i.e. interstellar) extinction was applied on top of our models, and has the effect of further lowering the near-IR flux levels and slightly deepening the 10 μm silicate feature.

4.1.2. Dust model

One of the foremost practical constraints on the choice of dust model is the need for opacities at short wavelengths (λ ≤ 0.1 μm). Since the luminosity of this source is on the order of 10⁵ L_⊙ and the amount of ionizing flux (inferred from the radio continuum observations) is equivalent to at least an O8 star star, the vast majority of the stellar radiation is emitted at UV and optical wavelengths.

A number of attractive dust models for conditions expected in protostellar cores, envelopes and protoplanetary disks have been developed (for example by Ossenkopf & Henning 1994; Henning & Stognienko 1996; Semenov et al. 2003). However, a common shortcoming of these models is the lack of opacity information at short (UV) wavelengths, due to convergence problems of the grain models (Semenov, priv. comm.). On the other hand, the empirically derived opacities of Laor & Draine (1993) and subsequent works cover wavelengths from 100 Å through 1 mm, with the caveats that the real physical makeup of the dust is unknown, and that at far-IR and (sub)mm wavelengths, where extinction is effectively impossible to measure, the dust opacities are not well known.

We therefore adopt a multi-component, temperature-dependent dust model, guided by the presence of known spectral features in the observations, and limited by the available coverage and reliability of different opacity models at different wavelength regimes. The transition between different regimes, and therefore the presence of different species throughout the surrounding material, are governed by the temperature structure iteratively solved by MoDust.

Based on the presence of the 3.1 μm water ice absorption feature in the ISO SWS spectrum, we infer that at least some of the amorphous silicate grains are surrounded by ice mantles, formed in the cold material of the parent molecular cloud. To include this in our model, we use opacities from the coagulation model of Ossenkopf & Henning (1994). Following the suggestions in that paper, we interpolate between Cols. 1 and 2 of their Table 1 using the depth of the 3.1 μm ice feature. From the ISO SWS spectra, we measure τ_3.1 = 0.4, and use the resulting interpolated opacities for regions where  K.

For grains with  K, we use the “astronomical silicates” and graphite grains of Laor & Draine (1993) with the “MRN” distribution of grain sizes (n(a) ∝ a^-3.5 for grains between 5 nm and 250 nm, where n is the number of grains, and a is the grain radius), which is generally appropriate for diffuse ISM conditions. We adopt a ratio of 60:40 for the silicate to graphite mass abundance in this temperature regime, although we note this is not directly constrained by observations. For the graphite grains, we adopt a destruction temperature of 950 K, which corresponds to the process of OH sputtering (Duschl et al. 1996). For the silicate grains, we do not enforce a destruction temperature; at the inner edge of our envelope models, the silicate grains typically reach a temperature of about 2000 K.

The dust model therefore consists of a region of “cold” (K) amorphous silicate and carbon grains with ice mantles, as well as “warm” (130 K K for graphites, and 130 K  for silicates) amorphous grains similar to the “standard” diffuse ISM mix.

4.1.3. Density structure

For the density structure of the dust in our model, we adopt a piecewise power law consisting of concentric, spherical shells, with no gaps in between the shells. In order to simultaneously reproduce the observed visibility levels, far-IR/sub-mm spatial distributions and SED with this model, we find it necessary to include at least three shells, with a jump in density at the interface between the first (innermost) shell and the second shell. At the interface between the second and third shells, however, we were able to achieve acceptable fits without a density jump, thereby avoiding the need for an additional model parameter describing this. We therefore characterize the density of material in the outer shells by the density ρ₂₃ at the contact point r₂₃ of the second and third shells. The precise structure used, therefore, has in principle nine free parameters (r₁, r₂, r₂₃, r_o, ρ₁, ρ₂₃, p₁, p₂, p₃) and follows the form (1)

4.1.4. Results

We present two fits to the observational data using the three-shell model: one for the “near” distance of 3.5 kpc, and one for the “far” distance of 5.3 kpc. The parameters are summarized in the first two lines of Table 3; comparisons with the observations are shown in Figs. 9 and 10: the top panel shows the modeled flux levels from the photometric data from Table 2; the center panel shows the modeled MIDI visibilities at three wavelengths; and the lower panel shows the radial distributions from Herschel (70 and 160 μm) and APEX (870 μm).

We defer a detailed discussion of the results to Sect. 5, but note the most prominent points of the model fits here. Both models reproduce the SED reasonably well at wavelengths longer than about 8 μm, but significantly underestimate the near-infrared flux levels. Both models also reproduce the overall visibility levels, as well as the spatial frequency of the zero crossing, however the far model provides a much better match, particularly at the two shorter wavelengths (red and green points in the middle panels of Figs. 9 and 10). Finally, despite the good agreement of the radial profiles at wavelengths of 70 and 870 μm, both models overestimate the extent of emission at 160 μm, although the near model is in better agreement with the observed radial profile. The 160 μm filter might be contaminated by strong [C II] emission at 157 μm, detected in emission in the ISO LWS spectrum, which could alter the spatial profile of the brightness distribution seen in this broadband filter. However, since no spatially-resolved [C II] spectral map is available for AFGL 4176, this ad-hoc explanation is not beyond doubt.

4.2. Two-dimensional geometric modeling

Fig. 7 J-, H- and Ks-band images from the VVV survey. AFGL 4176, indicated with a white cross, is saturated in both the H and Ks images. The circle in the lower left corner of each image indicates the seeing, as measured from the guide star by the autoguider.

While the one-dimensional approaches used above can provide useful insight, the deviations from spherical symmetry observed in the MIDI data should be explored. Such deviations can become especially clear in cases where multiple uv points have been measured with similar projected baselines (i.e., spatial frequencies), but differing position angles. Any significant variations in the visibility amplitude for a given spatial frequency will thus require deviations from a circularly-symmetric intensity distribution to be explained.

We see clear signs of such behavior for projected baselines longer than about 25 m, particularly for spatial frequencies larger than  ~10 fringe cycles/arcsec. For small variations in spatial frequency, we see visibility levels varying by over a factor of 5. Furthermore, there is a clear dependence on position angle: for measurements with V < 0.2, the visibility amplitude decreases with increasing position angle, while for measurements with V > 0.2, the visibility amplitude decreases with increasing position angle.

This behavior (visibility amplitude increasing with position angle for V < 0.2, and decreasing with position angle for V > 0.2) cannot be modeled with a simple elongated Gaussian. However, such a trend is possible in the case of a source distribution which has a well defined “edge” in its intensity distribution, causing a zero-crossing in the complex visibility with substantial power in the first side-lobe, for example a disk model.

Fig. 8 Spectrum of HCO+ from the MALT90 survey. The observations are shown with the gray line, while a Gaussian fit is shown with the black line.

Fig. 9 Model fits of the observational data for the near model. Top: the SED, where the aperture-corrected model flux levels are shown as crosses and the dotted line, and the observed flux values are shown as filled squares and the solid line. Middle: the modeled (lines) and observed (points) MIDI visibilities for three wavelengths. Bottom: the modeled (blue points) and observed (green points) radial intensity distributions from the Hi-GAL (70 and 160 μm) and ATLASGAL (870 μm) images.

Fig. 10 Model fits of the observational data for the far model. Figure notation is the same as for Fig. 9.

Fig. 11 Geometric fit of a disk to the MIDI visibilities at 10.6 μm. The color coding indicates the position angle of the measurements; circles show the observed values of the visibility amplitude, while triangles show the model values at the corresponding coordinates in uv space.

To model this asymmetric behavior, we use a simple parameterized disk to match the dependence on position angle seen in our MIDI visibility measurements. We approximate the disk with a geometrically thin, optically thick blackbody emitter, with a temperature depending on radius in the form of a power law in the form of (2)where T_f is the temperature at a fix-point r_f, and r_i and r_o are the inner and outer disk radii, respectively. We adopt T_f = 1500 K, i.e. the approximate evaporation temperature of silicate dust, and let , i.e. the radius at which a blackbody of luminosity L will have a temperature of 1500 K. For the luminosity of our far model, r_f = 37 AU. We find that models in which the inner disk radius is at r_f provide poor fits to our MIDI data, and a substantially smaller inner radius is required. We adopt r_i = 1 AU, which is well below the spatial frequencies probed by our MIDI measurements. We thus implicitly assume the inner disk, consisting of gas, remains optically thick. This is in line with interferometric measurements of high-mass young stars, whose disk emission appears consistently more compact than expected for dusty, disk-only models (e.g. Monnier et al. 2005).

Fig. 12 Image of the fit of the disk model to the mid-infrared visibilities at λ = 10.6 μm.

Fig. 13 Fit to model by Gürtler et al. (1991). Figure notation is the same as for Fig. 9, but in the top panel of the SED we distinguish between flux values available today (filled diamonds) and those available at the time of the Gürtler et al. (1991) study (filled squares).

Finally, to account for the more symmetric behavior at shorter baselines, we include a one-dimensional component in the image intensity distribution in the form of a Gaussian centered on the disk. Thus, as fit parameters, we use the outer radius of the disk r_o, the disk inclination θ and position angle φ, the exponent p of the temperature power-law, the flux ratio F_G/F_D of the Gaussian component to the disk component, and the FWHM θ_G of the Gaussian component.

4.2.1. Results

The parameters that best fit the observed MIDI visibilities were derived from a brute-force grid search of parameter space, where we further refine the parameters of the best fit from the grid search by using the downhill simplex algorithm. We find that this disk model produces a very good fit to our data, and in Table 4, we summarize the parameters of our best-fit model to the MIDI visibilities at a wavelength of 10.6 μm. We show the resulting model visibilities as triangles in Fig. 11, and the model image in Fig. 12.

The observations are best described by a disk spanning up to several hundred AU, with the disk temperature decreasing gradually as a function of radius (T ∝ r^-0.47). The disk is inclined (θ = 60°), and the semi-major axis is oriented 158° east of north. Besides the disk component, a substantial halo is also present, with a FWHM size equal to about half the diameter of the disk, and a total flux at λ = 10.6 μm of 1.6 times the disk component.

5. Discussion

5.1. Comparison with previous models

To gauge the importance of including spatial information in the modeling process, we reconstructed the one-dimensional radiative transfer model for this object presented by Gürtler et al. (1991), using their dust composition and density structure together with the MoDust code used for our models. The model from Gürtler et al. (1991), which we refer to as G91, is a two-zone power-law model, and can be described using the density prescription from Eq. (1) if the parameters r₁, ρ₁ and p₁ are omitted. In Table 3, we show the parameters for the G91 model using this prescription.

In Fig. 13 we show the resulting SED, mid-infrared visibilities, and far-infrared/sub-mm radial intensity profiles. Note that Gürtler et al. (1991) fit only the SED to the extent it was known at the time, and that they did not have detailed spatial information (i.e. mid-IR interferometry and far-IR imaging) available. We indicate the measurements they had available by squares in Fig. 13, while the more recent Hi-GAL and ATLASGAL measurements, as well as the 1.2 mm measurement from Beltrán et al. (2006), are indicated by diamonds. Regarding the flux values available at the time of their study, we note in particular that their 1.3 mm measurement of 131 mJy was probably in error, as the position indicated in their Table 2 is more than 50′′ from the infrared source. More recently, for example, Beltrán et al. (2006) reported a 1.2 mm flux value of 3.86 Jy for the source, which seems to be in good agreement with the ATLASGAL flux level of 11.56 Jy at 870 μm which we present here.

In Fig. 14, we show the density laws used for both our near and far models, as well as the G91 model. To account for the different distances to the source used in each model (see Table 3), we express the dust density as a function of angular position on the sky. We also overlay the effective angular scales probed by our MIDI measurements, where have used the same color coding and grouping by projected baseline as in Fig. 1.

For the near and far models, we see that the inner shell with enhanced density is probed primarily by the longest baseline measurements, while shorter baselines probe successively larger scales, out to about 01. Because this inner region is completely absent in the G91 model, the mid-IR visibilities at longer baselines are vastly under-predicted by this model. This is despite the relatively good fit to the SED (with the exception of the sub-mm/mm regime), as well as the 70 and 160 μm radial profiles, which only emphasizes the importance of including spatial information at different scales into such radiative transfer models.

5.2. The disk around AFGL 4176

The mid-infrared interferometric observations presented in this paper provide the first evidence for any signs of asymmetry in this deeply-embedded MYSO to date. Using a model consisting of a symmetric halo and an optically thick, geometrically thin disk, we find that a disk with an inclination angle of 60° and a position angle of the semi-major axis of 158° east of north provides an excellent fit to our MIDI data. The best-fit temperature profile of the disk is T ∝ r^-0.47. At longer wavelengths, the deviation from spherical symmetry in our MIDI data is much less pronounced, which suggests that the emission from the disk is dominated by a more spherical halo component at these wavelengths.

Nevertheless, because of the wavelength range, mid-infrared observations probe primarily cooler material (~200–600 K), expected to be well away from the hotter, inner circumstellar regions where signs of accretion onto the central star may be present. Furthermore, the inner radius of the disk around AFGL 4176, adopted in this work to be 1 AU, cannot be determined by our observations, as it remains unresolved by the MIDI data. In this context, near-infrared interferometric observations which probe these inner regions would make a superb compliment to the current data set.

On the other hand, the lack of detection of an outflow component, or any other larger-scale asymmetries, remains an open question. This may, for example, suggest that the object is still at a very early stage in the formation process, and that large-scale outflows have not had time to form. In any case, far-IR/sub-mm interferometric observations with ALMA may provide the key to bridging the symmetric structure on scales of tens to hundreds of thousands of AU to the asymmetries observed at tens to hundreds of AU.

5.3. Comparison with other MYSOs

Fig. 14 Density laws used for the three radiative transfer models, expressed as a function of angular distance on the sky to account for the difference in distance to the source each model. The colored rectangles show the spatial scales λ/2B covered by the MIDI observations, and the color coding matches the average baselines for the groups indicated in Fig. 1.

To date, very few MYSOs have been studied on scales of tens to hundreds of AUs. Here we compare our results with similar studies of five other MYSOs. In Table 5 we list the luminosities, the depth of the 10 μm silicate absorption feature, whether or not any compact H II region is present, and whether or not any outflows have been detected, which can be useful for independently constraining disk geometry. As previously noted, however, in the case of AFGL 4176, no outflow components have been detected (De Buizer 2003; De Buizer et al. 2009). For a broader discussion of these and other sources, we refer to the work by Grellmann et al. (2011).

The sources W33A and NGC 3603 IRS 9A have similar luminosities to AFGL 4176, on the order of 10⁵   L_⊙, and both have associated hyper-compact H II regions. However, IRS 9A seems to be much less embedded: besides having an optical component, the N-band spectrum shows no 10 μm silicate feature at all. At mid-infrared wavelengths, the source is very extended, and completely overresolved in MIDI observations probing  ~200 AU scales (Vehoff et al. 2010). W33A, on the other hand, could be even more embedded than AFGL 4176 (A_V ≳ 200 mag, de Wit et al. 2010), but the N-band visibility levels are poorly determined due to the low amount of (correlated) flux in the deep 10 μm silicate absorption feature. The source was modeled by de Wit et al. (2010), who concluded that most of the mid-infrared emission on scales of 10² AU is dominated by collimated outflow cones, without compelling evidence for significant emission from a dusty disk.

At an estimated luminosity of  ~5 × 10⁴ L_⊙, the SED of AFGL 2136 closely resembles that of AFGL 4176, including a similar depth of the 10 μm silicate feature (Smith et al. 2000). Using visibility measurements obtained with MIDI, both de Wit et al. (2011) and Follert (2011) reported mid-IR emission scales of  ~200 AU, while Follert (2011) found an elongated structure at these scales, parallel to earlier predictions of disk orientation.

Despite having a lower luminosity, on the order of 3 × 10⁴   L_⊙, the SED of IRAS13481-6124 is somewhat similar to that of AFGL 4176, although the 10 μm silicate feature is much shallower in IRAS13481-6124. Kraus et al. (2010) discovered a molecular outflow, and used near-IR interferometric observations to resolve an elongated structure oriented perpendicular to the outflow direction. The extended outflow shocks detected at near-IR wavelengths at large distances from the source (Stecklum et al. 2010³), as well as the radiative transfer models presented by Kraus et al. (2010), suggest the disk is strongly inclined. However, the 20 M_⊙ (gas + dust mass) disk at scales of 5–150 AU proposed by Kraus et al. (2010) is clearly a fundamentally different model than the structures we present here for AFGL 4176. For comparison, the dust mass enclosed in the near and far models for r < 150 AU is 8 × 10^-6 M_⊙ and 4 × 10^-5 M_⊙, respectively, although a direct comparison between studies using near- and mid-IR interferometric measurements is difficult.

Finally, the source with the lowest luminosity, M8E-IR, has an overall extinction of A_V ≈ 30 mag, as determined from the depth of the 10 μm silicate feature (Linz et al. 2009). Outflow signatures have been detected for this source (Mitchell et al. 1988). No cm continuum emission has been detected directly coinciding with the IR object. This lack of an ultra-compact H II region might be related to the hypothesis of a bloated, cooler central object, as suggested by Linz et al. (2009).

Given the small selection of sources for which this information is available, it is not possible to infer clear, empirical correlations between the properties listed in Table 5. Of the sources presented, none of them seem to be direct analogues of AFGL 4176, however W33A might be the closest match in terms of luminosity, SED and emission scales at mid-infrared wavelengths. Suggestions by de Wit et al. (2011) and Linz et al. (2011) indicate that mid-IR interferometric observations of at least some MYSOs may be heavily influenced by outflow components. On the other hand, the mid-IR visibilities presented by Follert (2011) apparently trace the disk structure of AFGL 2136. Taken together, it becomes clear that the mid-IR emission structure of MYSOs at scales of tens to hundreds of AUs can be complex, and caution should be exercised in interpreting data sets with limited uv coverage.

6. Summary and conclusions

In this paper, we reported on spatially resolved observations and modeling of the massive young stellar object AFGL 4176. We presented the largest amount of N-band visibilities for any MYSO to date, obtained with MIDI on the ESO/VLTI, together with maps of extended emission at 70 and 160 μm from the Hi-GAL survey on the Herschel Space Observatory and 870 μm imaging from the ATLASGAL survey with the APEX telescope.

With the exception of the mid-infrared interferometric measurements presented in this paper, the observational data available for this deeply-embedded object do not provide indications of asymmetry or preferred geometry. In view of this, we used spherically symmetric radiative transfer models, consisting of an envelope with a density enhancement at its inner edge, to model the source. Using these one-dimensional models, we were able to satisfactorily reproduce, simultaneously, the SED and spatial structure at mid-IR through mm wavelengths.

However, the one-dimensional models by definition cannot reproduce the asymmetries present in the MIDI data. We therefore modeled these data separately, using a multiple-component geometric analysis. We interpret the MIDI observations in terms of a circumstellar disk around AFGL 4176. If this disk hypothesis is correct, it would make AFGL 4176 the most luminous (massive) young star to have a disk detected around it to date. However, beyond initial detection, much remains to be clarified regarding the nature and extent of the disk, making this object a top priority for additional interferometric imaging, at both near-IR and far-IR/sub-mm wavelengths.

Acknowledgments

We would like to thank our referee, Dr. Koji Murakawa, for insightful critique during the review process of this manuscript. We thank Christoph Leinert, Svitlana Zhukovska, Dima Semenov and Keiichi Ohnaka for useful discussions. We also express our gratitude to Frederic Schuller for generously making the reduced ATLASGAL observations available to us, and Şeyma Çalişkan for supporting work during an earlier phase of this project. A.M.S. was partly supported by Russian state contract No. 16.518.11.7074 and the Russian Foundation for Basic Research (grants 11-02-01332 and 11-02-97124-p).

References

All Tables

All Figures

	Fig. 1 uv coverage of MIDI observations. North is up, east is left. The positions have been color–coded by projected baseline length, and the dashed lines indicate the mean projected baselines for each color (orange – 8.1 m; blue – 15.2 m; green – 28.5 m; red – 42.5 m; purple – 61.7 m).
In the text

	Fig. 2 Above: acquisition image of AFGL 4176 from MIDI obtained on the 8.2 m UT3 telescope using the SiC filter, which has a central wavelength of 11.8 μm and a width of 2.32 μm. Below: comparison of a cut through the acquisition image of AFGL 4176, and an analogous image of the unresolved calibrator star HD 102461.
In the text

	Fig. 3 Average N-band spectrum measured with MIDI on the UTs.
In the text

	Fig. 4 Complete set of visibility amplitudes measured with MIDI. The plots are sorted by projected baseline length, indicated as B in the inset text, together with the position angle θ. The color coding is the same as for Fig. 1.
In the text

	Fig. 4 continued.
In the text

	Fig. 5 Differential phases measured with MIDI for uv points with projected baselines longer than 36 m. The vertical dashed line for the two measurements with B = 61.6 and 61.8 m shows the location of the phase flip, indicative of a zero crossing in the visibility amplitude.
In the text

	Fig. 6 The 70 and 160 μm images obtained from the Hi-GAL survey with the Herschel Space Observatory, and the 870 μm image from the ATLASGAL survey with the APEX telescope. The circle in the lower left corner of each image indicates the beam size. The two Hi-GAL images are shown with logarithmic scaling of the intensity, while the ATLASGAL image is shown with linear scaling.
In the text

	Fig. 7 J-, H- and Ks-band images from the VVV survey. AFGL 4176, indicated with a white cross, is saturated in both the H and Ks images. The circle in the lower left corner of each image indicates the seeing, as measured from the guide star by the autoguider.
In the text

	Fig. 8 Spectrum of HCO+ from the MALT90 survey. The observations are shown with the gray line, while a Gaussian fit is shown with the black line.
In the text

Fig. 9 Model fits of the observational data for the near model. Top: the SED, where the aperture-corrected model flux levels are shown as crosses and the dotted line, and the observed flux values are shown as filled squares and the solid line. Middle: the modeled (lines) and observed (points) MIDI visibilities for three wavelengths. Bottom: the modeled (blue points) and observed (green points) radial intensity distributions from the Hi-GAL (70 and 160 μm) and ATLASGAL (870 μm) images.

In the text

	Fig. 10 Model fits of the observational data for the far model. Figure notation is the same as for Fig. 9.
In the text

	Fig. 11 Geometric fit of a disk to the MIDI visibilities at 10.6 μm. The color coding indicates the position angle of the measurements; circles show the observed values of the visibility amplitude, while triangles show the model values at the corresponding coordinates in uv space.
In the text

	Fig. 12 Image of the fit of the disk model to the mid-infrared visibilities at λ = 10.6 μm.
In the text

	Fig. 13 Fit to model by Gürtler et al. (1991). Figure notation is the same as for Fig. 9, but in the top panel of the SED we distinguish between flux values available today (filled diamonds) and those available at the time of the Gürtler et al. (1991) study (filled squares).
In the text

	Fig. 14 Density laws used for the three radiative transfer models, expressed as a function of angular distance on the sky to account for the difference in distance to the source each model. The colored rectangles show the spatial scales λ/2B covered by the MIDI observations, and the color coding matches the average baselines for the groups indicated in Fig. 1.
In the text